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 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 = ( 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()