Merge pull request #1475 from dlakaplan/add_site_clock_to_simulation

Update of #1165 (add site clock to simulation) to fix BIPM

CHANGELOG.md

@@ -9,6 +9,7 @@ and this project, at least loosely, adheres to [Semantic Versioning](https://sem
 ### Added
 - method to identify mask parameters with no TOAs and optionally freeze them
 ### Fixed
+- Creating fake TOAs properly handles site clock corrections
 - corrected a precision issue with reading ASCII representations of pulse profiles
 - Fixed matplotlib 3.6 import issue in pintk
 ### Removed

src/pint/simulation.py

@@ -64,10 +64,7 @@ def zero_residuals(ts, model, maxiter=10, tolerance=None):
     ts.compute_pulse_numbers(model)
     maxresid = None
     if tolerance is None:
-        if pint.utils.check_longdouble_precision():
-            tolerance = 1 * u.ns
-        else:
-            tolerance = 5 * u.us
+        tolerance = 1 * u.ns if pint.utils.check_longdouble_precision() else 5 * u.us
     for i in range(maxiter):
         r = pint.residuals.Residuals(ts, model, track_mode="use_pulse_numbers")
         resids = r.calc_time_resids(calctype="taylor")
@@ -85,13 +82,11 @@ def zero_residuals(ts, model, maxiter=10, tolerance=None):
         )
 
 
-def update_fake_toa_clock(ts, model, include_bipm=False, include_gps=True):
-    """Update the clock settings (corrections, etc) for fake TOAs
+def get_fake_toa_clock_versions(model, include_bipm=False, include_gps=True):
+    """Get the clock settings (corrections, etc) for fake TOAs
 
     Parameters
     ----------
-    ts : pint.toa.TOAs
-        Input TOAs (modified in-place)
     model : pint.models.timing_model.TimingModel
         current model
     include_bipm : bool, optional
@@ -111,7 +106,6 @@ def update_fake_toa_clock(ts, model, include_bipm=False, include_gps=True):
             if len(clk) == 2:
                 ctype, cvers = clk
                 if ctype == "TT" and cvers.startswith("BIPM"):
-                    include_bipm = True
                     if bipm_version is None:
                         bipm_version = cvers
                         log.info(f"Using CLOCK = {bipm_version} from the given model")
@@ -120,19 +114,19 @@ def update_fake_toa_clock(ts, model, include_bipm=False, include_gps=True):
                         f'CLOCK = {model["CLOCK"].value} is not implemented. '
                         f"Using TT({bipm_default}) instead."
                     )
+                include_bipm = True
         else:
             log.warning(
                 f'CLOCK = {model["CLOCK"].value} is not implemented. '
                 f"Using TT({bipm_default}) instead."
             )
+            include_bipm = True
 
-    ts.clock_corr_info.update(
-        {
-            "include_bipm": include_bipm,
-            "bipm_version": bipm_version,
-            "include_gps": include_gps,
-        }
-    )
+    return {
+        "include_bipm": include_bipm,
+        "bipm_version": bipm_version,
+        "include_gps": include_gps,
+    }
 
 
 def make_fake_toas(ts, model, add_noise=False, name="fake"):
@@ -257,17 +251,23 @@ def make_fake_toas_uniform(
         pint.toa.TOA(t.value, obs=obs, freq=f, scale=get_observatory(obs).timescale)
         for t, f in zip(times, freq_array)
     ]
-    ts = pint.toa.TOAs(toalist=t1)
-    ts.planets = model["PLANET_SHAPIRO"].value
-    ts.ephem = model["EPHEM"].value
-    update_fake_toa_clock(ts, model, include_bipm=include_bipm, include_gps=include_gps)
+    clk_version = get_fake_toa_clock_versions(
+        model, include_bipm=include_bipm, include_gps=include_gps
+    )
+    ts = pint.toa.get_TOAs_list(
+        toa_list=t1,
+        ephem=model["EPHEM"].value,
+        include_bipm=clk_version["include_bipm"],
+        bipm_version=clk_version["bipm_version"],
+        include_gps=clk_version["include_gps"],
+        planets=model["PLANET_SHAPIRO"].value,
+    )
+
     ts.table["error"] = error
     if dm is not None:
         for f in ts.table["flags"]:
             f["pp_dm"] = str(dm.to_value(pint.dmu))
             f["pp_dme"] = str(dm_error.to_value(pint.dmu))
