From e91908d0a985d130c662e85cfd1bd80195b3e872 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Fri, 1 Sep 2023 16:03:14 -0500 Subject: [PATCH 1/7] 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 2/7] 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 0cfaa5aa088f51ba2b92681a19a2041af2824ea3 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Tue, 5 Sep 2023 12:18:19 -0500 Subject: [PATCH 3/7] 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 4/7] 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 5/7] 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 6/7] 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 7/7] 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]