Merge pull request #1861 from kerrm/fix_orbwave_units

Fix WAVE_OM units

CHANGELOG-unreleased.md

@@ -16,6 +16,7 @@ the released changes.
 - When TCB->TDB conversion info is missing, will print parameter name
 - `add_param` returns the name of the parameter (useful for numbered parameters)
 ### Fixed
+- Changed WAVE_OM units from 1/d to rad/d.
 - When EQUAD is created from TNEQ, has proper TCB->TDB conversion info
 - TOA selection masks will work when only TOA is the first one
 ### Removed

src/pint/models/wave.py

@@ -15,9 +15,8 @@ class Wave(PhaseComponent):
     sine/cosine components.
 
     For consistency with the implementation in tempo2, this signal is treated
-    as a time series, but trivially converted into phase by multiplication by
-    F0, which could makes changes to PEPOCH fragile if there is strong spin
-    frequency evolution.
+    as a time series and trivially converted into phase by multiplication by
+    F0.
 
     Parameters supported:
 
@@ -42,7 +41,7 @@ def __init__(self):
             floatParameter(
                 name="WAVE_OM",
                 description="Base frequency of wave solution",
-                units="1/d",
+                units="rad/d",
                 convert_tcb2tdb=False,
             )
         )

src/pint/utils.py

@@ -1745,11 +1745,8 @@ def _translate_wave_freqs(om: Union[float, u.Quantity], k: int) -> u.Quantity:
     astropy.units.Quantity
         WXFREQ_ quantity in units 1/d that can be used in WaveX model
     """
-    if isinstance(om, u.quantity.Quantity):
-        om.to(u.d**-1)
-    else:
-        om *= u.d**-1
-    return (om * (k + 1)) / (2.0 * np.pi)
+    om <<= u.rad / u.d
+    return (om * (k + 1)) / (2.0 * np.pi * u.rad)
 
 
 def _translate_wavex_freqs(wxfreq: Union[float, u.Quantity], k: int) -> u.Quantity:
@@ -1769,13 +1766,12 @@ def _translate_wavex_freqs(wxfreq: Union[float, u.Quantity], k: int) -> u.Quanti
     astropy.units.Quantity
         WAVEOM quantity in units 1/d that can be used in Wave model
     """
-    if isinstance(wxfreq, u.quantity.Quantity):
-        wxfreq.to(u.d**-1)
-    else:
-        wxfreq *= u.d**-1
+    wxfreq <<= u.d**-1
     if len(wxfreq) == 1:
