model renaming, putting vesc,v0 only in halo

tests/test_halo.py

@@ -1,11 +1,11 @@
 import pandas as pd
-from wimprates import j2000, standard_halo_model
+from wimprates import j2000, StandardHaloModel
 import numericalunits as nu
 import numpy as np
 
 
 def test_shm_values():
-    halo_model = standard_halo_model()
+    halo_model = StandardHaloModel()
     assert np.abs(halo_model.v_0 /(nu.km/nu.s) - 220.)<1e-6
     assert np.abs(halo_model.v_esc /(nu.km/nu.s) - 544.)<1e-6
 

tests/test_wimprates.py

@@ -23,7 +23,7 @@ def test_elastic():
     isclose(wr.rate_wimp_std(1, **opts), ref)
 
     # Test numericalunits.reset_units() does not affect results
-    nu.reset_units("SI")
+    nu.reset_units(123)
     isclose(wr.rate_wimp_std(1, **opts), ref)
 
     # Test vectorized call
@@ -63,5 +63,5 @@ def test_halo_scaling():
     #check that passing rho multiplies the rate correctly:
     ref = 33.19098343826968
 
-    isclose(wr.rate_wimp_std(1,halo_model = wr.standard_halo_model(rho_dm_value = wr.rho_dm()) , **opts), ref)
+    isclose(wr.rate_wimp_std(1,halo_model = wr.StandardHaloModel(rho_dm = 0.3 * nu.GeV/nu.c0**2 / nu.cm**3) , **opts), ref)
 

wimprates/bremsstrahlung.py

@@ -102,8 +102,8 @@ def integrand(erec):
 @export
 @wr.vectorize_first
 def rate_bremsstrahlung(w, mw, sigma_nucleon, interaction='SI',
-                        m_med=float('inf'), t=None, 
-                        halo_model = None,**kwargs):
+                        m_med=float('inf'), t=None,
+                        halo_model=None, **kwargs):
     """Differential rate per unit detector mass and recoil energy of
     Bremsstrahlung elastic WIMP-nucleus scattering.
 
@@ -113,7 +113,8 @@ def rate_bremsstrahlung(w, mw, sigma_nucleon, interaction='SI',
     :param m_med: Mediator mass. If not given, assumed very heavy.
     :param t: A J2000.0 timestamp. If not given,
     a conservative velocity distribution is used.
-    :param halo_model: class (default to standard halo model) containing velocity distribution
+    :param halo_model: class (default to standard halo model)
+    containing velocity distribution
     :param interaction: string describing DM-nucleus interaction.
     See sigma_erec for options
     :param progress_bar: if True, show a progress bar during evaluation
@@ -122,7 +123,7 @@ def rate_bremsstrahlung(w, mw, sigma_nucleon, interaction='SI',
     Further kwargs are passed to scipy.integrate.quad numeric integrator
     (e.g. error tolerance).
     """
-    halo_model = wr.standard_halo_model() if halo_model is None else halo_model
+    halo_model = wr.StandardHaloModel() if halo_model is None else halo_model
     vmin = vmin_w(w, mw)
 
     if vmin >= wr.v_max(t, halo_model.v_esc):

wimprates/elastic_nr.py

@@ -178,8 +178,8 @@ def vmin_elastic(erec, mw):
 @export
 @wr.vectorize_first
 def rate_elastic(erec, mw, sigma_nucleon, interaction='SI',
-        m_med=float('inf'), t=None, 
-        halo_model = None, **kwargs):
+                 m_med=float('inf'), t=None,
+                 halo_model=None, **kwargs):
     """Differential rate per unit detector mass and recoil energy of
     elastic WIMP scattering
 
@@ -191,7 +191,8 @@ def rate_elastic(erec, mw, sigma_nucleon, interaction='SI',
     :param m_med: Mediator mass. If not given, assumed very heavy.
     :param t: A J2000.0 timestamp.
     If not given, conservative velocity distribution is used.
-    :param halo_model: class (default to standard halo model) containing velocity distribution
+    :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)
 
@@ -200,14 +201,16 @@ def rate_elastic(erec, mw, sigma_nucleon, interaction='SI',
 
     Analytic expressions are known for this rate, but they are not used here.
     """
