From d767cfd853ecaf7d7d1135e3e576e75430300f31 Mon Sep 17 00:00:00 2001 From: Dotan Gazith Date: Sun, 25 Feb 2024 10:13:09 +0200 Subject: [PATCH 1/5] correcting ELL1 model references --- src/pint/models/binary_ell1.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/pint/models/binary_ell1.py b/src/pint/models/binary_ell1.py index f10c64765..bc62934fd 100644 --- a/src/pint/models/binary_ell1.py +++ b/src/pint/models/binary_ell1.py @@ -1,4 +1,5 @@ """Approximate binary model for small eccentricity.""" + import astropy.units as u import numpy as np from astropy.time import Time @@ -81,8 +82,8 @@ class BinaryELL1(PulsarBinary): - Zhu et al. (2019), MNRAS, 482 (3), 3249-3260 [2]_ - Fiore et al. (2023), arXiv:2305.13624 [astro-ph.HE] [3]_ - .. [1] https://ui.adsabs.harvard.edu/abs/2019MNRAS.482.3249Z/abstract - .. [2] https://ui.adsabs.harvard.edu/abs/2001MNRAS.326..274L/abstract + .. [1] https://ui.adsabs.harvard.edu/abs/2001MNRAS.326..274L/abstract + .. [2] https://ui.adsabs.harvard.edu/abs/2019MNRAS.482.3249Z/abstract .. [3] https://arxiv.org/abs/2305.13624 Notes From 7079f00bb4cf94ab7cd42515103b00f613cad8d4 Mon Sep 17 00:00:00 2001 From: Dotan Gazith Date: Tue, 5 Mar 2024 16:33:49 +0200 Subject: [PATCH 2/5] fixing units in binary delay derivatives when deriving with respect to parameters with unit "deg" now returns correct units --- .../stand_alone_psr_binaries/BT_model.py | 23 +++++++++++++------ .../stand_alone_psr_binaries/DD_model.py | 8 +++---- 2 files changed, 19 insertions(+), 12 deletions(-) diff --git a/src/pint/models/stand_alone_psr_binaries/BT_model.py b/src/pint/models/stand_alone_psr_binaries/BT_model.py index 2e8a121f1..0fba9fe71 100644 --- a/src/pint/models/stand_alone_psr_binaries/BT_model.py +++ b/src/pint/models/stand_alone_psr_binaries/BT_model.py @@ -143,17 +143,17 @@ def BTdelay(self): """Full BT model delay""" return (self.delayL1() + self.delayL2()) * self.delayR() - # NOTE: Below, OMEGA is supposed to be in RADIANS! - # TODO: Fix UNITS!!! def d_delayL1_d_E(self): a1 = self.a1() / c.c - return -a1 * np.sin(self.omega()) * np.sin(self.E()) + return -a1 * np.sin(self.omega()) * np.sin(self.E()) / self.E().unit def d_delayL2_d_E(self): a1 = self.a1() / c.c return ( - a1 * np.cos(self.omega()) * np.sqrt(1 - self.ecc() ** 2) + self.GAMMA - ) * np.cos(self.E()) + (a1 * np.cos(self.omega()) * np.sqrt(1 - self.ecc() ** 2) + self.GAMMA) + * np.cos(self.E()) + / self.E().unit + ) def d_delayL1_d_A1(self): return np.sin(self.omega()) * (np.cos(self.E()) - self.ecc()) / c.c @@ -171,7 +171,12 @@ def d_delayL2_d_A1DOT(self): def d_delayL1_d_OM(self): a1 = self.a1() / c.c - return a1 * np.cos(self.omega()) * (np.cos(self.E()) - self.ecc()) + return ( + a1 + * np.cos(self.omega()) + * (np.cos(self.E()) - self.ecc()) + / self.omega().unit + ) def d_delayL1_d_OMDOT(self): return self.tt0 * self.d_delayL1_d_OM() @@ -179,7 +184,11 @@ def d_delayL1_d_OMDOT(self): def d_delayL2_d_OM(self): a1 = self.a1() / c.c return ( - -a1 * np.sin(self.omega()) * np.sqrt(1 - self.ecc() ** 2) * np.sin(self.E()) + -a1 + * np.sin(self.omega()) + * np.sqrt(1 - self.ecc() ** 2) + * np.sin(self.E()) + / self.omega().unit ) def d_delayL2_d_OMDOT(self): diff --git a/src/pint/models/stand_alone_psr_binaries/DD_model.py b/src/pint/models/stand_alone_psr_binaries/DD_model.py index e9b82bdc3..2b6414f65 100644 --- a/src/pint/models/stand_alone_psr_binaries/DD_model.py +++ b/src/pint/models/stand_alone_psr_binaries/DD_model.py @@ -1,4 +1,5 @@ """Damour and Deruelle binary model.""" + import astropy.constants as c import astropy.units as u import numpy as np @@ -301,7 +302,7 @@ def d_beta_d_par(self, par): d_eTheta_d_par = self.d_eTheta_d_par(par) d_beta_d_a1 = 1.0 / c.c * (1 - eTheta**2) ** 0.5 * cosOmg - d_beta_d_omega = -a1 / c.c * (1 - eTheta**2) ** 0.5 * sinOmg + d_beta_d_omega = -a1 / c.c * (1 - eTheta**2) ** 0.5 * sinOmg / omega.unit d_beta_d_eTheta = a1 / c.c * (-eTheta) / np.sqrt(1 - eTheta**2) * cosOmg with u.set_enabled_equivalencies(u.dimensionless_angles()): return ( @@ -756,10 +757,7 @@ def d_delayS_d_par(self, par): -2 * TM2 / logNum - * ( - e * sE - - self.SINI * (np.sqrt(1 - e**2) * cE * cOmega - sE * sOmega) - ) + * (e * sE - self.SINI * (np.sqrt(1 - e**2) * cE * cOmega - sE * sOmega)) ) domega_dpar = self.prtl_der("omega", par) dsDelay_domega = ( From d531d64f4913beff3d816024c5a9d8a9061e3a2e Mon Sep 17 00:00:00 2001 From: Dotan Gazith Date: Wed, 6 Mar 2024 08:49:15 +0200 Subject: [PATCH 3/5] Revert "fixing units in binary delay derivatives" This reverts commit 7079f00bb4cf94ab7cd42515103b00f613cad8d4. --- .../stand_alone_psr_binaries/BT_model.py | 23 ++++++------------- .../stand_alone_psr_binaries/DD_model.py | 8 ++++--- 2 files changed, 12 insertions(+), 19 deletions(-) diff --git a/src/pint/models/stand_alone_psr_binaries/BT_model.py b/src/pint/models/stand_alone_psr_binaries/BT_model.py index 0fba9fe71..2e8a121f1 100644 --- a/src/pint/models/stand_alone_psr_binaries/BT_model.py +++ b/src/pint/models/stand_alone_psr_binaries/BT_model.py @@ -143,17 +143,17 @@ def BTdelay(self): """Full BT model delay""" return (self.delayL1() + self.delayL2()) * self.delayR() + # NOTE: Below, OMEGA is supposed to be in RADIANS! + # TODO: Fix UNITS!!! def d_delayL1_d_E(self): a1 = self.a1() / c.c - return -a1 * np.sin(self.omega()) * np.sin(self.E()) / self.E().unit + return -a1 * np.sin(self.omega()) * np.sin(self.E()) def d_delayL2_d_E(self): a1 = self.a1() / c.c return ( - (a1 * np.cos(self.omega()) * np.sqrt(1 - self.ecc() ** 2) + self.GAMMA) - * np.cos(self.E()) - / self.E().unit - ) + a1 * np.cos(self.omega()) * np.sqrt(1 - self.ecc() ** 2) + self.GAMMA + ) * np.cos(self.E()) def d_delayL1_d_A1(self): return np.sin(self.omega()) * (np.cos(self.E()) - self.ecc()) / c.c @@ -171,12 +171,7 @@ def d_delayL2_d_A1DOT(self): def d_delayL1_d_OM(self): a1 = self.a1() / c.c - return ( - a1 - * np.cos(self.omega()) - * (np.cos(self.E()) - self.ecc()) - / self.omega().unit - ) + return a1 * np.cos(self.omega()) * (np.cos(self.E()) - self.ecc()) def d_delayL1_d_OMDOT(self): return self.tt0 * self.d_delayL1_d_OM() @@ -184,11 +179,7 @@ def d_delayL1_d_OMDOT(self): def d_delayL2_d_OM(self): a1 = self.a1() / c.c return ( - -a1 - * np.sin(self.omega()) - * np.sqrt(1 - self.ecc() ** 2) - * np.sin(self.E()) - / self.omega().unit + -a1 * np.sin(self.omega()) * np.sqrt(1 - self.ecc() ** 2) * np.sin(self.E()) ) def d_delayL2_d_OMDOT(self): diff --git a/src/pint/models/stand_alone_psr_binaries/DD_model.py b/src/pint/models/stand_alone_psr_binaries/DD_model.py index 2b6414f65..e9b82bdc3 100644 --- a/src/pint/models/stand_alone_psr_binaries/DD_model.py +++ b/src/pint/models/stand_alone_psr_binaries/DD_model.py @@ -1,5 +1,4 @@ """Damour and Deruelle binary model.""" - import astropy.constants as c import astropy.units as u import numpy as np @@ -302,7 +301,7 @@ def d_beta_d_par(self, par): d_eTheta_d_par = self.d_eTheta_d_par(par) d_beta_d_a1 = 1.0 / c.c * (1 - eTheta**2) ** 0.5 * cosOmg - d_beta_d_omega = -a1 / c.c * (1 - eTheta**2) ** 0.5 * sinOmg / omega.unit + d_beta_d_omega = -a1 / c.c * (1 - eTheta**2) ** 0.5 * sinOmg d_beta_d_eTheta = a1 / c.c * (-eTheta) / np.sqrt(1 - eTheta**2) * cosOmg with u.set_enabled_equivalencies(u.dimensionless_angles()): return ( @@ -757,7 +756,10 @@ def d_delayS_d_par(self, par): -2 * TM2 / logNum - * (e * sE - self.SINI * (np.sqrt(1 - e**2) * cE * cOmega - sE * sOmega)) + * ( + e * sE + - self.SINI * (np.sqrt(1 - e**2) * cE * cOmega - sE * sOmega) + ) ) domega_dpar = self.prtl_der("omega", par) dsDelay_domega = ( From fa37e9f21afbf7dd5ed940cd2ac86a30c8b5cec8 Mon Sep 17 00:00:00 2001 From: Dotan Gazith Date: Sat, 9 Mar 2024 17:57:11 +0200 Subject: [PATCH 4/5] Fixing units of derivatives wrt 'rad' parameters in BT and DD models. And adding exact derivatives to BT model. --- .../stand_alone_psr_binaries/BT_model.py | 110 ++++++++++++++++-- .../stand_alone_psr_binaries/DD_model.py | 8 +- 2 files changed, 106 insertions(+), 12 deletions(-) diff --git a/src/pint/models/stand_alone_psr_binaries/BT_model.py b/src/pint/models/stand_alone_psr_binaries/BT_model.py index 2e8a121f1..8dc3da717 100644 --- a/src/pint/models/stand_alone_psr_binaries/BT_model.py +++ b/src/pint/models/stand_alone_psr_binaries/BT_model.py @@ -143,17 +143,17 @@ def BTdelay(self): """Full BT model delay""" return (self.delayL1() + self.delayL2()) * self.delayR() - # NOTE: Below, OMEGA is supposed to be in RADIANS! - # TODO: Fix UNITS!!! def d_delayL1_d_E(self): a1 = self.a1() / c.c - return -a1 * np.sin(self.omega()) * np.sin(self.E()) + return -a1 * np.sin(self.omega()) * np.sin(self.E()) / u.rad def d_delayL2_d_E(self): a1 = self.a1() / c.c return ( - a1 * np.cos(self.omega()) * np.sqrt(1 - self.ecc() ** 2) + self.GAMMA - ) * np.cos(self.E()) + (a1 * np.cos(self.omega()) * np.sqrt(1 - self.ecc() ** 2) + self.GAMMA) + * np.cos(self.E()) + / u.rad + ) def d_delayL1_d_A1(self): return np.sin(self.omega()) * (np.cos(self.E()) - self.ecc()) / c.c @@ -171,7 +171,7 @@ def d_delayL2_d_A1DOT(self): def d_delayL1_d_OM(self): a1 = self.a1() / c.c - return a1 * np.cos(self.omega()) * (np.cos(self.E()) - self.ecc()) + return a1 * np.cos(self.omega()) * (np.cos(self.E()) - self.ecc()) / u.rad def d_delayL1_d_OMDOT(self): return self.tt0 * self.d_delayL1_d_OM() @@ -179,7 +179,11 @@ def d_delayL1_d_OMDOT(self): def d_delayL2_d_OM(self): a1 = self.a1() / c.c return ( - -a1 * np.sin(self.omega()) * np.sqrt(1 - self.ecc() ** 2) * np.sin(self.E()) + -a1 + * np.sin(self.omega()) + * np.sqrt(1 - self.ecc() ** 2) + * np.sin(self.E()) + / u.rad ) def d_delayL2_d_OMDOT(self): @@ -207,6 +211,8 @@ def d_delayL1_d_GAMMA(self): def d_delayL2_d_GAMMA(self): return np.sin(self.E()) + # NOTE: The two derivatives below are taken care by lines 222, 237 + # TODO: Remove if indeed redundant def d_delayL1_d_T0(self): return self.d_delayL1_d_E() * self.d_E_d_T0() @@ -245,3 +251,93 @@ def d_delayL2_d_par(self, par): def d_BTdelay_d_par(self, par): return self.delayR() * (self.d_delayL1_d_par(par) + self.d_delayL2_d_par(par)) + + +class BTmodel_extended_derivatives(BTmodel): + """ + This class extends the BTmodel class to include the exact first derivatives and some second derivatives. + """ + + def __init__(self, t=None, input_params=None): + super().__init__(t, input_params) + + def d_delayR_d_E(self): + a1 = self.a1() / c.c + omega = self.omega() + ecc = self.ecc() + E = self.E() + num = np.sqrt(1 - ecc**2) * np.cos(omega) * np.sin(E) + ( + np.cos(E) - ecc + ) * np.sin(omega) + den = (1 - ecc * np.cos(E)) ** 2 + return 2 * np.pi * a1 * num / (den * self.pb().to(u.second)) + + def d_delayR_d_A1(self): + omega = self.omega() + ecc = self.ecc() + E = self.E() + num = np.cos(omega) * np.sqrt(1 - ecc**2) * np.cos(E) - np.sin( + omega + ) * np.sin(E) + den = 1.0 - ecc * np.cos(E) + return 1.0 - 2 * np.pi * num / (den * self.pb().to(u.second)) + + def d_delayR_d_A1DOT(self): + return self.tt0 * self.d_delayR_d_A1() + + def d_delayR_d_OM(self): + a1 = self.a1() / c.c + omega = self.omega() + ecc = self.ecc() + E = self.E() + num = -np.sin(omega) * np.sqrt(1 - ecc**2) * np.cos(E) - np.cos( + omega + ) * np.sin(E) + den = 1.0 - ecc * np.cos(E) + return -2 * np.pi * a1 * num / (den * self.pb().to(u.second)) + + def d_delayR_d_OMDOT(self): + return self.tt0 * self.d_delayR_d_OM() + + def d_delayR_d_ECC(self): + a1 = self.a1() / c.c + omega = self.omega() + ecc = self.ecc() + E = self.E() + num = (ecc - np.cos(E)) * np.cos(omega) + np.sqrt(1 - ecc**2) * np.sin( + E + ) * np.sin(omega) + den = np.sqrt(1 - ecc**2) * (1 - ecc * np.cos(E)) ** 2 + return ( + 2 * np.pi * a1 * np.cos(E) * num / (den * self.pb().to(u.second)) + + self.d_delayR_d_E() * self.d_E_d_ECC() + ) + + def d_delayR_d_EDOT(self): + return self.tt0 * self.d_delayR_d_ECC() + + def d_delayR_d_GAMMA(self): + return np.zeros(len(self.t)) * u.second / u.second + + def d_delayR_d_T0(self): + return self.d_delayR_d_E() * self.d_E_d_T0() + + def d_delayR_d_par(self, par): + if par not in self.binary_params: + errorMesg = f"{par} is not in binary parameter list." + raise ValueError(errorMesg) + + par_obj = getattr(self, par) + if not hasattr(self, f"d_delayR_d_{par}"): + return ( + self.d_delayR_d_E() * self.d_E_d_par(par) + if par in self.orbits_cls.orbit_params + else np.zeros(len(self.t)) * u.second / par_obj.unit + ) + func = getattr(self, f"d_delayR_d_{par}") + return func() + + def d_BTdelay_d_par(self, par): + return self.delayR() * ( + self.d_delayL1_d_par(par) + self.d_delayL2_d_par(par) + ) + (self.delayL1() + self.delayL2()) * self.d_delayR_d_par(par) diff --git a/src/pint/models/stand_alone_psr_binaries/DD_model.py b/src/pint/models/stand_alone_psr_binaries/DD_model.py index e9b82bdc3..b2116f1ad 100644 --- a/src/pint/models/stand_alone_psr_binaries/DD_model.py +++ b/src/pint/models/stand_alone_psr_binaries/DD_model.py @@ -1,4 +1,5 @@ """Damour and Deruelle binary model.""" + import astropy.constants as c import astropy.units as u import numpy as np @@ -301,7 +302,7 @@ def d_beta_d_par(self, par): d_eTheta_d_par = self.d_eTheta_d_par(par) d_beta_d_a1 = 1.0 / c.c * (1 - eTheta**2) ** 0.5 * cosOmg - d_beta_d_omega = -a1 / c.c * (1 - eTheta**2) ** 0.5 * sinOmg + d_beta_d_omega = -a1 / c.c * (1 - eTheta**2) ** 0.5 * sinOmg / u.rad d_beta_d_eTheta = a1 / c.c * (-eTheta) / np.sqrt(1 - eTheta**2) * cosOmg with u.set_enabled_equivalencies(u.dimensionless_angles()): return ( @@ -756,10 +757,7 @@ def d_delayS_d_par(self, par): -2 * TM2 / logNum - * ( - e * sE - - self.SINI * (np.sqrt(1 - e**2) * cE * cOmega - sE * sOmega) - ) + * (e * sE - self.SINI * (np.sqrt(1 - e**2) * cE * cOmega - sE * sOmega)) ) domega_dpar = self.prtl_der("omega", par) dsDelay_domega = ( From 8479a2d232c8c71d933b229009be94ce2added9e Mon Sep 17 00:00:00 2001 From: Dotan Gazith Date: Thu, 14 Mar 2024 19:41:32 +0200 Subject: [PATCH 5/5] fix units of generic d_E/d_T0 to be rad/s (instead of 1/s) --- src/pint/models/stand_alone_psr_binaries/binary_generic.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pint/models/stand_alone_psr_binaries/binary_generic.py b/src/pint/models/stand_alone_psr_binaries/binary_generic.py index c62f174a3..35ff04c7f 100644 --- a/src/pint/models/stand_alone_psr_binaries/binary_generic.py +++ b/src/pint/models/stand_alone_psr_binaries/binary_generic.py @@ -502,7 +502,7 @@ def d_E_d_T0(self): def d_E_d_ECC(self): E = self.E() - return np.sin(E) / (1.0 - self.ecc() * np.cos(E)) + return np.sin(E) / (1.0 - self.ecc() * np.cos(E)) * u.rad def d_E_d_EDOT(self): return self.tt0 * self.d_E_d_ECC()