From 6b286367a3da0e6cfeb3213d881431fd6dd35111 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Thu, 31 Aug 2023 10:47:19 -0500 Subject: [PATCH 01/14] is_wideband --- src/pint/toa.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/pint/toa.py b/src/pint/toa.py index 014403918..cc5742878 100644 --- a/src/pint/toa.py +++ b/src/pint/toa.py @@ -1631,7 +1631,14 @@ 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 != [] + + if len(valid_data) not in [0, len(self)]: + log.warning( + "Only a subset of the TOAs have wideband DM information." + "A TOAs object will be treated wideband only of ALL TOAs have DM information." + ) + + return len(valid_data) == len(self) def get_all_flags(self): """Return a list of all the flags used by any TOA.""" From fd72267d057bdb061d9c0b5b2a239d9c7cf95710 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Thu, 31 Aug 2023 10:54:44 -0500 Subject: [PATCH 02/14] test_is_wideband --- src/pint/toa.py | 4 ++-- tests/test_toa.py | 14 ++++++++++++++ 2 files changed, 16 insertions(+), 2 deletions(-) diff --git a/src/pint/toa.py b/src/pint/toa.py index cc5742878..210a167d2 100644 --- a/src/pint/toa.py +++ b/src/pint/toa.py @@ -1634,8 +1634,8 @@ def is_wideband(self): if len(valid_data) not in [0, len(self)]: log.warning( - "Only a subset of the TOAs have wideband DM information." - "A TOAs object will be treated wideband only of ALL TOAs have DM information." + "Only a subset of the TOAs have wideband DM information. " + "A TOAs object will be treated as wideband only of *ALL* TOAs have DM information." ) return len(valid_data) == len(self) diff --git a/tests/test_toa.py b/tests/test_toa.py index 20303f174..b262e198b 100644 --- a/tests/test_toa.py +++ b/tests/test_toa.py @@ -182,3 +182,17 @@ 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_is_wideband(): + 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 + """ + ) + ) + assert not t1.is_wideband() + assert t1[0].is_wideband() + assert not t1[1].is_wideband() From 74e63bf9f1a547ce405a4ec218e58c5d6b5382da Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Thu, 31 Aug 2023 11:04:48 -0500 Subject: [PATCH 03/14] changelog --- CHANGELOG-unreleased.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index ac9bc2495..e998a4d5b 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -11,6 +11,7 @@ the released changes. ### Changed - 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()`). +- `TOAs.is_wideband` returns True only if *ALL* TOAs have the -pp_dm flag. Emits a warning if only some TOAs have this flag. ### Added ### Fixed - Fixed RTD by specifying theme explicitly. From 0fff91943458920a67388202e3103e088a473043 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Thu, 31 Aug 2023 12:25:43 -0500 Subject: [PATCH 04/14] raise exception on read --- src/pint/toa.py | 16 +++++++++------- tests/test_toa.py | 20 +++++++++----------- 2 files changed, 18 insertions(+), 18 deletions(-) diff --git a/src/pint/toa.py b/src/pint/toa.py index 210a167d2..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,13 +1639,6 @@ 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) - - if len(valid_data) not in [0, len(self)]: - log.warning( - "Only a subset of the TOAs have wideband DM information. " - "A TOAs object will be treated as wideband only of *ALL* TOAs have DM information." - ) - return len(valid_data) == len(self) def get_all_flags(self): @@ -2993,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_toa.py b/tests/test_toa.py index b262e198b..7b367ba6d 100644 --- a/tests/test_toa.py +++ b/tests/test_toa.py @@ -184,15 +184,13 @@ def test_merge_toas(): assert np.all(toas_out.table == toas_outb.table) -def test_is_wideband(): - 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 - """ +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 + """ + ) ) - ) - assert not t1.is_wideband() - assert t1[0].is_wideband() - assert not t1[1].is_wideband() From d32e9954b30b260d742fffa28b1729e46adadd8c Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Thu, 31 Aug 2023 12:27:53 -0500 Subject: [PATCH 05/14] CHANGELOG --- CHANGELOG-unreleased.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index e998a4d5b..9c0ab243d 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -11,7 +11,6 @@ the released changes. ### Changed - 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()`). -- `TOAs.is_wideband` returns True only if *ALL* TOAs have the -pp_dm flag. Emits a warning if only some TOAs have this flag. ### Added ### Fixed - Fixed RTD by specifying theme explicitly. @@ -19,4 +18,5 @@ the released changes. - Setting `model.PARAM1 = model.PARAM2` no longer overrides the name of `PARAM1` - 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 From 0a4c316593102909af18de5d59e779ab04d84b18 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Thu, 31 Aug 2023 15:28:14 -0500 Subject: [PATCH 06/14] FIX TEST --- tests/test_wideband_dm_data.py | 46 +++++++++++++--------------------- 1 file changed, 18 insertions(+), 28 deletions(-) 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) From e91908d0a985d130c662e85cfd1bd80195b3e872 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Fri, 1 Sep 2023 16:03:14 -0500 Subject: [PATCH 07/14] initial version --- src/pint/models/timing_model.py | 78 +++++++++++++++++++++++++++++++++ 1 file changed, 78 insertions(+) diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index 1f02b8c75..ea511defc 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -799,6 +799,84 @@ 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 radial_position(self, barytimes): + """Return line-of-sight position 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 position + """ + # this should also updaate 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 ( + bbi.a1() * np.sin(psi) * (1 - bbi.ecc() ** 2) / (1 + bbi.ecc() * np.cos(nu)) + ) + + def radial_velocity(self, barytimes): + """Return line-of-sight velocity 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 + ----- + See [1]_ + + .. [1] Lorimer & Kramer, 2008, "The Handbook of Pulsar Astronomy", Eqn. 8.24 + """ + # this should also updaate 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 conjunction(self, baryMJD): """Return the time(s) of the first superior conjunction(s) after baryMJD. From 74506a917b92aa2ed8cd7104387ff98fac405765 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Fri, 1 Sep 2023 16:04:39 -0500 Subject: [PATCH 08/14] added some explanation --- src/pint/models/timing_model.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index ea511defc..9183a4baa 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -834,7 +834,7 @@ def radial_position(self, barytimes): ) def radial_velocity(self, barytimes): - """Return line-of-sight velocity at barycentric MJD times. + """Return line-of-sight velocity of the pulsar at barycentric MJD times. Parameters ---------- @@ -858,6 +858,9 @@ def radial_velocity(self, barytimes): Notes ----- + This is the radial velocity of the pulsar. For the radial velocity of the companion, + this must be multiplied by -1 times the mass of the pulsar divided by the mass of the companion. + See [1]_ .. [1] Lorimer & Kramer, 2008, "The Handbook of Pulsar Astronomy", Eqn. 8.24 From 25f38cafe8cd685b414b9741890e5448010e40ef Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Tue, 5 Sep 2023 09:16:25 -0500 Subject: [PATCH 09/14] Update CHANGELOG-unreleased.md --- CHANGELOG-unreleased.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index f0082c65d..726dfaef5 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -9,17 +9,17 @@ 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 `WaveX` model as a `DelayComponent` with Fourier amplitudes as fitted parameters ### 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 From 0cfaa5aa088f51ba2b92681a19a2041af2824ea3 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Tue, 5 Sep 2023 12:18:19 -0500 Subject: [PATCH 10/14] added tests --- tests/test_rv.py | 43 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) create mode 100644 tests/test_rv.py diff --git a/tests/test_rv.py b/tests/test_rv.py new file mode 100644 index 000000000..78867705e --- /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.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.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) From 7156eff57db4ee241304208f5de6526595db51e3 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Tue, 5 Sep 2023 12:18:57 -0500 Subject: [PATCH 11/14] changelog --- CHANGELOG-unreleased.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index f0082c65d..c53d67941 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -14,6 +14,7 @@ the released changes. - 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 method for binary models ### Fixed - Wave model `validate()` can correctly use PEPOCH to assign WAVEEPOCH parameter - Fixed RTD by specifying theme explicitly. From fd90fcc8ba8c48260bbe99b5b3ae5855b60de2d2 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Tue, 5 Sep 2023 14:32:50 -0500 Subject: [PATCH 12/14] added method for companion as well --- CHANGELOG-unreleased.md | 2 +- src/pint/models/timing_model.py | 43 ++++++++++++++++++++++++++++++--- tests/test_rv.py | 4 +-- 3 files changed, 42 insertions(+), 7 deletions(-) diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index 8ef61d444..f00179747 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -14,7 +14,7 @@ the released changes. - 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 method for binary models +- Added radial velocity methods for binary models ### Fixed - Wave model `validate()` can correctly use PEPOCH to assign WAVEEPOCH parameter - Fixed RTD by specifying theme explicitly. diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index 9183a4baa..f2b3dcff7 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -833,8 +833,8 @@ def radial_position(self, barytimes): bbi.a1() * np.sin(psi) * (1 - bbi.ecc() ** 2) / (1 + bbi.ecc() * np.cos(nu)) ) - def radial_velocity(self, barytimes): - """Return line-of-sight velocity of the pulsar at barycentric MJD times. + def pulsar_radial_velocity(self, barytimes): + """Return line-of-sight velocity of the pulsar relative to the system barycenter at barycentric MJD times. Parameters ---------- @@ -858,8 +858,7 @@ def radial_velocity(self, barytimes): Notes ----- - This is the radial velocity of the pulsar. For the radial velocity of the companion, - this must be multiplied by -1 times the mass of the pulsar divided by the mass of the companion. + This is the radial velocity of the pulsar. See [1]_ @@ -880,6 +879,42 @@ def radial_velocity(self, barytimes): * (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/tests/test_rv.py b/tests/test_rv.py index 78867705e..6ea9d4293 100644 --- a/tests/test_rv.py +++ b/tests/test_rv.py @@ -20,7 +20,7 @@ def setup_method(self): 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.radial_velocity(self.ts) + 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) @@ -38,6 +38,6 @@ def test_rv_basemodel(self): @pytest.mark.parametrize("othermodel", ["ELL1", "ELL1H", "DDS", "BT"]) def test_rv_othermodels(self, othermodel): mc = pint.binaryconvert.convert_binary(self.m, othermodel) - vc = mc.radial_velocity(self.ts) + 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) From a7416b4fd0befdf66a199627c5767d2d09f214cf Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Tue, 5 Sep 2023 14:34:22 -0500 Subject: [PATCH 13/14] removed untested function --- src/pint/models/timing_model.py | 34 --------------------------------- 1 file changed, 34 deletions(-) diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index f2b3dcff7..81943395f 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -799,40 +799,6 @@ 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 radial_position(self, barytimes): - """Return line-of-sight position 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 position - """ - # this should also updaate 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 ( - bbi.a1() * np.sin(psi) * (1 - bbi.ecc() ** 2) / (1 + bbi.ecc() * np.cos(nu)) - ) - def pulsar_radial_velocity(self, barytimes): """Return line-of-sight velocity of the pulsar relative to the system barycenter at barycentric MJD times. From cb8cdfdbe2a1e380ec983eec934023d5ab6fb2c1 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Thu, 7 Sep 2023 11:39:06 -0500 Subject: [PATCH 14/14] typo --- src/pint/models/timing_model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index 81943395f..6cc845412 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -830,7 +830,7 @@ def pulsar_radial_velocity(self, barytimes): .. [1] Lorimer & Kramer, 2008, "The Handbook of Pulsar Astronomy", Eqn. 8.24 """ - # this should also updaate the binary instance + # 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]