diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index f5ce22c22..819c1da50 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -9,17 +9,19 @@ the released changes. ## Unreleased ### Changed -- WAVE parameters can be added to a Wave model with `add_wave_component()` in wave.py +- `WAVE` parameters can be added to a `Wave` model with `add_wave_component()` in `wave.py` - Moved design matrix normalization code from `pint.fitter` to the new `pint.utils.normalize_designmatrix()` function. - Made `Residuals` independent of `GLSFitter` (GLS chi2 is now computed using the new function `Residuals._calc_gls_chi2()`). ### Added - Added WaveX model as DelayComponent with wave amplitudes as fitted parameters +- Added radial velocity methods for binary models - Support for wideband data in `pint.bayesian` (no correlated noise). ### Fixed - Wave model `validate()` can correctly use PEPOCH to assign WAVEEPOCH parameter - Fixed RTD by specifying theme explicitly. - `.value()` now works for pairParameters - Setting `model.PARAM1 = model.PARAM2` no longer overrides the name of `PARAM1` -- Fixed an incorrect docstring in pbprime() functions. +- Fixed an incorrect docstring in `pbprime()` functions. - Fix ICRS -> ECL conversion when parameter uncertainties are not set. +- `get_TOAs` raises an exception upon finding mixed narrowband and wideband TOAs in a tim file. `TOAs.is_wideband` returns True only if *ALL* TOAs have the -pp_dm flag. ### Removed diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index 1f02b8c75..6cc845412 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -799,6 +799,88 @@ def orbital_phase(self, barytimes, anom="mean", radians=True): # return with radian units or return as unitless cycles from 0-1 return anoms * u.rad if radians else anoms / (2 * np.pi) + def pulsar_radial_velocity(self, barytimes): + """Return line-of-sight velocity of the pulsar relative to the system barycenter at barycentric MJD times. + + Parameters + ---------- + barytimes: Time, TOAs, array-like, or float + MJD barycentric time(s). The times to compute the + orbital phases. Needs to be a barycentric time in TDB. + If a TOAs instance is passed, the barycentering will happen + automatically. If an astropy Time object is passed, it must + be in scale='tdb'. If an array-like object is passed or + a simple float, the time must be in MJD format. + + Raises + ------ + ValueError + If an astropy Time object is passed with scale!="tdb". + + Returns + ------- + array + The line-of-sight velocity + + Notes + ----- + This is the radial velocity of the pulsar. + + See [1]_ + + .. [1] Lorimer & Kramer, 2008, "The Handbook of Pulsar Astronomy", Eqn. 8.24 + """ + # this should also update the binary instance + nu = self.orbital_phase(barytimes, anom="true") + b = self.components[ + [x for x in self.components.keys() if x.startswith("Binary")][0] + ] + bbi = b.binary_instance # shorthand + psi = nu + bbi.omega() + return ( + 2 + * np.pi + * bbi.a1() + / (bbi.pb() * np.sqrt(1 - bbi.ecc() ** 2)) + * (np.cos(psi) + bbi.ecc() * np.cos(bbi.omega())) + ).cgs + + def companion_radial_velocity(self, barytimes, massratio): + """Return line-of-sight velocity of the companion relative to the system barycenter at barycentric MJD times. + + Parameters + ---------- + barytimes: Time, TOAs, array-like, or float + MJD barycentric time(s). The times to compute the + orbital phases. Needs to be a barycentric time in TDB. + If a TOAs instance is passed, the barycentering will happen + automatically. If an astropy Time object is passed, it must + be in scale='tdb'. If an array-like object is passed or + a simple float, the time must be in MJD format. + massratio : float + Ratio of pulsar mass to companion mass + + + Raises + ------ + ValueError + If an astropy Time object is passed with scale!="tdb". + + Returns + ------- + array + The line-of-sight velocity + + Notes + ----- + This is the radial velocity of the companion. + + See [1]_ + + .. [1] Lorimer & Kramer, 2008, "The Handbook of Pulsar Astronomy", Eqn. 8.24 + """ + return -self.pulsar_radial_velocity(barytimes) * massratio + def conjunction(self, baryMJD): """Return the time(s) of the first superior conjunction(s) after baryMJD. diff --git a/src/pint/toa.py b/src/pint/toa.py index 014403918..682ca60ce 100644 --- a/src/pint/toa.py +++ b/src/pint/toa.py @@ -321,6 +321,14 @@ def get_TOAs( if "pulse_number" in t.table.colnames and not include_pn: log.warning("'pulse_number' column exists but not being read in") t.remove_pulse_numbers() + + dm_data, valid_data = t.get_flag_value("pp_dm", as_type=float) + if len(valid_data) not in [0, len(t)]: + raise ValueError( + "Mixing narrowband and wideband toas in a TOAs object is not allowed. " + "Make sure that either all or no TOAs have the -pp_dm flag." + ) + return t @@ -1631,7 +1639,7 @@ def is_wideband(self): # there may be a more elegant way to do this dm_data, valid_data = self.get_flag_value("pp_dm", as_type=float) - return valid_data != [] + return len(valid_data) == len(self) def get_all_flags(self): """Return a list of all the flags used by any TOA.""" @@ -2986,4 +2994,5 @@ def get_TOAs_array( t.compute_TDBs(method=tdb_method, ephem=ephem) if "ssb_obs_pos" not in t.table.colnames: t.compute_posvels(ephem, planets) + return t diff --git a/tests/test_rv.py b/tests/test_rv.py new file mode 100644 index 000000000..6ea9d4293 --- /dev/null +++ b/tests/test_rv.py @@ -0,0 +1,43 @@ +import pytest +import numpy as np +from astropy import units as u, constants as c + +from pint.models import get_model +import pint.simulation +import pint.binaryconvert +from pinttestdata import datadir +import os + + +class TestRV: + def setup_method(self): + self.m = get_model( + os.path.join(datadir, "B1855+09_NANOGrav_dfg+12_modified_DD.par") + ) + # set the eccentricity to nonzero + # but not too high so that the models with various approximations won't fail + self.m.ECC.value = 0.01 + self.ts = self.m.T0.value + np.linspace(0, 2 * self.m.PB.value, 101) + self.ts = pint.simulation.make_fake_toas_fromMJDs(self.ts, self.m) + + self.v = self.m.pulsar_radial_velocity(self.ts) + + def test_rv_basemodel(self): + nu = self.m.orbital_phase(self.ts, anom="true", radians=True) + # HBOPA 8.24 + v = ( + (2 * np.pi / self.m.PB.quantity) + * (self.m.A1.quantity / np.sqrt(1 - self.m.ECC.quantity**2)) + * ( + np.cos(nu + self.m.OM.quantity) + + self.m.ECC.quantity * np.cos(self.m.OM.quantity) + ) + ) + assert np.allclose(self.v, v) + + @pytest.mark.parametrize("othermodel", ["ELL1", "ELL1H", "DDS", "BT"]) + def test_rv_othermodels(self, othermodel): + mc = pint.binaryconvert.convert_binary(self.m, othermodel) + vc = mc.pulsar_radial_velocity(self.ts) + # have a generous tolerance here since some of the models don't work well for high ECC + assert np.allclose(vc, self.v, atol=20 * u.km / u.s, rtol=1e-2) diff --git a/tests/test_toa.py b/tests/test_toa.py index 20303f174..7b367ba6d 100644 --- a/tests/test_toa.py +++ b/tests/test_toa.py @@ -182,3 +182,15 @@ def test_merge_toas(): toas_out = pint.toa.merge_TOAs([toas, toas2]) toas_outb = toas + toas2 assert np.all(toas_out.table == toas_outb.table) + + +def test_mix_nb_wb(): + with pytest.raises(ValueError): + t1 = pint.toa.get_TOAs( + io.StringIO( + """ + fake.ff 1430.000000 53393.561383615118386 0.178 ao -fe L-wide -be ASP -pp_dm 10.0 + fake.ff 1430.000000 53394.561383615118386 0.178 ao -fe L-wide -be ASP + """ + ) + ) diff --git a/tests/test_wideband_dm_data.py b/tests/test_wideband_dm_data.py index c62d5d9f3..32df2af5d 100644 --- a/tests/test_wideband_dm_data.py +++ b/tests/test_wideband_dm_data.py @@ -60,19 +60,6 @@ def wb_model(tmpdir): return get_model(str(parfile)) -@pytest.fixture -def wb_toas(wb_model): - toas = get_TOAs(io.StringIO(tim)) - for i in range(9): - r = Residuals(toas, wb_model) - if np.all(r.time_resids < 1 * u.ns): - break - toas.adjust_TOAs(TimeDelta(-r.time_resids)) - else: - raise ValueError - return toas - - @pytest.fixture def wb_toas_all(wb_model): toas = get_TOAs(io.StringIO(tim_all)) @@ -162,14 +149,18 @@ def test_dmjump_derivative(self): assert self.model.dm_derivs[dmj_param.name] == [self.model.d_dm_d_dmjump] -def test_wideband_residuals(wb_model, wb_toas): - r = WidebandTOAResiduals(wb_toas, wb_model, dm_resid_args=dict(subtract_mean=False)) - assert len(r.toa.time_resids) == len(wb_toas) - assert len(r.dm.dm_data) < len(wb_toas) +def test_wideband_residuals(wb_model, wb_toas_all): + r = WidebandTOAResiduals( + wb_toas_all, wb_model, dm_resid_args=dict(subtract_mean=False) + ) + assert len(r.toa.time_resids) == len(wb_toas_all) + assert len(r.dm.dm_data) == len(wb_toas_all) -def test_wideband_residuals_dmjump(wb_model, wb_toas): - r = WidebandTOAResiduals(wb_toas, wb_model, dm_resid_args=dict(subtract_mean=False)) +def test_wideband_residuals_dmjump(wb_model, wb_toas_all): + r = WidebandTOAResiduals( + wb_toas_all, wb_model, dm_resid_args=dict(subtract_mean=False) + ) model = deepcopy(wb_model) assert wb_model.DMJUMP1.value == 0 model.DMJUMP1.value = 10 @@ -178,10 +169,17 @@ def test_wideband_residuals_dmjump(wb_model, wb_toas): model.DMJUMP0 with pytest.raises(AttributeError): model.DMJUMP2 - r2 = WidebandTOAResiduals(wb_toas, model, dm_resid_args=dict(subtract_mean=False)) + r2 = WidebandTOAResiduals( + wb_toas_all, model, dm_resid_args=dict(subtract_mean=False) + ) assert 0 < np.sum(r.dm.resids_value != r2.dm.resids_value) < len(r.dm.resids_value) +def test_read_mixed_timfile(): + with pytest.raises(ValueError): + get_TOAs(io.StringIO(tim)) + + def test_wideband_residuals_dof(wb_model, wb_toas_all): wb_model.free_params = ["DMJUMP1"] r = WidebandTOAResiduals( @@ -192,14 +190,6 @@ def test_wideband_residuals_dof(wb_model, wb_toas_all): assert_allclose(r.reduced_chi2, r.chi2 / r.dof) -@pytest.mark.xfail(reason="All TOAs must have DMs, currently") -def test_wideband_fit_dmjump(wb_model, wb_toas): - wb_model.free_params = ["DMJUMP1"] - fitter = WidebandTOAFitter(wb_toas, wb_model) - fitter.fit_toas() - assert_allclose(fitter.model.DMJUMP1.value, -10, atol=1e-3) - - def test_wideband_fit_dmjump_all(wb_model, wb_toas_all): wb_model.free_params = ["DMJUMP1"] fitter = WidebandTOAFitter(wb_toas_all, wb_model)