trying to fix wimprates

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()
+    nu.reset_units("SI")
     isclose(wr.rate_wimp_std(1, **opts), ref)
 
     # Test vectorized call
@@ -58,3 +58,10 @@ def test_dme():
                     mw=nu.GeV/nu.c0**2, sigma_dme=4e-44 * nu.cm**2)
             * nu.kg * nu.keV * nu.day,
     2.232912243660405e-06)
+
+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)
+

wimprates/elastic_nr.py

@@ -206,8 +206,7 @@ def rate_elastic(erec, mw, sigma_nucleon, interaction='SI',
         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

@@ -86,11 +86,11 @@ def v_min_dme(eb, erec, q, mw, halo_model = wr.standard_halo_model()):
 
 
 # Precompute velocity integrals for t=None
-_v_mins = np.linspace(0, 1, 1000) * wr.v_max(None, halo_model.v_esc)
+_v_mins = np.linspace(0, 1, 1000) * wr.v_max(None, wr.v_esc())
 _ims = np.array([
-    quad(lambda v: 1 / v * halo_model.velocity_dist(v),
+    quad(lambda v: 1 / v * wr.observed_speed_dist(v),
          _v_min,
-         wr.v_max(None, halo_model.v_esc))[0]
+         wr.v_max(None, wr.v_esc() ))[0]
     for _v_min in _v_mins])
 
 # Store interpolator in km/s rather than unit-dependent numbers

wimprates/halo.py

@@ -136,17 +136,19 @@ def observed_speed_dist(v, t=None, v_0_value=v_0(), v_esc_value=v_esc()):
 
     # 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)
+    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)
                       / (2 * v_earth_t * v))
+    #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)
 
     # Zero if v > v_max
     try:
@@ -173,11 +175,11 @@ class used to pass a halo model to the rate computation
         The standard halo model also allows variation of v_0
         :param v_0 
     """
-    def __init__(self, v_0_value = v_0(), v_esc_value=v_esc(),rho_dm_value=rho_dm())
+    def __init__(self, v_0_value = v_0(), v_esc_value=v_esc(),rho_dm_value=rho_dm()):
         self.v_0 = v_0_value
         self.v_esc = v_esc_value
         self.rho_dm = rho_dm_value
     def velocity_dist(self,v,t):
         #in units of per velocity, 
         #v is in units of velocity 
-        return observed_speed_dist(v, t=None, self.v_0, self.v_esc)
+        return observed_speed_dist(v, t, self.v_0, self.v_esc)