-    halo_model = wr.standard_halo_model() if halo_model is None else halo_model
+    halo_model = wr.StandardHaloModel() if halo_model is None else halo_model
     v_min = vmin_elastic(erec, mw)
 
     if v_min >= wr.v_max(t, halo_model.v_esc):
         return 0
 
     def integrand(v):
-        return (sigma_erec(erec, v, mw, sigma_nucleon, interaction, m_med) * v * halo_model.velocity_dist(v, t))
+        return (sigma_erec(erec, v, mw, sigma_nucleon,
+                           interaction, m_med) * v
+                * halo_model.velocity_dist(v, t))
 
     return halo_model.rho_dm / mw * (1 / mn()) * quad(
         integrand,

wimprates/electron.py

@@ -85,21 +85,27 @@ def v_min_dme(eb, erec, q, mw):
 
 
 # Precompute velocity integrals for t=None
-_v_mins = np.linspace(0, 1, 1000) * wr.v_max(None, wr.v_esc())
-_ims = np.array([
-    quad(lambda v: 1 / v * wr.observed_speed_dist(v),
-         _v_min,
-         wr.v_max(None, wr.v_esc() ))[0]
-    for _v_min in _v_mins])
-
-# Store interpolator in km/s rather than unit-dependent numbers
-# so we don't have to recalculate them when nu.reset_units() is called
-inverse_mean_speed_kms = interp1d(
-    _v_mins / (nu.km/nu.s),
-    _ims * (nu.km/nu.s),
-    # If we don't have 0 < v_min < v_max, we want to return 0
-    # so the integrand vanishes
-    fill_value=0, bounds_error=False)
+@export 
+def velocity_integral_without_time(halo_model=None):
+    halo_model = wr.StandardHaloModel() if halo_model is None else halo_model
+    _v_mins = np.linspace(0, 1, 1000) * wr.v_max(None, halo_model.v_esc)
+    _ims = np.array([
+        quad(lambda v: 1 / v * halo_model.velocity_dist(v,None),
+            _v_min,
+             wr.v_max(None, halo_model.v_esc ))[0]
+        for _v_min in _v_mins])
+
+    # Store interpolator in km/s rather than unit-dependent numbers
+    # so we don't have to recalculate them when nu.reset_units() is called
+    inverse_mean_speed_kms = interp1d(
+        _v_mins / (nu.km/nu.s),
+        _ims * (nu.km/nu.s),
+        # If we don't have 0 < v_min < v_max, we want to return 0
+        # so the integrand vanishes
+        fill_value=0, bounds_error=False)
+    return inverse_mean_speed_kms
+
+inverse_mean_speed_kms = velocity_integral_without_time()
 
 
 @export
@@ -118,7 +124,7 @@ def rate_dme(erec, n, l, mw, sigma_dme,
     If not given, a conservative velocity distribution is used.
     :param halo_model: class (default to standard halo model) containing velocity distribution
     """
-    halo_model = wr.standard_halo_model() if halo_model is None else halo_model
+    halo_model = wr.StandardHaloModel() if halo_model is None else halo_model
     shell = shell_str(n, l)
     eb = binding_es_for_dme(n, l)
 

wimprates/halo.py

@@ -10,26 +10,6 @@
 export, __all__ = wr.exporter()
 
 
-@export
-def rho_dm():
-    """Local dark matter density"""
-    return 0.3 * nu.GeV/nu.c0**2 / nu.cm**3
-
-
-@export
-def v_0():
-    """Most common velocity of WIMPs in the halo,
-    relative to galactic center (asymptotic)
-    """
-    return 220 * nu.km/nu.s
-
-
-@export
-def v_esc():
-    """Galactic escape velocity"""
-    return 544 * nu.km/nu.s
-
-
 # J2000.0 epoch conversion (converts datetime to days since epoch)
 # Zero of this convention is defined as 12h Terrestrial time on 1 January 2000
 # This is similar to UTC or GMT with negligible error (~1 minute).
@@ -110,17 +90,19 @@ def v_earth(t=None):
 
 
 @export
-def v_max(t=None,v_esc_value = None):
+def v_max(t=None, v_esc=None):
     """Return maximum observable dark matter velocity on Earth."""
-    v_esc_value = v_esc() if v_esc_value is None else v_esc_value #default args do not change value when you do a reset_unit so this is necessary to avoid errors
+    v_esc = 544 * nu.km/nu.s if v_esc is None else v_esc  # default
+    # args do not change value when you do a
+    # reset_unit so this is necessary to avoid errors
     if t is None:
-        return v_esc_value + v_earth(t)
+        return v_esc + v_earth(t)
     else:
-        return v_esc_value + np.sum(earth_velocity(t) ** 2) ** 0.5
+        return v_esc + np.sum(earth_velocity(t) ** 2) ** 0.5
 
 
 @export
-def observed_speed_dist(v, t=None, v_0_value=None, v_esc_value=None):
+def observed_speed_dist(v, t=None, v_0=None, v_esc=None):
     """Observed distribution of dark matter particle speeds on earth
     under the standard halo model.
 
@@ -131,58 +113,63 @@ def observed_speed_dist(v, t=None, v_0_value=None, v_esc_value=None):
     Optionally supply J2000.0 time t to take into account Earth's orbital
     velocity.
 
-    Further inputs: scale velocity v_0_value and escape velocity v_esc_value
+    Further inputs: scale velocity v_0 and escape velocity v_esc_value
     """
-    v_0_value = v_0() if v_0_value is None else v_0_value
-    v_esc_value = v_esc() if v_esc_value is None else v_esc_value
+    v_0 = 220 * nu.km/nu.s if v_0 is None else v_0
+    v_esc = 544 * nu.km/nu.s if v_esc is None else v_esc
     v_earth_t = v_earth(t)
 
     # Normalization constant, see Lewin&Smith appendix 1a
-    _w = v_esc_value/v_0_value
-    k = erf(_w) - 2/np.pi**0.5 * _w * np.exp(-_w**2) #unitless
+    _w = v_esc/v_0
+    k = erf(_w) - 2/np.pi**0.5 * _w * np.exp(-_w**2)  # unitless
 
     # Maximum cos(angle) for this velocity, otherwise v0
     xmax = np.minimum(1,
-                      (v_esc_value**2 - v_earth_t**2 - v**2)
+                      (v_esc**2 - v_earth_t**2 - v**2)
                       / (2 * v_earth_t * v))
-    #unitless
+    # unitless
 
-    y = (k * v / (np.pi**0.5 * v_0_value * v_earth_t)
-         * (np.exp(-((v-v_earth_t)/v_0_value)**2)
-            - np.exp(-1/v_0_value**2 * (v**2 + v_earth_t**2
-                                    + 2 * v * v_earth_t * xmax))))
-            #units / (velocity)
+    y = (k * v / (np.pi**0.5 * v_0 * v_earth_t)
+         * (np.exp(-((v-v_earth_t)/v_0)**2)
+         - np.exp(-1/v_0**2 * (v**2 + v_earth_t**2
+                  + 2 * v * v_earth_t * xmax))))
+    # units / (velocity)
 
     # Zero if v > v_max
     try:
         len(v)
     except TypeError:
         # Scalar argument
-        if v > v_max(t, v_esc_value):
+        if v > v_max(t, v_esc):
             return 0
         else:
             return y
     else:
         # Array argument
-        y[v > v_max(t, v_esc_value)] = 0
+        y[v > v_max(t, v_esc)] = 0
         return y
 
-@export 
-class standard_halo_model:
+
+@export
+class StandardHaloModel:
     """
         class used to pass a halo model to the rate computation
-        must contain: 
+        must contain:
         :param v_esc -- escape velocity
-        :function velocity_dist -- function taking v,t giving normalised valocity distribution in earth rest-frame. 
+        :function velocity_dist -- function taking v,t
+        giving normalised valocity distribution in earth rest-frame.
         :param rho_dm -- density in mass/volume of dark matter at the Earth
         The standard halo model also allows variation of v_0
-        :param v_0 
+        :param v_0
     """
-    def __init__(self, v_0_value =None, v_esc_value=None,rho_dm_value=None):
-        self.v_0 = v_0() if v_0_value is None else v_0_value
-        self.v_esc = v_esc() if v_esc_value is None else v_esc_value
-        self.rho_dm = rho_dm() if rho_dm_value is None else rho_dm_value
-    def velocity_dist(self,v,t):
-        #in units of per velocity, 
-        #v is in units of velocity 
+
+    def __init__(self, v_0=None, v_esc=None, rho_dm=None):
+        self.v_0 = 220 * nu.km/nu.s if v_0 is None else v_0
+        self.v_esc = 544 * nu.km/nu.s if v_esc is None else v_esc
+        self.rho_dm = 0.3 * nu.GeV/nu.c0**2 / nu.cm**3 if rho_dm is None else rho_dm
+
+    def velocity_dist(self, v, t):
+        # in units of per velocity,
+        # v is in units of velocity
         return observed_speed_dist(v, t, self.v_0, self.v_esc)
+

wimprates/migdal.py

@@ -41,7 +41,7 @@ def vmin_migdal(w, erec, mw):
 @wr.vectorize_first
 def rate_migdal(w, mw, sigma_nucleon, interaction='SI', m_med=float('inf'),
                 include_approx_nr=False,
-                t=None, halo_model = None, **kwargs):
+                t=None, halo_model=None, **kwargs):
     """Differential rate per unit detector mass and deposited ER energy of
     Migdal effect WIMP-nucleus scattering
 
@@ -58,14 +58,15 @@ def rate_migdal(w, mw, sigma_nucleon, interaction='SI', m_med=float('inf'),
         presented the Migdal spectra.
     :param t: A J2000.0 timestamp.
     If not given, conservative velocity distribution is used.
-    :param halo_model: class (default to standard halo model) containing velocity distribution
+    :param halo_model: class (default to standard halo model)
+    containing velocity distribution
     :param progress_bar: if True, show a progress bar during evaluation
     (if w is an array)
 
     Further kwargs are passed to scipy.integrate.quad numeric integrator
     (e.g. error tolerance).
     """
-    halo_model = wr.standard_halo_model() if halo_model is None else halo_model
+    halo_model = wr.StandardHaloModel() if halo_model is None else halo_model
     include_approx_nr = 1 if include_approx_nr else 0
 
     result = 0

wimprates/summary.py

@@ -10,7 +10,7 @@
 @export
 def rate_wimp(es, mw, sigma_nucleon, interaction='SI',
               detection_mechanism='elastic_nr', m_med=float('inf'),
-              t=None, halo_model = None ,
+              t=None, halo_model=None,
               **kwargs):
     """Differential rate per unit time, unit detector mass
     and unit recoil energy of WIMP-nucleus scattering.
@@ -43,7 +43,7 @@ def rate_wimp(es, mw, sigma_nucleon, interaction='SI',
     Further kwargs are passed to scipy.integrate.quad numeric integrator
     (e.g. error tolerance).
     """
-    halo_model = wr.standard_halo_model() if halo_model is None else halo_model
+    halo_model = wr.StandardHaloModel() if halo_model is None else halo_model
     dmechs = dict(elastic_nr=wr.rate_elastic,
                   bremsstrahlung=wr.rate_bremsstrahlung,
                   migdal=wr.rate_migdal)
@@ -57,23 +57,24 @@ def rate_wimp(es, mw, sigma_nucleon, interaction='SI',
 
 @export
 def rate_wimp_std(es, mw, sigma_nucleon, m_med=float('inf'),
-        t=None, halo_model = None , **kwargs):
+                  t=None, halo_model=None, **kwargs):
     """Differential rate per (ton year keV) of WIMP-nucleus scattering.
     :param es: Recoil energies in keV
     :param mw: WIMP mass in GeV/c^2
     :param sigma_nucleon: WIMP-nucleon cross-section in cm^2
     :param m_med: Medator mass in GeV/c^2. If not given, assumed very heavy.
     :param t: A J2000.0 timestamp. If not given,
     conservative velocity distribution is used.
-    :function halo_model : class (similar to the standard halo model) giving velocity distribution and dark matter density
+    :function halo_model : class (similar to the standard
+    halo model) giving velocity distribution and dark matter density
     :returns: numpy array of same length as es
 
     Further arguments are as for rate_wimp; see docstring of rate_wimp.
     """
-    halo_model = wr.standard_halo_model() if halo_model is None else halo_model
+    halo_model = wr.StandardHaloModel() if halo_model is None else halo_model
     return (rate_wimp(es=es * nu.keV,
                       mw=mw * nu.GeV/nu.c0**2,
                       sigma_nucleon=sigma_nucleon * nu.cm**2,
-                      m_med=m_med * nu.GeV/nu.c0**2, 
+                      m_med=m_med * nu.GeV/nu.c0**2,
                       t=t, halo_model=halo_model, **kwargs)
             * (nu.keV * (1000 * nu.kg) * nu.year))