Merge pull request #1614 from dlakaplan/simplanets

Simulating TOAs from tim file when PLANET_SHAPIRO is true now works

CHANGELOG-unreleased.md

@@ -18,6 +18,7 @@ the released changes.
 - `pint.print_info()` function for bug reporting
 ### Fixed
 - Deleting JUMP1 from flag tables will not prevent fitting
+- Simulating TOAs from tim file when PLANET_SHAPIRO is true now works
 - Docstrings for `get_toas()` and `get_model_and_toas()`
 - Set `DelayComponent_list` and `NoiseComponent_list` to empty list if such components are absent
 - Fix invalid access of `PLANET_SHAPIRO` in models without `Astrometry`

src/pint/simulation.py

@@ -426,7 +426,7 @@ def make_fake_toas_fromtim(timfile, model, add_noise=False, name="fake"):
     --------
     :func:`make_fake_toas`
     """
-    input_ts = pint.toa.get_TOAs(timfile)
+    input_ts = pint.toa.get_TOAs(timfile, planets=model.PLANET_SHAPIRO.value)
 
     if input_ts.is_wideband():
         dm_errors = input_ts.get_dm_errors()

tests/test_fake_toas.py

@@ -259,7 +259,7 @@ def test_fake_uniform(t1, t2):
     assert np.isclose(r.calc_time_resids().std(), 1 * u.us, rtol=0.2)
 
 
-def test_fake_from_timfile():
+def test_fake_from_toas():
     # FIXME: this file is unnecessarily huge
     m, t = get_model_and_toas(
         pint.config.examplefile("B1855+09_NANOGrav_9yv1.gls.par"),
@@ -279,6 +279,29 @@ def test_fake_from_timfile():
     )
 
 
+@pytest.mark.parametrize("planets", (True, False))
+def test_fake_from_timfile(planets):
+    m = get_model(pint.config.examplefile("NGC6440E.par.good"))
+    t = get_TOAs(pint.config.examplefile("NGC6440E.tim"), planets=planets)
+
+    m.PLANET_SHAPIRO.value = planets
+
+    r = pint.residuals.Residuals(t, m)
+
+    t_sim = pint.simulation.make_fake_toas_fromtim(
+        pint.config.examplefile("NGC6440E.tim"), m, add_noise=True
+    )
+    r_sim = pint.residuals.Residuals(t_sim, m)
+
+    m, t = get_model_and_toas(
+        pint.config.examplefile("B1855+09_NANOGrav_9yv1.gls.par"),
+        pint.config.examplefile("B1855+09_NANOGrav_9yv1.tim"),
+    )
+    assert np.isclose(
+        r.calc_time_resids().std(), r_sim.calc_time_resids().std(), rtol=2
+    )
+
+
 def test_fake_highF1():
     m = get_model(os.path.join(datadir, "ngc300nicer.par"))
     m.F1.quantity *= 10