Test if both j2000 implementations are equal

tests/test_halo.py

@@ -1,5 +1,5 @@
 import pandas as pd
-from wimprates import j2000, StandardHaloModel, j200_from_ymd
+from wimprates import j2000, StandardHaloModel, j2000_from_ymd
 import numericalunits as nu
 import numpy as np
 
@@ -11,20 +11,39 @@ def test_shm_values():
 
 
 def test_j2000():
-    assert j200_from_ymd(2009, 1, 31.75) == 3318.25
+    assert j2000_from_ymd(2009, 1, 31.75) == 3318.25
 
 
 def test_j2000_datetime():
     date = pd.datetime(year=2009, month=1, day=31, hour=18)
-    assert j2000(date=date) == 3318.25
+    assert j2000(date) == 3318.25
 
 
 def test_j2000_ns_int():
     date = pd.datetime(year=2009, month=1, day=31, hour=18)
-    assert j2000(date=pd.to_datetime(date).value) == 3318.25
+    value = pd.to_datetime(date).value
+    assert isinstance(value, int)
+    assert j2000(value) == 3318.25
 
 
 def test_j2000_ns_array():
-    date = pd.datetime(year=2009, month=1, day=31, hour=18)
-    dates = np.array([pd.to_datetime(date).value] * 3)
-    np.testing.assert_array_equal(j2000(date=dates), np.array([3318.25] * 3))
+    # Generate some test cases and compare the two implementations
+    years = np.arange(2008, 2020)
+    months = np.arange(1, 12 + 1)
+    days = np.arange(1, 12 + 1)
+    hours = np.arange(1, 12 + 1)
+
+    dates = np.zeros(12)
+    j2000_ymd = np.zeros(12)
+    for i in range(len(j2000_ymd)):
+        print(years[i], months[i], days[i], hours[i])
+        date = pd.datetime(year=years[i],
+                           month=months[i],
+                           day=days[i],
+                           hour=hours[i])
+        dates[i] = pd.to_datetime(date).value
+        # Compute according to j2000_from_ymd
+        j2000_ymd[i] = j2000_from_ymd(years[i],
+                                      months[i],
+                                      days[i] + hours[i] / 24)
+    np.testing.assert_array_almost_equal(j2000(dates), j2000_ymd, decimal=6)

wimprates/halo.py

@@ -29,7 +29,7 @@ def j2000(date):
     """
     zero = pd.to_datetime('2000-01-01T12:00')
     nanoseconds_per_day = 1e9 * 3600 * 24
-    if isinstance(date, pd.datetime):
+    if isinstance(date, datetime):
         # pd.datetime refers to datetime.datetime
         # make it into a pd.Timestamp
         # Timestamp.value gives timestamp in ns
@@ -40,7 +40,7 @@ def j2000(date):
 
 
 @export
-def j200_from_ymd(year, month, day_of_month):
+def j2000_from_ymd(year, month, day_of_month):
     """"Returns the fractional number of days since J2000.0 epoch.
     :param year: Year
     :param month: Month (January = 1)
@@ -73,7 +73,7 @@ def earth_velocity(t):
     e_2 = np.array([-0.0670, 0.4927, -0.8676])
     # t1 is the time of the vernal equinox, March 21. Does it matter what
     # year? Precession of equinox takes 25800 years so small effect.
-    t1 = j2000(2000, 3, 21)
+    t1 = j2000_from_ymd(2000, 3, 21)
     # Angular frequency
     omega = 2 * np.pi / 365.25
     phi = omega * (t - t1)