PEP 8 changes

wimprates/elastic_nr.py

@@ -7,18 +7,26 @@
 from scipy.integrate import quad
 
 import wimprates as wr
+
 export, __all__ = wr.exporter()
 
 
 @export
 def an(material='Xe'):
-    """Standard atomic weight of target (averaged across all isotopes)"""
+    """
+    Standard atomic weight of target (averaged across all isotopes)
+    :param material: The target material
+    :return: atomic weight of specifeikd
+    """
     if material is 'Xe':
         return 131.293
     if material is 'Ar':
-        return 39.948 
+        return 39.948
     if material is 'Ge':
         return 72.64
+    else:
+        raise NotImplementedError("unknown material %s" % str(material))
+
 
 @export
 def mn(material='Xe'):
@@ -29,8 +37,8 @@ def mn(material='Xe'):
 spin_isotopes = [
     # A, mass, J (nuclear spin), abundance
     # Data from Wikipedia (Jelle, 12 January 2018)
-    (129, 128.9047794, 1/2, 26.401e-2),
-    (131, 130.9050824, 3/2, 21.232e-2),
+    (129, 128.9047794, 1 / 2, 26.401e-2),
+    (131, 130.9050824, 3 / 2, 21.232e-2),
 ]
 
 # Load spin-dependent structure functions
@@ -65,13 +73,13 @@ def e_max(mw, v, m_nucleus=None):
     """
     if m_nucleus is None:
         m_nucleus = mn()
-    return 2 * reduced_mass(mw, m_nucleus)**2 * v**2 / m_nucleus
+    return 2 * reduced_mass(mw, m_nucleus) ** 2 * v ** 2 / m_nucleus
 
 
 @export
 def spherical_bessel_j1(x):
     """Spherical Bessel function j1 according to Wolfram Alpha"""
-    return np.sin(x)/x**2 + - np.cos(x)/x
+    return np.sin(x) / x ** 2 + - np.cos(x) / x
 
 
 @export
@@ -83,6 +91,7 @@ def helm_form_factor_squared(erec, anucl=None, material='Xe'):
 
     :param erec: nuclear recoil energy
     :param anucl: Nuclear mass number
+    :param material: name of the detection material (default is 'Xe')
     """
     if anucl is None:
         anucl = an(material)
@@ -94,31 +103,31 @@ def helm_form_factor_squared(erec, anucl=None, material='Xe'):
     #  and hardcoded constants...
 
     # First we get rn squared, in fm
-    mnucl = nu.amu/(nu.GeV/nu.c0**2)    # Mass of a nucleon, in GeV/c^2
+    mnucl = nu.amu / (nu.GeV / nu.c0 ** 2)  # Mass of a nucleon, in GeV/c^2
     pi = np.pi
-    c = 1.23*anucl**(1/3)-0.60
+    c = 1.23 * anucl ** (1 / 3) - 0.60
     a = 0.52
     s = 0.9
-    rn_sq = c**2 + (7.0/3.0) * pi**2 * a**2 - 5 * s**2
+    rn_sq = c ** 2 + (7.0 / 3.0) * pi ** 2 * a ** 2 - 5 * s ** 2
     rn = np.sqrt(rn_sq)  # units fm
     mass_kev = anucl * mnucl * 1e6
     hbarc_kevfm = 197327  # hbar * c in keV *fm (from Wolfram alpha)
 
     # E in units keV, rn in units fm, hbarc_kev units keV.fm
     # Formula is spherical bessel fn of Q=sqrt(E*2*Mn_keV)*rn
-    q = np.sqrt(en*2.*mass_kev)
-    qrn_over_hbarc = q*rn/hbarc_kevfm
+    q = np.sqrt(en * 2. * mass_kev)
+    qrn_over_hbarc = q * rn / hbarc_kevfm
     sph_bess = spherical_bessel_j1(qrn_over_hbarc)