-        return (2.0 * np.pi * wxfreq) / (k + 1.0)
-    wave_om = [((2.0 * np.pi * wxfreq[i]) / (k[i] + 1.0)) for i in range(len(wxfreq))]
+        return (2.0 * np.pi * u.rad * wxfreq) / (k + 1.0)
+    wave_om = [
+        ((2.0 * np.pi * u.rad * wxfreq[i]) / (k[i] + 1.0)) for i in range(len(wxfreq))
+    ]
     return (
         sum(wave_om) / len(wave_om)
         if np.allclose(wave_om, wave_om[0], atol=1e-3)

tests/datafile/vela_wave.par

@@ -28,45 +28,18 @@ TZRMJD         55896.55312516020475
 TZRFRQ         1370.2919919999999365    
 TZRSITE        pks   
 TRES           310.453      
-EPHVER         2                           
 CLK            TT(TAI)
 MODE 1
 UNITS          TDB
 TIMEEPH        FB90
 DILATEFREQ     N
 PLANET_SHAPIRO N
-T2CMETHOD      TEMPO
 NE_SW          9.961
 CORRECT_TROPOSPHERE  N
 EPHEM          DE405
 NITS           1
 NTOA           339
 CHI2R          328088.1294 298
-JUMP -B 20CM 0 0
-JUMP -B 40CM 0 0
-JUMP -B 50CM 0 0
-JUMP -dfb3_J0437_55319_56160 1 2.2e-07 0
-JUMP -dfb3_J0437_56160_60000 1 4.5e-07 0
-JUMP -pdfb1_128_ch 1 -3.5547e-06 0
-JUMP -pdfb1_2048_ch 1 -4.68829e-05 0
-JUMP -pdfb1_512_ch 1 -1.22813e-05 0
-JUMP -pdfb1_post_2006 1 -1.3e-07 0
-JUMP -pdfb1_pre_2006 1 -1.13e-06 0
-JUMP -pdfb2_1024_MHz 1 -5.435e-06 0
-JUMP -pdfb2_256MHz_1024_ch 1 -1.1395e-05 0
-JUMP -pdfb2_256MHz_512ch 1 -4.75e-06 0
-JUMP -pdfb3_1024_256_512 1 2.45e-06 0
-JUMP -pdfb3_1024_MHz 1 1.03e-06 0
-JUMP -pdfb3_256MHz_1024ch 1 4.295e-06 0
-JUMP -pdfb3_64MHz_1024ch 1 1.494e-05 0
-JUMP -pdfb3_64MHz_512ch 1 8.9e-06 0
-JUMP -pdfb4.*_1024_[1,2]... 1 2.23e-06 0
-JUMP -pdfb4_256MHz_1024ch 1 5.05e-06 0
-JUMP -pdfb4_55319_56055_cals 1 9.27e-07 0
-JUMP -pdfb4_56110_56160_cals 1 5.41e-07 0
-JUMP -pdfb4_56160_60000_cals 1 4.25e-07 0
-JUMP -wbb256_512_128_3p_b 1 -6.2e-07 0
-JUMP -wbb_c_config 1 3.8e-07 0
 WAVEEPOCH 55305
 WAVE_OM 0.0015182579855022 0
 WAVE1 -0.21573979911255 -0.049063841960712

tests/test_model_wave.py

@@ -1,28 +1,81 @@
+from io import StringIO
 import os
 import pytest
 
 import astropy.units as u
+import numpy as np
 
-import pint.fitter
-import pint.models
+from pint.models import get_model
+from pint import toa
 import pint.residuals
-import pint.toa
 from pinttestdata import datadir
 
-# Not included in the test here, but as a sanity check I used this same
-# ephemeris to phase up Fermi data, and it looks good.
+par_nowave = """
+    PSRJ           J0835-4510
+    RAJ            08:35:20.61149
+    DECJ           -45:10:34.8751
+    F0             11.18965156782
+    PEPOCH         55305                       
+    DM             67.99
+    UNITS          TDB
+"""
 
-parfile = os.path.join(datadir, "vela_wave.par")
-timfile = os.path.join(datadir, "vela_wave.tim")
+wave_terms = """
+    WAVEEPOCH 55305
+    WAVE_OM 0.0015182579855022
+    WAVE1 -0.21573979911255 -0.049063841960712
+    WAVE2 0.62795320246729 -0.11984954655989
+    WAVE3 0.099618608456845 0.28530756908788
+    WAVE4 -0.21537340649058 0.18849486610628
+    WAVE5 0.021980474493165 -0.23566696662127
+"""
 
 
-class TestWave:
-    @classmethod
-    def setup_class(cls):
-        cls.m = pint.models.get_model(parfile)
-        cls.t = pint.toa.get_TOAs(timfile, ephem="DE405", include_bipm=False)
+def test_wave_ephem():
+    parfile = os.path.join(datadir, "vela_wave.par")
+    timfile = os.path.join(datadir, "vela_wave.tim")
+    m = get_model(parfile)
+    t = toa.get_TOAs(timfile, ephem="DE405", include_bipm=False)
+    rs = pint.residuals.Residuals(t, m).time_resids
+    assert rs.std() < 350.0 * u.us
 
-    def test_vela_rms_is_small_enough(self):
-        rs = pint.residuals.Residuals(self.t, self.m).time_resids
-        rms = rs.to(u.us).std()
-        assert rms < 350.0 * u.us
+
+def test_wave_construction():
+    m = get_model(StringIO(par_nowave + wave_terms))
+    assert np.allclose(m.WAVE_OM.quantity, 0.0015182579855022 * u.rad / u.day)
+    assert np.allclose(m.WAVE1.quantity[0], -0.21573979911255 * u.s)
+    assert np.allclose(m.WAVE2.quantity[1], -0.1198495465598 * u.s)
+
+
+def test_wave_computation():
+    m0 = get_model(StringIO(par_nowave))
+    m1 = get_model(StringIO(par_nowave + wave_terms))
+    # make some barycentric TOAs
+    tdbs = np.linspace(54500, 60000, 10)
+    ts = toa.TOAs(
+        toalist=[
+            toa.TOA(t, obs="@", freq=np.inf, error=1 * u.ms, scale="tdb") for t in tdbs
+        ]
+    )
+    ts.compute_TDBs(ephem="DE421")
+    ts.compute_posvels(ephem="DE421")
+    ph0 = m0.phase(ts)
+    ph1 = m1.phase(ts)
+    dphi = (ph1.int - ph0.int) + (ph1.frac - ph0.frac)
+    test_phase = np.zeros(len(tdbs))
+    WAVEEPOCH = 55305
+    WAVE_OM = 0.0015182579855022
+    WAVE = [
+        [-0.21573979911255, -0.049063841960712],
+        [0.62795320246729, -0.11984954655989],
+        [0.099618608456845, 0.28530756908788],
+        [-0.21537340649058, 0.18849486610628],
+        [0.021980474493165, -0.23566696662127],
+    ]
+    ph = (tdbs - WAVEEPOCH) * WAVE_OM
+    for i in range(5):
+        test_phase += WAVE[i][0] * np.sin((i + 1) * ph) + WAVE[i][1] * np.cos(
+            (i + 1) * ph
+        )
+    test_phase *= m0.F0.quantity.to(u.Hz).value
+    assert np.allclose(test_phase, dphi)