-    ts.compute_TDBs()
-    ts.compute_posvels()
     return make_fake_toas(ts, model=model, add_noise=add_noise, name=name)
 
 
@@ -340,17 +340,22 @@ def make_fake_toas_fromMJDs(
         pint.toa.TOA(t.value, obs=obs, freq=f, scale=get_observatory(obs).timescale)
         for t, f in zip(times, freq_array)
     ]
-    ts = pint.toa.TOAs(toalist=t1)
-    ts.planets = model["PLANET_SHAPIRO"].value
-    ts.ephem = model["EPHEM"].value
-    update_fake_toa_clock(ts, model, include_bipm=include_bipm, include_gps=include_gps)
+    clk_version = get_fake_toa_clock_versions(
+        model, include_bipm=include_bipm, include_gps=include_gps
+    )
+    ts = pint.toa.get_TOAs_list(
+        toa_list=t1,
+        ephem=model["EPHEM"].value,
+        include_bipm=clk_version["include_bipm"],
+        bipm_version=clk_version["bipm_version"],
+        include_gps=clk_version["include_gps"],
+        planets=model["PLANET_SHAPIRO"].value,
+    )
     ts.table["error"] = error
     if dm is not None:
         for f in ts.table["flags"]:
             f["pp_dm"] = str(dm.to_value(pint.dmu))
             f["pp_dme"] = str(dm_error.to_value(pint.dmu))
-    ts.compute_TDBs()
-    ts.compute_posvels()
     return make_fake_toas(ts, model=model, add_noise=add_noise, name=name)
 
 

tests/test_fake_toas.py

@@ -10,6 +10,68 @@
 import pint.config
 from pint.fitter import GLSFitter
 from pinttestdata import datadir
+import pytest
+
+
+def roundtrip(toas, model):
+    with tempfile.NamedTemporaryFile("wt") as f:
+        toas.write_TOA_file(f)
+        f.flush()
+        toas2 = get_TOAs(f.name, model=model)
+    return toas2
+
+
+NGC6440E = """
+PSR              1748-2021E
+RAJ       17:48:52.75  1
+DECJ      -20:21:29.0  1
+F0       61.485476554  1
+F1         -1.181D-15  1
+PEPOCH        53750.000000
+POSEPOCH      53750.000000
+DM              223.9  1
+SOLARN0               0.00
+EPHEM               DE421
+UNITS               TDB
+TIMEEPH             FB90
+T2CMETHOD           TEMPO
+CORRECT_TROPOSPHERE N
+DILATEFREQ          N
+TZRMJD  53801.38605118223
+TZRFRQ            1949.609
+TZRSITE                  1
+"""
+
+
+@pytest.mark.parametrize(
+    "clock, planet",
+    [
+        ("UTC(NIST)", "Y"),
+        ("UTC(NIST)", "N"),
+        ("TT(BIPM2021)", "Y"),
+        ("TT(BIPM2021)", "N"),
+        ("TT(TAI)", "Y"),
+        ("TT(TAI)", "N"),
+    ],
+)
+def test_roundtrip(clock, planet):
+    # test for roundtrip accuracy at high precision
+    partext = f"{NGC6440E}\nCLK {clock}\nPLANET_SHAPIRO {planet}"
+    model = get_model(io.StringIO(partext))
+
+    toas = pint.simulation.make_fake_toas_uniform(
+        startMJD=56000,
+        endMJD=56400,
+        ntoas=100,
+        model=model,
+        obs="gbt",
+        error=1 * u.microsecond,
+        freq=1400 * u.MHz,
+    )
+    res = pint.residuals.Residuals(toas, model)
+    toas2 = roundtrip(toas, model)
+    res2 = pint.residuals.Residuals(toas2, model)
+    assert np.allclose(res.time_resids, res2.time_resids)
 
 
 def test_noise_addition():

tests/test_fitter_error_checking.py

@@ -232,4 +232,6 @@ def test_update_model_sets_things(Fitter):
     assert re.search(r"EPHEM *DE421", par_out)
     assert re.search(r"DMDATA *N", par_out)
     assert re.search(r"START *58000.0", par_out)
-    assert re.search(r"FINISH *59000.0", par_out)
+    assert re.search(r"FINISH *59000.0", par_out) or re.search(
+        r"FINISH *58999.9+", par_out
+    )