-    retval = 9. * sph_bess * sph_bess / (qrn_over_hbarc*qrn_over_hbarc)
-    qs_over_hbarc = q*s/hbarc_kevfm
-    retval *= np.exp(-qs_over_hbarc*qs_over_hbarc)
+    retval = 9. * sph_bess * sph_bess / (qrn_over_hbarc * qrn_over_hbarc)
+    qs_over_hbarc = q * s / hbarc_kevfm
+    retval *= np.exp(-qs_over_hbarc * qs_over_hbarc)
     return retval
 
 
 @export
 def sigma_erec(erec, v, mw, sigma_nucleon,
-               interaction='SI', m_med=float('inf'), 
-               material = 'Xe'):
+               interaction='SI', m_med=float('inf'),
+               material='Xe'):
     """Differential elastic WIMP-nucleus cross section
     (dependent on recoil energy and wimp-earth speed v)
 
@@ -129,12 +138,13 @@ def sigma_erec(erec, v, mw, sigma_nucleon,
     :param interaction: string describing DM-nucleus interaction.
     See rate_wimps for options.
     :param m_med: Mediator mass. If not given, assumed much heavier than mw.
+    :param material: name of the detection material (default is 'Xe')
     """
     if interaction == 'SI':
         sigma_nucleus = (sigma_nucleon
                          * (mu_nucleus(mw, material) / reduced_mass(
-                            nu.amu, mw))**2
-                         * an(material)**2)
+                            nu.amu, mw)) ** 2
+                         * an(material) ** 2)
         result = (sigma_nucleus
                   / e_max(mw, v, mn(material))
                   * helm_form_factor_squared(erec,
@@ -143,7 +153,7 @@ def sigma_erec(erec, v, mw, sigma_nucleon,
 
     elif interaction.startswith('SD'):
         if material is not 'Xe':
-            raise not NotImplementedError("SI for %s-detector not available"%(
+            raise not NotImplementedError("SI for %s-detector not available" % (
                 material))
         _, coupling, s_assumption = interaction.split('_')
 
@@ -156,8 +166,8 @@ def sigma_erec(erec, v, mw, sigma_nucleon,
             # then divide by it in the next line.
             # Obviously there's no point to this, so let's not.
             x = (sigma_nucleon * 4 * np.pi
-                 * reduced_mass(mw, mn_isotope)**2
-                 / (3 * reduced_mass(mw, nu.mp)**2 * (2 * J + 1)))
+                 * reduced_mass(mw, mn_isotope) ** 2
+                 / (3 * reduced_mass(mw, nu.mp) ** 2 * (2 * J + 1)))
             result += (abundance
                        * x / e_max(mw, v, mn_isotope)
                        * s(erec / nu.keV))
@@ -173,17 +183,18 @@ def sigma_erec(erec, v, mw, sigma_nucleon,
 def mediator_factor(erec, m_med):
     if m_med == float('inf'):
         return 1
-    q = (2 * mn() * erec)**0.5
-    return m_med**4 / (m_med**2 + (q/nu.c0)**2)**2
+    q = (2 * mn() * erec) ** 0.5
+    return m_med ** 4 / (m_med ** 2 + (q / nu.c0) ** 2) ** 2
 
 
 @export
-def vmin_elastic(erec, mw, material = 'Xe'):
+def vmin_elastic(erec, mw, material='Xe'):
     """Minimum WIMP velocity that can produce a recoil of energy erec
     :param erec: recoil energy
     :param mw: Wimp mass
+    :param material: name of the detection material (default is 'Xe')
     """
-    return np.sqrt(mn(material) * erec / (2 * mu_nucleus(mw, material)**2))
+    return np.sqrt(mn(material) * erec / (2 * mu_nucleus(mw, material) ** 2))
 
 
 @export
@@ -204,8 +215,7 @@ def rate_elastic(erec, mw, sigma_nucleon, interaction='SI',
     If not given, conservative velocity distribution is used.
     :param halo_model: class (default to standard halo model)
     containing velocity distribution
-    :param progress_bar: if True, show a progress bar during evaluation
-    (if erec is an array)
+    :param material: name of the detection material (default is 'Xe')
 
     Further kwargs are passed to scipy.integrate.quad numeric integrator
     (e.g. error tolerance).