diff --git a/CHANGELOG.md b/CHANGELOG.md index 74b1c36f8..6f667ba1e 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/src/pint/simulation.py b/src/pint/simulation.py index 65bba9be4..f14d38b1e 100644 --- a/src/pint/simulation.py +++ b/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) diff --git a/tests/test_fake_toas.py b/tests/test_fake_toas.py index 982db2526..503a66ec5 100644 --- a/tests/test_fake_toas.py +++ b/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(): diff --git a/tests/test_fitter_error_checking.py b/tests/test_fitter_error_checking.py index f1ec1e4eb..331bcd3d0 100644 --- a/tests/test_fitter_error_checking.py +++ b/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 + )