Remove hardcoded v_0 from earth_velocity (#14)

* remove hardcoded v0 * Update halo.py * fix references * make the example * update defaults * ah none defaults * remove messy compare * remove unused arg * pin numpy * fix syntax * Update pytest.yml * Update test_wimprates.py * Update test_wimprates.py --------- Co-authored-by: Joran Angevaare <[email protected]>
JelleAalbers · Feb 13, 2023 · e3c42bb · e3c42bb
1 parent e3aa729
commit e3c42bb
Show file tree

Hide file tree

Showing 6 changed files with 150 additions and 83 deletions.
diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml
@@ -36,6 +36,7 @@ jobs:
         run: |
           # Requirements for running the notebooks as a pytest
           pip install cython ipython
+          pip install "numpy<1.23"
           pip install nbconvert==6.4.3 nbmake pytest-xdist pytest coverage coveralls pytest-cov pytest-notebook ipython_genutils
           # Several optional packages that are imported in the notebooks
           pip install git+https://github.com/XENON1T/laidbax

diff --git a/setup.py b/setup.py
@@ -19,7 +19,7 @@
     package_dir={'wimprates': 'wimprates'},
     package_data={'wimprates': [
         'data/bs/*', 'data/migdal/*', 'data/sd/*', 'data/dme/*']},
-    tests_require=requirements + ['pytest'],
+    tests_require=requirements + ['pytest', 'unittest'],
     keywords='wimp,spin-dependent,spin-independent,bremsstrahlung,migdal',
     classifiers=['Intended Audience :: Science/Research',
                  'Development Status :: 3 - Alpha',

diff --git a/tests/test_halo.py b/tests/test_halo.py
@@ -8,7 +8,7 @@
 
 def test_shm_values():
     halo_model = StandardHaloModel()
-    assert np.abs(halo_model.v_0 /(nu.km/nu.s) - 220.)<1e-6
+    assert np.abs(halo_model.v_0 /(nu.km/nu.s) - 238.)<1e-6
     assert np.abs(halo_model.v_esc /(nu.km/nu.s) - 544.)<1e-6
 
 

diff --git a/tests/test_wimprates.py b/tests/test_wimprates.py
@@ -8,60 +8,99 @@
 import numpy as np
 
 import wimprates as wr
-
-
-opts = dict(mw=50, sigma_nucleon=1e-45)
-
-
-def isclose(x, y):
-    assert np.abs(x - y)/x < 1e-5
-
-
-def test_elastic():
-    ref = 33.19098343826968
-
-    isclose(wr.rate_wimp_std(1, **opts), ref)
-
-    # Test numericalunits.reset_units() does not affect results
-    nu.reset_units(123)
-    isclose(wr.rate_wimp_std(1, **opts), ref)
-
-    # Test vectorized call
-    energies = np.linspace(0.01, 40, 100)
-    dr = wr.rate_wimp_std(energies, **opts)
-    assert dr[0] == wr.rate_wimp_std(0.01, **opts)
-
-
-def test_lightmediator():
-    isclose(wr.rate_wimp_std(1, m_med=1e-3, **opts),
-            0.0005502663384403058)
-
-
-def test_spindependent():
-    isclose(wr.rate_wimp_std(1, interaction='SD_n_central', **opts),
-            0.00021779266679860948)
-
-
-def test_migdal():
-    isclose(wr.rate_wimp_std(1, detection_mechanism='migdal', **opts),
-            0.2610240963512907)
-
-
-def test_brems():
-    isclose(wr.rate_wimp_std(1, detection_mechanism='bremsstrahlung', **opts),
-            0.00017062652972332665)
-
-
-def test_dme():
-    isclose(
-        wr.rate_dme(100* nu.eV, 4, 'd',
-                    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.StandardHaloModel(rho_dm = 0.3 * nu.GeV/nu.c0**2 / nu.cm**3) , **opts), ref)
-
+import unittest
+
+
+class TestBenchmarks(unittest.TestCase):
+    opts = dict(mw=50,
+                sigma_nucleon=1e-45,
+                )
+    def test_elastic(self):
+        ref = 30.39515403337126
+
+        self.assertAlmostEqual(wr.rate_wimp_std(1, **self.opts), ref)
+
+        # Test numericalunits.reset_units() does not affect results
+        nu.reset_units(123)
+        self.assertAlmostEqual(wr.rate_wimp_std(1, **self.opts), ref)
+
+        # Test vectorized call
+        energies = np.linspace(0.01, 40, 100)
+        dr = wr.rate_wimp_std(energies, **self.opts)
+        self.assertEqual(dr[0], wr.rate_wimp_std(0.01, **self.opts))
+
+
+    def test_lightmediator(self):
+        self.assertAlmostEqual(wr.rate_wimp_std(1, m_med=1e-3, **self.opts),
+                0.0005039148255734496)
+
+
+    def test_spindependent(self):
+        self.assertAlmostEqual(wr.rate_wimp_std(1, interaction='SD_n_central', **self.opts),
+                0.00019944698779638946)
+
+
+    def test_migdal(self):
+        self.assertAlmostEqual(wr.rate_wimp_std(1, detection_mechanism='migdal', **self.opts),
+                0.27459766238555017)
+
+
+    def test_brems(self):
+        self.assertAlmostEqual(wr.rate_wimp_std(1, detection_mechanism='bremsstrahlung', **self.opts),
+                0.00017949417705393473)
+
+
+    def test_dme(self):
+        self.assertAlmostEqual(
+            wr.rate_dme(100* nu.eV, 4, 'd',
+                        mw=nu.GeV/nu.c0**2, sigma_dme=4e-44 * nu.cm**2)
+                * nu.kg * nu.keV * nu.day,
+        2.87553027086139e-06)
+
+    def test_halo_scaling(self):
+        #check that passing rho multiplies the rate correctly:
+        ref = 30.39515403337113
+        halo_model = wr.StandardHaloModel(rho_dm=0.3 * nu.GeV / nu.c0 ** 2 / nu.cm ** 3)
+        self.assertAlmostEqual(wr.rate_wimp_std(1,
+                                                halo_model=halo_model,
+                                                **self.opts,
+                                                ), ref)
+
+
+    def test_v_earth_old(self):
+        """
+        Check that v_0 = 220 gives 232 km/s for v_earth (old convention)
+        https://arxiv.org/pdf/hep-ph/0504010.pdf
+        """
+        kms = nu.km/nu.s
+        v_0_old = 220 * kms
+
+        # Unfortunately, we actually don't get the right answer in this test.
+        # the reason is that in the old convention, the years average 
+        # was approximated by a different date from the current convention.
+        # See https://github.com/JelleAalbers/wimprates/pull/14
+        # We're now going to to just add the factor that we are of by, for
+        # bookkeeping
+        wrong_by_this_much = 2.40803047754608
+
+        self.assertAlmostEqual(232 + wrong_by_this_much,
+                               wr.v_earth(t=None, v_0=v_0_old)/kms,
+                               )
+
+    def test_average_v_earth(self):
+        """
+        In v_earth a day is set that should return ~ the annual average, test
+        this is true
+        """
+        kms = nu.km / nu.s
+        v_0_default = 238 * kms
+        # Lazy way of averaging over one year
+        v_earth_average = np.mean(
+            [wr.v_earth(t, v_0=v_0_default)
+             for t in np.linspace(0, 365.25, int(1e4))]
+        ) / kms
+        self.assertAlmostEqual(
+            v_earth_average, wr.v_earth(t=None, v_0=v_0_default) / kms,
+            # places=1 means that we get the same results at the first decimal (fine for 500.0<?>)
+            places=1
+        )
diff --git a/wimprates/__init__.py b/wimprates/__init__.py
@@ -1,9 +1,18 @@
 __version__ = '0.4.1'
 
+from packaging import version
+if version.parse('0.4.1') < version.parse(__version__) < version.parse('0.6.0'):
+    # Remove this warning at some point
+    import warnings
+    warnings.warn(
+        'Default WIMP parameters are changed in accordance with '
+        'https://arxiv.org/abs/2105.00599 (github.com/JelleAalbers/wimprates/pull/14)')
+
 from .utils import *
 from .halo import *
 from .elastic_nr import *
 from .bremsstrahlung import *
 from .migdal import *
 from .electron import *
 from .summary import *
+
diff --git a/wimprates/halo.py b/wimprates/halo.py
@@ -7,10 +7,20 @@
 from scipy.special import erf
 
 
+
 import wimprates as wr
 export, __all__ = wr.exporter()
 
 
+# See https://arxiv.org/abs/2105.00599 and references therein
+_HALO_DEFAULTS = dict(
+    rho_dm = 0.3, # GeV / c2 / cm3
+    v_esc = 544,  # km/s
+    v_orbit = 29.8,  # km/s
+    v_pec = (11.1, 12.2, 7.3),  # km/s
+    v_0 = 238,  # km/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).
@@ -59,14 +69,18 @@ def j2000_from_ymd(year, month, day_of_month):
 
 
 @export
-def earth_velocity(t):
+def earth_velocity(t, v_0 = None):
     """Returns 3d velocity of earth, in the galactic rest frame,
     in galactic coordinates.
     :param t: J2000.0 timestamp
+    :param v_0: Local standard of rest velocity
 
     Values and formula from https://arxiv.org/abs/1209.3339
     Assumes earth circular orbit.
     """
+    if v_0 is None :
+        v_0 = _HALO_DEFAULTS['v_0'] * nu.km/nu.s
+
     # e_1 and e_2 are the directions of earth's velocity at t1
     # and t1 + 0.25 year.
     e_1 = np.array([0.9931, 0.1170, -0.01032])
@@ -79,41 +93,45 @@ def earth_velocity(t):
     phi = omega * (t - t1)
 
     # Mean orbital velocity of the Earth (Lewin & Smith appendix B)
-    v_orbit = 29.79 * nu.km / nu.s
+    v_orbit = _HALO_DEFAULTS['v_orbit'] * nu.km / nu.s
 
     v_earth_sun = v_orbit * (e_1 * np.cos(phi) + e_2 * np.sin(phi))
 
     # Velocity of Local Standard of Rest
-    v_lsr = np.array([0, 220, 0]) * nu.km/nu.s
+    v_lsr = np.array([0, v_0, 0])
     # Solar peculiar velocity
-    v_pec = np.array([11, 12, 7]) * nu.km/nu.s
+    v_pec = np.array(_HALO_DEFAULTS['v_pec']) * nu.km/nu.s
 
     return v_lsr + v_pec + v_earth_sun
 
 
 @export
-def v_earth(t=None):
+def v_earth(t=None, v_0=None):
     """Return speed of earth relative to galactic rest frame
+    Velocity of earth/sun relative to gal. center (eccentric orbit, so not
+    equal to v_0).
+
     :param t: J2000 timestamp or None
+    :param v_0: Local standard of rest velocity
     """
     if t is None:
-        # Velocity of earth/sun relative to gal. center
-        # (eccentric orbit, so not equal to v_0)
-        return 232 * nu.km / nu.s
-    else:
-        return np.sum(earth_velocity(t) ** 2) ** 0.5
+        # This day (Feb 29 2000) gives ~ the annual average speed
+        t = 59.37
+    return np.sum(earth_velocity(t, v_0=v_0) ** 2) ** 0.5
 
 
 @export
-def v_max(t=None, v_esc=None):
+def v_max(t=None, v_esc=None, v_0=None):
     """Return maximum observable dark matter velocity on Earth."""
-    v_esc = 544 * nu.km/nu.s if v_esc is None else v_esc  # default
+    # defaults
+    v_esc = _HALO_DEFAULTS['v_esc'] * nu.km/nu.s if v_esc is None else v_esc
+    v_0 = _HALO_DEFAULTS['v_0'] * nu.km / nu.s if v_0 is None else v_0
     # 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 + v_earth(t)
+        return v_esc + v_earth(t, v_0=v_0)
     else:
-        return v_esc + np.sum(earth_velocity(t) ** 2) ** 0.5
+        return v_esc + np.sum(earth_velocity(t, v_0=v_0) ** 2) ** 0.5
 
 
 @export
@@ -130,9 +148,9 @@ def observed_speed_dist(v, t=None, v_0=None, v_esc=None):
 
     Further inputs: scale velocity v_0 and escape velocity 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)
+    v_0 = _HALO_DEFAULTS['v_0'] * nu.km/nu.s if v_0 is None else v_0
+    v_esc = _HALO_DEFAULTS['v_esc'] * nu.km/nu.s if v_esc is None else v_esc
+    v_earth_t = v_earth(t, v_0=v_0)
 
     # Normalization constant, see Lewin&Smith appendix 1a
     _w = v_esc/v_0
@@ -155,13 +173,13 @@ def observed_speed_dist(v, t=None, v_0=None, v_esc=None):
         len(v)
     except TypeError:
         # Scalar argument
-        if v > v_max(t, v_esc):
+        if v > v_max(t, v_esc, v_0=v_0):
             return 0
         else:
             return y
     else:
         # Array argument
-        y[v > v_max(t, v_esc)] = 0
+        y[v > v_max(t, v_esc, v_0=v_0)] = 0
         return y
 
 
@@ -175,16 +193,16 @@ class used to pass a halo model to the rate computation
         giving normalised velocity 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: Local standard of rest 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
+        self.v_0 = _HALO_DEFAULTS['v_0'] * nu.km/nu.s if v_0 is None else v_0
+        self.v_esc = _HALO_DEFAULTS['v_esc'] * nu.km/nu.s if v_esc is None else v_esc
+        self.rho_dm = _HALO_DEFAULTS['rho_dm'] * 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)
+        return observed_speed_dist(v, t, v_0=self.v_0, v_esc=self.v_esc)