From 6175ec7041f61183b466bce30602f7fde47f2674 Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Fri, 29 Sep 2023 08:53:05 +0200 Subject: [PATCH 01/14] First version of MockPulsar --- enterprise/pulsar.py | 59 +++++++++++++++++++++++++++++++++++++ enterprise/signals/utils.py | 39 ++++++++++++++++++++++++ 2 files changed, 98 insertions(+) diff --git a/enterprise/pulsar.py b/enterprise/pulsar.py index b8f23cd5..e44edc68 100644 --- a/enterprise/pulsar.py +++ b/enterprise/pulsar.py @@ -602,6 +602,65 @@ def destroy(psr): # pragma: py-lt-38 psr._deflated = "destroyed" + +class MockPulsar(BasePulsar): + """Class to allow mock pulsars to be used with Enterprise""" + # TODO: allow units? + # test utils.get_psrname_from_raj_decj + # utils.get_psrname_from_pos + # utils.create_spindown_timing_model + + def __init__(self, obs_times_mjd, elong=None, elat=None, raj=None, + decj=None, ssbfreqs=1440.0, residuals=None, toaerrs=1e-6, + sort=True, telescope='GBT', spindown_order=2): + + self.name = utils.get_psrname_from_pos( + elong=elong, elat=elat, raj=raj, decj=decj) + + if elong is not None and elat is not None: + ec = ephem.Ecliptic(elong*np.pi/180.0, elat*np.pi/180.0) + eq = ephem.Equatorial(ec, epoch=ephem.J2000) + raj, decj = np.double(eq.ra), np.double(eq.dec) + + self._raj, self._decj = raj, decj + + self._sort = sort + self.planets = False + + self._toas = np.double(obs_times_mjd) * 86400 + self._stoas = np.double(obs_times_mjd) * 86400 + self._residuals = residuals if residuals else np.zeros_like(obs_times_mjd) + self._toaerrs = np.ones_like(obs_times_mjd) * toaerrs + + self._designmatrix, self.fitpars = utils.create_spindown_timing_model( + self._toas, order=spindown_order=2) + self._ssbfreqs = np.ones_like(self._toas) * ssbfreqs / 1e6 + self._telescope = telescope + + # set parameters + self.setpars = [fp for fp in self.fitpars] + + flags = {} + + # new-style storage of flags as a numpy record array (previously, psr._flags = flags) + self._flags = np.zeros(len(self._toas), dtype=[(key, val.dtype) for key, val in flags.items()]) + for key, val in flags.items(): + self._flags[key] = val + + self._pdist = self._get_pdist() + + self._pos = self._get_pos() + self._planetssb = None + self._sunssb = None + + self._pos_t = np.tile(self._pos, (len(self._toas), 1)) + + self.sort_data() + + def set_residuals(self, residuals): + self._residuals = residuals + + def Pulsar(*args, **kwargs): ephem = kwargs.get("ephem", None) clk = kwargs.get("clk", None) diff --git a/enterprise/signals/utils.py b/enterprise/signals/utils.py index df949af0..a2b45ea9 100644 --- a/enterprise/signals/utils.py +++ b/enterprise/signals/utils.py @@ -14,6 +14,7 @@ from scipy.integrate import odeint from scipy.interpolate import interp1d from sksparse.cholmod import cholesky +from ephem import Ecliptic, Equatorial import enterprise from enterprise import constants as const @@ -31,6 +32,44 @@ logger = logging.getLogger(__name__) +def get_psrname_from_raj_decj(raj, decj): + """Get the pulsar name from RAJ and DECJ values""" + + raj_fhr = raj * 12 / np.pi + raj_hr = int(raj_fhr) + raj_min = int((raj_fhr-raj_hr)*60) + + decj_fdeg = decj * 180 / np.pi + decj_deg = int(decj_fdeg) + decj_min = int(np.abs((decj_fdeg-decj_deg)*60)) + + sign = "+" if np.sign(decj) > 0 else "-" + pos_str = f"J{raj_hr:02}{raj_min:02}{sign}{int(np.abs(decj_deg)):02}{decj_min:02}" + + return pos_str + +def get_psrname_from_pos(elong=None, elat=None, raj=None, decj=None): + """Get the pulsar name from position parameters""" + + if elong is not None and elat is not None: + ec = ephem.Ecliptic(elong*np.pi/180.0, elat*np.pi/180.0) + eq = ephem.Equatorial(ec, epoch=ephem.J2000) + raj, decj = float(eq.ra), float(eq.dec) + elif raj is None or decj is None: + raise ValueError("Need to provide either raj/decj or elong/elat!") + + return get_psrname_from_raj_decj(raj, decj) + +def create_spindown_timing_model(toas, order=2): + """Create a spindown-only timing model with some order (default 2 = QSD)""" + + avetoas = (toas-np.mean(toas)) / np.mean(toas) + designmatrix = np.vstack([avetoas**ii for ii in range(1+order)]).T + parameter_names = ['Offset'] + ['F{ii}'.format(ii) for ii in range(order)] + + return designmatrix, parameter_names + + class ConditionalGP: def __init__(self, pta, phiinv_method="cliques"): """This class allows the computation of conditional means and From 68181e25eb2fd2bf454f164ef97938163614d959 Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Fri, 29 Sep 2023 10:38:22 +0200 Subject: [PATCH 02/14] Small typo --- enterprise/pulsar.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/enterprise/pulsar.py b/enterprise/pulsar.py index e44edc68..cdee2ce4 100644 --- a/enterprise/pulsar.py +++ b/enterprise/pulsar.py @@ -633,7 +633,7 @@ def __init__(self, obs_times_mjd, elong=None, elat=None, raj=None, self._toaerrs = np.ones_like(obs_times_mjd) * toaerrs self._designmatrix, self.fitpars = utils.create_spindown_timing_model( - self._toas, order=spindown_order=2) + self._toas, order=spindown_order) self._ssbfreqs = np.ones_like(self._toas) * ssbfreqs / 1e6 self._telescope = telescope From 160f4762ee1c3b91735745e3ad09dfdd8fb695b9 Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Fri, 24 Nov 2023 11:17:33 +0100 Subject: [PATCH 03/14] Ansatz for getting a simple Astrometry model coded up for use in MockPulsar --- enterprise/signals/utils.py | 143 ++++++++++++++++++++++++++++++++++++ 1 file changed, 143 insertions(+) diff --git a/enterprise/signals/utils.py b/enterprise/signals/utils.py index a2b45ea9..c38b1024 100644 --- a/enterprise/signals/utils.py +++ b/enterprise/signals/utils.py @@ -69,6 +69,149 @@ def create_spindown_timing_model(toas, order=2): return designmatrix, parameter_names +def ssb_to_earth_vector(mjd_timestamps): + """Create the ssb_to_earth 3D vector""" + + from astropy.time import Time + from astropy.coordinates import get_body_barycentric, solar_system_ephemeris, ICRS, GeocentricTrueEcliptic + + # Set solar system ephemeris to 'builtin' for offline calculations + solar_system_ephemeris.set('builtin') + + # Convert MJD timestamps to Astropy Time objects + times = Time(mjd_timestamps, format='mjd') + + # Calculate Earth's position w.r.t. SSB at each timestamp + earth_positions_icrs = get_body_barycentric('earth', times) + + # Extract x, y, and z coordinates + ssb_obs = np.array([earth_positions_icrs.x.value, + earth_positions_icrs.y.value, + earth_positions_icrs.z.value]).T + + return ssb_obs + +def ssb_to_pulsar_vector(ra_radians, dec_radians, distance_parsecs): + """Create the ssb_to_pulsar 3D vector""" + + from astropy.coordinates import SkyCoord + import astropy.units as u + + # Create a SkyCoord object with the given RA, DEC, and distance + pulsar_coord = SkyCoord(ra=ra_radians*u.radian, + dec=dec_radians*u.radian, + distance=distance_parsecs*u.parsec, + frame='icrs') + + # Convert to Cartesian coordinates (x, y, z) + pulsar_vector = pulsar_coord.cartesian.xyz.value + + return pulsar_vector + +def d_delay_astrometry_d_PX(ssb_to_earth_vector, ssb_to_pulsar_vector): + """Calculate the derivative wrt PX + + Roughly following Smart, 1977, chapter 9. + + px_r: Extra distance to Earth, wrt SSB, from pulsar + r_e: Position of earth (vector) wrt SSB + u_p: Unit vector from SSB pointing to pulsar + t_d: Parallax delay + c: Speed of light + delta: Parallax + + The parallax delay is due to a distance orthogonal to the line of sight + to the pulsar from the SSB: + + px_r = sqrt( r_e**2 - (r_e.u_p)**2 ), + + with delay + + t_d = 0.5 * px_r * delta'/ c, and delta = delta' * px_r / (1 AU) + + """ + import astropy.units as u + import astropy.constants as const + + ssb_obs_r = np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) + in_psr_obs = np.sum(ssb_to_earth_vector * ssb_to_pulsar_vector, axis=1) + + px_r = np.sqrt(ssb_obs_r ** 2 - in_psr_obs ** 2) + dd_dpx = 0.5 * (px_r**2 / (u.AU * const.c)) * (u.mas / u.radian) + + return dd_dpx.decompose(u.si.bases) / u.mas + + def d_delay_astrometry_d_RAJ(ssb_to_earth_vector, ra_radians, dec_radians): + """Calculate the derivative wrt RAJ + + For the RAJ and DEC derivatives, use the following approximate model for + the pulse delay. (Inner-product between two Cartesian vectors): + + - de = Earth declination (wrt SSB) + - ae = Earth right ascension + - dp = pulsar declination + - aa = pulsar right ascension + - r = distance from SSB to Earh + - c = speed of light + + delay = r*[cos(de)*cos(dp)*cos(ae-aa)+sin(de)*sin(dp)]/c + """ + + earth_ra, earth_dec = np.arctan2(ssb_to_earth_vector[:,1], ssb_to_earth_vector[:,0]), np.arcsin(ssb_to_earth_vector[:,2] / np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1))) + + geom = ( + np.cos(earth_dec) * np.cos(dec_radians) * np.sin(ra_radians - earth_ra) + ) + dd_draj = np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) * geom / (const.c * u.radian) + + return dd_draj.decompose(u.si.bases) + +def d_delay_astrometry_d_DECJ(ssb_to_earth_vector, ra_radians, dec_radians): + """Calculate the derivative wrt DECJ + + Definitions as in d_delay_d_RAJ + """ + + earth_ra, earth_dec = np.arctan2(ssb_to_earth_vector[:,1], ssb_to_earth_vector[:,0]), np.arcsin(ssb_to_earth_vector[:,2] / np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1))) + + geom = np.cos(earth_dec) * np.sin(dec_radians) * np.cos(ra_radians - earth_ra) - np.sin(earth_dec) * np.cos(dec_radians) + dd_ddecj = np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) * geom / (const.c * u.radian) + + return dd_ddecj.decompose(u.si.bases) + +def d_delay_astrometry_d_PMRA(mjd_timestamps, ssb_to_earth_vector, ra_radians, posepoch_mjd): + """Calculate the derivative wrt PMRA + + Definitions as in d_delay_d_RAJ. Now we have a derivative in mas/yr for + the pulsar RA + """ + + earth_ra = np.arctan2(ssb_to_earth_vector[:,1], ssb_to_earth_vector[:,0]) + + te = Time(mjd_timestamps, format='mjd') - Time(posepoch_mjd, format='mjd') + geom = np.cos(np.arcsin(ssb_to_earth_vector[:,2] / np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)))) * np.sin(ra_radians - earth_ra) + + deriv = np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) * geom * te.jd / (const.c * u.radian) + dd_dpmra = deriv * u.mas / u.year + + return dd_dpmra.decompose(u.si.bases) / (u.mas / u.year) + +def d_delay_astrometry_d_PMDEC(mjd_timestamps, ssb_to_earth_vector, ra_radians, dec_radians, posepoch_mjd): + """Calculate the derivative wrt PMDEC + + Definitions as in d_delay_d_RAJ. Now we have a derivative in mas/yr for + the pulsar DEC + """ + + earth_ra, earth_dec = np.arctan2(ssb_to_earth_vector[:,1], ssb_to_earth_vector[:,0]), np.arcsin(ssb_to_earth_vector[:,2] / np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1))) + + te = Time(mjd_timestamps, format='mjd') - Time(posepoch_mjd, format='mjd') + geom = np.cos(earth_dec) * np.sin(dec_radians) * np.cos(ra_radians - earth_ra) - np.cos(dec_radians) * np.sin(earth_dec) + + deriv = np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) * geom * te.jd / (const.c * u.radian) + dd_dpmdec = deriv * u.mas / u.year + + return dd_dpmdec.decompose(u.si.bases) / (u.mas / u.year) class ConditionalGP: def __init__(self, pta, phiinv_method="cliques"): From 06ef3f8541fa969f4ae1eff25406058dbc8edc63 Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Fri, 24 Nov 2023 15:03:04 +0100 Subject: [PATCH 04/14] Added docstrings and helper functions for astrometry model --- enterprise/signals/utils.py | 192 ++++++++++++++++++++++++++---------- 1 file changed, 141 insertions(+), 51 deletions(-) diff --git a/enterprise/signals/utils.py b/enterprise/signals/utils.py index fe15d0c3..8440debb 100644 --- a/enterprise/signals/utils.py +++ b/enterprise/signals/utils.py @@ -33,8 +33,14 @@ def get_psrname_from_raj_decj(raj, decj): - """Get the pulsar name from RAJ and DECJ values""" + """ + Get the pulsar name from RAJ and DECJ values + + :param raj: Right ascension of the pulsar [rad] + :param decj: Declination of the pulsar [rad] + :returns: A string of the pulsar name in canonical format [str] + """ raj_fhr = raj * 12 / np.pi raj_hr = int(raj_fhr) raj_min = int((raj_fhr-raj_hr)*60) @@ -49,7 +55,16 @@ def get_psrname_from_raj_decj(raj, decj): return pos_str def get_psrname_from_pos(elong=None, elat=None, raj=None, decj=None): - """Get the pulsar name from position parameters""" + """ + Get the pulsar name from position parameters + + :param elong: Ecliptic longitude of the pulsar [rad] + :param elag: Ecliptic lattitute of the pulsar [rad] + :param raj: Right ascension of the pulsar [rad] + :param decj: Declination of the pulsar [rad] + + :returns: A string of the pulsar name in canonical format [str] + """ if elong is not None and elat is not None: ec = ephem.Ecliptic(elong*np.pi/180.0, elat*np.pi/180.0) @@ -60,17 +75,14 @@ def get_psrname_from_pos(elong=None, elat=None, raj=None, decj=None): return get_psrname_from_raj_decj(raj, decj) -def create_spindown_timing_model(toas, order=2): - """Create a spindown-only timing model with some order (default 2 = QSD)""" - - avetoas = (toas-np.mean(toas)) / np.mean(toas) - designmatrix = np.vstack([avetoas**ii for ii in range(1+order)]).T - parameter_names = ['Offset'] + ['F{ii}'.format(ii) for ii in range(order)] +def ssb_to_earth_vector(mjd_timestamps): + """ + Create the ssb_to_earth 3D vector - return designmatrix, parameter_names + :param mjd_timestamps: Timestamps [MJD] -def ssb_to_earth_vector(mjd_timestamps): - """Create the ssb_to_earth 3D vector""" + :returns: 2D array with the position of the Earth wrt the SSB + """ from astropy.time import Time from astropy.coordinates import get_body_barycentric, solar_system_ephemeris, ICRS, GeocentricTrueEcliptic @@ -92,7 +104,15 @@ def ssb_to_earth_vector(mjd_timestamps): return ssb_obs def ssb_to_pulsar_vector(ra_radians, dec_radians, distance_parsecs): - """Create the ssb_to_pulsar 3D vector""" + """ + Create the ssb_to_pulsar 3D vector + + :param ra_radians: Right ascension of the pulsar [rad] + :param dec_radians: Declination of the pulsar [rad] + :param distance_parsecs: Distance to the pulsar [pcs] + + :returns: 3D vector of the pulsar position in ICRS coordinates + """ from astropy.coordinates import SkyCoord import astropy.units as u @@ -109,26 +129,13 @@ def ssb_to_pulsar_vector(ra_radians, dec_radians, distance_parsecs): return pulsar_vector def d_delay_astrometry_d_PX(ssb_to_earth_vector, ssb_to_pulsar_vector): - """Calculate the derivative wrt PX - - Roughly following Smart, 1977, chapter 9. - - px_r: Extra distance to Earth, wrt SSB, from pulsar - r_e: Position of earth (vector) wrt SSB - u_p: Unit vector from SSB pointing to pulsar - t_d: Parallax delay - c: Speed of light - delta: Parallax - - The parallax delay is due to a distance orthogonal to the line of sight - to the pulsar from the SSB: - - px_r = sqrt( r_e**2 - (r_e.u_p)**2 ), - - with delay + """ + Calculate the derivative wrt PX - t_d = 0.5 * px_r * delta'/ c, and delta = delta' * px_r / (1 AU) + :param ssb_to_earth_vector: Position of Earth wrt SSB per timestamp + :param ssb_to_pulsar_vector: Direction of Pulsar wrt SSB + :returns: d_delay / d_PX """ import astropy.units as u import astropy.constants as const @@ -141,21 +148,18 @@ def d_delay_astrometry_d_PX(ssb_to_earth_vector, ssb_to_pulsar_vector): return dd_dpx.decompose(u.si.bases) / u.mas - def d_delay_astrometry_d_RAJ(ssb_to_earth_vector, ra_radians, dec_radians): - """Calculate the derivative wrt RAJ - - For the RAJ and DEC derivatives, use the following approximate model for - the pulse delay. (Inner-product between two Cartesian vectors): +def d_delay_astrometry_d_RAJ(ssb_to_earth_vector, ra_radians, dec_radians): + """ + Calculate the derivative wrt RAJ - - de = Earth declination (wrt SSB) - - ae = Earth right ascension - - dp = pulsar declination - - aa = pulsar right ascension - - r = distance from SSB to Earh - - c = speed of light + :param ssb_to_earth_vector: Position of Earth wrt SSB per timestamp + :param ra_radians: Right ascension of the Pulsar [rad] + :param dec_radians: Declination of the Pulsar [rad] - delay = r*[cos(de)*cos(dp)*cos(ae-aa)+sin(de)*sin(dp)]/c + :returns: d_delay / d_RAJ """ + import astropy.units as u + import astropy.constants as const earth_ra, earth_dec = np.arctan2(ssb_to_earth_vector[:,1], ssb_to_earth_vector[:,0]), np.arcsin(ssb_to_earth_vector[:,2] / np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1))) @@ -169,8 +173,14 @@ def d_delay_astrometry_d_RAJ(ssb_to_earth_vector, ra_radians, dec_radians): def d_delay_astrometry_d_DECJ(ssb_to_earth_vector, ra_radians, dec_radians): """Calculate the derivative wrt DECJ - Definitions as in d_delay_d_RAJ + :param ssb_to_earth_vector: Position of Earth wrt SSB per timestamp + :param ra_radians: Right ascension of the Pulsar + :param dec_radians: Declination of the Pulsar + + :returns: d_delay / d_DECJ """ + import astropy.units as u + import astropy.constants as const earth_ra, earth_dec = np.arctan2(ssb_to_earth_vector[:,1], ssb_to_earth_vector[:,0]), np.arcsin(ssb_to_earth_vector[:,2] / np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1))) @@ -182,16 +192,22 @@ def d_delay_astrometry_d_DECJ(ssb_to_earth_vector, ra_radians, dec_radians): def d_delay_astrometry_d_PMRA(mjd_timestamps, ssb_to_earth_vector, ra_radians, posepoch_mjd): """Calculate the derivative wrt PMRA - Definitions as in d_delay_d_RAJ. Now we have a derivative in mas/yr for - the pulsar RA + :param mjd_timestamps: Timestamps [MJD] + :param ssb_to_earth_vector: Position of Earth wrt SSB per timestamp [rad] + :param ra_radians: Right ascension of the Pulsar [rad] + :param posepoch_mjd: Position epoch [MJD] + + :returns: d_delay / d_PMRA """ + import astropy.units as u + import astropy.constants as const earth_ra = np.arctan2(ssb_to_earth_vector[:,1], ssb_to_earth_vector[:,0]) - te = Time(mjd_timestamps, format='mjd') - Time(posepoch_mjd, format='mjd') + te = (mjd_timestamps - posepoch_mjd) * u.day geom = np.cos(np.arcsin(ssb_to_earth_vector[:,2] / np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)))) * np.sin(ra_radians - earth_ra) - deriv = np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) * geom * te.jd / (const.c * u.radian) + deriv = np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) * geom * te / (const.c * u.radian) dd_dpmra = deriv * u.mas / u.year return dd_dpmra.decompose(u.si.bases) / (u.mas / u.year) @@ -199,20 +215,94 @@ def d_delay_astrometry_d_PMRA(mjd_timestamps, ssb_to_earth_vector, ra_radians, p def d_delay_astrometry_d_PMDEC(mjd_timestamps, ssb_to_earth_vector, ra_radians, dec_radians, posepoch_mjd): """Calculate the derivative wrt PMDEC - Definitions as in d_delay_d_RAJ. Now we have a derivative in mas/yr for - the pulsar DEC + :param mjd_timestamps: Timestamps [MJD] + :param ssb_to_earth_vector: Position of Earth wrt SSB per timestamp + :param ra_radians: Right ascension of the Pulsar [rad] + :param dec_radians: Declination of the Pulsar [rad] + :param posepoch_mjd: Position epoch [MJD] + + :returns: d_delay / d_PMDEC """ + import astropy.units as u + import astropy.constants as const earth_ra, earth_dec = np.arctan2(ssb_to_earth_vector[:,1], ssb_to_earth_vector[:,0]), np.arcsin(ssb_to_earth_vector[:,2] / np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1))) - te = Time(mjd_timestamps, format='mjd') - Time(posepoch_mjd, format='mjd') + te = (mjd_timestamps - posepoch_mjd) * u.day geom = np.cos(earth_dec) * np.sin(dec_radians) * np.cos(ra_radians - earth_ra) - np.cos(dec_radians) * np.sin(earth_dec) - deriv = np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) * geom * te.jd / (const.c * u.radian) + deriv = np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) * geom * te / (const.c * u.radian) dd_dpmdec = deriv * u.mas / u.year return dd_dpmdec.decompose(u.si.bases) / (u.mas / u.year) +def create_astrometry_timing_model(toas, raj, decj, posepoch): + """ + Create an astrometry-only timing model for pulsar (raj, decj) + + :param toas: The TOAs of the pulsar [sec] + :param raj: Right ascension of the pulsar [rad] + :param decj: Declination of the pulsar [rad] + :param posepoch: Position epoch [sec] + + :returns: design matrix, parameter names + """ + + toas_mjds = toas * 86400.0 + posepoch_mjd = posepoch * 86400.0 + + parameter_names = ['RAJ', 'DECJ', 'PMRA', 'PMDEC', 'PX'] + designmatrix = np.zeros((len(toas_mjds), 5)) + + ssb_to_earth = ssb_to_earth_vector(toas_mjds) + ssb_to_pulsar = ssb_to_pulsar_vector(raj, decj, 1.0) + ssb_to_pulsar_norm = ssb_to_pulsar / np.linalg.norm(ssb_to_pulsar) + + designmatrix[:,0] = d_delay_astrometry_d_RAJ(ssb_to_earth, raj, decj).value + designmatrix[:,1] = d_delay_astrometry_d_DECJ(ssb_to_earth, raj, decj).value + designmatrix[:,2] = d_delay_astrometry_d_PMRA(toas_mjds, ssb_to_earth, raj, posepoch_mjd).value + designmatrix[:,3] = d_delay_astrometry_d_PMDEC(toas_mjds, ssb_to_earth, raj, decj, posepoch_mjd).value + designmatrix[:,4] = d_delay_astrometry_d_PX(ssb_to_earth, ssb_to_pulsar_norm).value + + return designmatrix, parameter_names + +def create_spindown_timing_model(toas, order=2): + """ + Create a spindown-only timing model with some order (default 2 = QSD) + + :param toas: The TOAs of the pulsar [sec] + :param order: Maximum order of spindown [default 2] + + :returns: design matrix, parameter names + """ + + avetoas = (toas-np.mean(toas)) / np.mean(toas) + designmatrix = np.vstack([avetoas**ii for ii in range(1+order)]).T + parameter_names = ['Offset'] + ['F{ii}'.format(ii) for ii in range(order)] + + return designmatrix, parameter_names + +def create_astrometry_spin_timing_model(toas, raj, decj, posepoch, spindown_order=2): + """ + Create a spindown + astrometry timing model + + :param toas: The TOAs of the pulsar [sec] + :param raj: Right ascension of the pulsar [rad] + :param decj: Declination of the pulsar [rad] + :param posepoch: Position epoch [sec] + :param spindown_order: Maximum order of spindown [default 2] + + :returns: design matrix, parameter names + """ + + designmatrix_qsd, parameter_names_qsd = create_spindown_timing_model(toas, order=spindown_order) + designmatrix_astro, parameter_names_astro = create_astrometry_timing_model(toas, raj, decj, posepoch) + + parameter_names = parameter_names_qsd + parameter_names_astro + designmatrix = np.hstack([designmatrix_qsd, designmatrix_astro]) + + return designmatrix, parameter_names + class ConditionalGP: def __init__(self, pta, phiinv_method="cliques"): """This class allows the computation of conditional means and From 3c299122d1549718180b0ff4f53ac81a960a09b0 Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Fri, 24 Nov 2023 15:19:16 +0100 Subject: [PATCH 05/14] Linting, and added astrometry option to MockPulsar --- enterprise/pulsar.py | 38 +++++++----- enterprise/signals/utils.py | 111 ++++++++++++++++++++---------------- 2 files changed, 87 insertions(+), 62 deletions(-) diff --git a/enterprise/pulsar.py b/enterprise/pulsar.py index cdee2ce4..bfc1c1b2 100644 --- a/enterprise/pulsar.py +++ b/enterprise/pulsar.py @@ -311,7 +311,6 @@ def telescope(self): class PintPulsar(BasePulsar): def __init__(self, toas, model, sort=True, drop_pintpsr=True, planets=True): - self._sort = sort self.planets = planets self.name = model.PSR.value @@ -448,7 +447,6 @@ def _get_sunssb(self, toas, model): class Tempo2Pulsar(BasePulsar): def __init__(self, t2pulsar, sort=True, drop_t2pulsar=True, planets=True): - self._sort = sort self.t2pulsar = t2pulsar self.planets = planets @@ -602,23 +600,33 @@ def destroy(psr): # pragma: py-lt-38 psr._deflated = "destroyed" - class MockPulsar(BasePulsar): """Class to allow mock pulsars to be used with Enterprise""" + # TODO: allow units? # test utils.get_psrname_from_raj_decj # utils.get_psrname_from_pos # utils.create_spindown_timing_model - def __init__(self, obs_times_mjd, elong=None, elat=None, raj=None, - decj=None, ssbfreqs=1440.0, residuals=None, toaerrs=1e-6, - sort=True, telescope='GBT', spindown_order=2): - - self.name = utils.get_psrname_from_pos( - elong=elong, elat=elat, raj=raj, decj=decj) + def __init__( + self, + obs_times_mjd, + elong=None, + elat=None, + raj=None, + decj=None, + ssbfreqs=1440.0, + residuals=None, + toaerrs=1e-6, + sort=True, + telescope="GBT", + spindown_order=2, + inc_astrometry=True, + ): + self.name = utils.get_psrname_from_pos(elong=elong, elat=elat, raj=raj, decj=decj) if elong is not None and elat is not None: - ec = ephem.Ecliptic(elong*np.pi/180.0, elat*np.pi/180.0) + ec = ephem.Ecliptic(elong * np.pi / 180.0, elat * np.pi / 180.0) eq = ephem.Equatorial(ec, epoch=ephem.J2000) raj, decj = np.double(eq.ra), np.double(eq.dec) @@ -631,9 +639,14 @@ def __init__(self, obs_times_mjd, elong=None, elat=None, raj=None, self._stoas = np.double(obs_times_mjd) * 86400 self._residuals = residuals if residuals else np.zeros_like(obs_times_mjd) self._toaerrs = np.ones_like(obs_times_mjd) * toaerrs + self._posepoch = np.mean(self._toas) - self._designmatrix, self.fitpars = utils.create_spindown_timing_model( - self._toas, order=spindown_order) + if inc_astrometry: + self._designmatrix, self.fitpars = utils.create_astrometry_spin_timing_model( + self._toas, self._raj, self._decj, self._posepoch, spindown_order=spindown_order + ) + else: + self._designmatrix, self.fitpars = utils.create_spindown_timing_model(self._toas, order=spindown_order) self._ssbfreqs = np.ones_like(self._toas) * ssbfreqs / 1e6 self._telescope = telescope @@ -714,7 +727,6 @@ def Pulsar(*args, **kwargs): return PintPulsar(toas, model, sort=sort, drop_pintpsr=drop_pintpsr, planets=planets) elif timing_package.lower() == "tempo2": - # hack to set maxobs maxobs = get_maxobs(reltimfile) + 100 t2pulsar = t2.tempopulsar(relparfile, reltimfile, maxobs=maxobs, ephem=ephem, clk=clk) diff --git a/enterprise/signals/utils.py b/enterprise/signals/utils.py index 8440debb..d2478891 100644 --- a/enterprise/signals/utils.py +++ b/enterprise/signals/utils.py @@ -35,7 +35,7 @@ def get_psrname_from_raj_decj(raj, decj): """ Get the pulsar name from RAJ and DECJ values - + :param raj: Right ascension of the pulsar [rad] :param decj: Declination of the pulsar [rad] @@ -43,21 +43,22 @@ def get_psrname_from_raj_decj(raj, decj): """ raj_fhr = raj * 12 / np.pi raj_hr = int(raj_fhr) - raj_min = int((raj_fhr-raj_hr)*60) + raj_min = int((raj_fhr - raj_hr) * 60) decj_fdeg = decj * 180 / np.pi decj_deg = int(decj_fdeg) - decj_min = int(np.abs((decj_fdeg-decj_deg)*60)) + decj_min = int(np.abs((decj_fdeg - decj_deg) * 60)) sign = "+" if np.sign(decj) > 0 else "-" pos_str = f"J{raj_hr:02}{raj_min:02}{sign}{int(np.abs(decj_deg)):02}{decj_min:02}" return pos_str + def get_psrname_from_pos(elong=None, elat=None, raj=None, decj=None): """ Get the pulsar name from position parameters - + :param elong: Ecliptic longitude of the pulsar [rad] :param elag: Ecliptic lattitute of the pulsar [rad] :param raj: Right ascension of the pulsar [rad] @@ -67,7 +68,7 @@ def get_psrname_from_pos(elong=None, elat=None, raj=None, decj=None): """ if elong is not None and elat is not None: - ec = ephem.Ecliptic(elong*np.pi/180.0, elat*np.pi/180.0) + ec = ephem.Ecliptic(elong * np.pi / 180.0, elat * np.pi / 180.0) eq = ephem.Equatorial(ec, epoch=ephem.J2000) raj, decj = float(eq.ra), float(eq.dec) elif raj is None or decj is None: @@ -75,6 +76,7 @@ def get_psrname_from_pos(elong=None, elat=None, raj=None, decj=None): return get_psrname_from_raj_decj(raj, decj) + def ssb_to_earth_vector(mjd_timestamps): """ Create the ssb_to_earth 3D vector @@ -88,21 +90,20 @@ def ssb_to_earth_vector(mjd_timestamps): from astropy.coordinates import get_body_barycentric, solar_system_ephemeris, ICRS, GeocentricTrueEcliptic # Set solar system ephemeris to 'builtin' for offline calculations - solar_system_ephemeris.set('builtin') + solar_system_ephemeris.set("builtin") # Convert MJD timestamps to Astropy Time objects - times = Time(mjd_timestamps, format='mjd') + times = Time(mjd_timestamps, format="mjd") # Calculate Earth's position w.r.t. SSB at each timestamp - earth_positions_icrs = get_body_barycentric('earth', times) + earth_positions_icrs = get_body_barycentric("earth", times) # Extract x, y, and z coordinates - ssb_obs = np.array([earth_positions_icrs.x.value, - earth_positions_icrs.y.value, - earth_positions_icrs.z.value]).T + ssb_obs = np.array([earth_positions_icrs.x.value, earth_positions_icrs.y.value, earth_positions_icrs.z.value]).T return ssb_obs + def ssb_to_pulsar_vector(ra_radians, dec_radians, distance_parsecs): """ Create the ssb_to_pulsar 3D vector @@ -118,16 +119,16 @@ def ssb_to_pulsar_vector(ra_radians, dec_radians, distance_parsecs): import astropy.units as u # Create a SkyCoord object with the given RA, DEC, and distance - pulsar_coord = SkyCoord(ra=ra_radians*u.radian, - dec=dec_radians*u.radian, - distance=distance_parsecs*u.parsec, - frame='icrs') + pulsar_coord = SkyCoord( + ra=ra_radians * u.radian, dec=dec_radians * u.radian, distance=distance_parsecs * u.parsec, frame="icrs" + ) # Convert to Cartesian coordinates (x, y, z) pulsar_vector = pulsar_coord.cartesian.xyz.value return pulsar_vector + def d_delay_astrometry_d_PX(ssb_to_earth_vector, ssb_to_pulsar_vector): """ Calculate the derivative wrt PX @@ -139,15 +140,16 @@ def d_delay_astrometry_d_PX(ssb_to_earth_vector, ssb_to_pulsar_vector): """ import astropy.units as u import astropy.constants as const - + ssb_obs_r = np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) in_psr_obs = np.sum(ssb_to_earth_vector * ssb_to_pulsar_vector, axis=1) - px_r = np.sqrt(ssb_obs_r ** 2 - in_psr_obs ** 2) + px_r = np.sqrt(ssb_obs_r**2 - in_psr_obs**2) dd_dpx = 0.5 * (px_r**2 / (u.AU * const.c)) * (u.mas / u.radian) return dd_dpx.decompose(u.si.bases) / u.mas + def d_delay_astrometry_d_RAJ(ssb_to_earth_vector, ra_radians, dec_radians): """ Calculate the derivative wrt RAJ @@ -160,16 +162,17 @@ def d_delay_astrometry_d_RAJ(ssb_to_earth_vector, ra_radians, dec_radians): """ import astropy.units as u import astropy.constants as const - - earth_ra, earth_dec = np.arctan2(ssb_to_earth_vector[:,1], ssb_to_earth_vector[:,0]), np.arcsin(ssb_to_earth_vector[:,2] / np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1))) - geom = ( - np.cos(earth_dec) * np.cos(dec_radians) * np.sin(ra_radians - earth_ra) + earth_ra, earth_dec = np.arctan2(ssb_to_earth_vector[:, 1], ssb_to_earth_vector[:, 0]), np.arcsin( + ssb_to_earth_vector[:, 2] / np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) ) + + geom = np.cos(earth_dec) * np.cos(dec_radians) * np.sin(ra_radians - earth_ra) dd_draj = np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) * geom / (const.c * u.radian) return dd_draj.decompose(u.si.bases) + def d_delay_astrometry_d_DECJ(ssb_to_earth_vector, ra_radians, dec_radians): """Calculate the derivative wrt DECJ @@ -181,14 +184,19 @@ def d_delay_astrometry_d_DECJ(ssb_to_earth_vector, ra_radians, dec_radians): """ import astropy.units as u import astropy.constants as const - - earth_ra, earth_dec = np.arctan2(ssb_to_earth_vector[:,1], ssb_to_earth_vector[:,0]), np.arcsin(ssb_to_earth_vector[:,2] / np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1))) - geom = np.cos(earth_dec) * np.sin(dec_radians) * np.cos(ra_radians - earth_ra) - np.sin(earth_dec) * np.cos(dec_radians) + earth_ra, earth_dec = np.arctan2(ssb_to_earth_vector[:, 1], ssb_to_earth_vector[:, 0]), np.arcsin( + ssb_to_earth_vector[:, 2] / np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) + ) + + geom = np.cos(earth_dec) * np.sin(dec_radians) * np.cos(ra_radians - earth_ra) - np.sin(earth_dec) * np.cos( + dec_radians + ) dd_ddecj = np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) * geom / (const.c * u.radian) return dd_ddecj.decompose(u.si.bases) + def d_delay_astrometry_d_PMRA(mjd_timestamps, ssb_to_earth_vector, ra_radians, posepoch_mjd): """Calculate the derivative wrt PMRA @@ -201,17 +209,20 @@ def d_delay_astrometry_d_PMRA(mjd_timestamps, ssb_to_earth_vector, ra_radians, p """ import astropy.units as u import astropy.constants as const - - earth_ra = np.arctan2(ssb_to_earth_vector[:,1], ssb_to_earth_vector[:,0]) - + + earth_ra = np.arctan2(ssb_to_earth_vector[:, 1], ssb_to_earth_vector[:, 0]) + te = (mjd_timestamps - posepoch_mjd) * u.day - geom = np.cos(np.arcsin(ssb_to_earth_vector[:,2] / np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)))) * np.sin(ra_radians - earth_ra) + geom = np.cos(np.arcsin(ssb_to_earth_vector[:, 2] / np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)))) * np.sin( + ra_radians - earth_ra + ) deriv = np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) * geom * te / (const.c * u.radian) dd_dpmra = deriv * u.mas / u.year return dd_dpmra.decompose(u.si.bases) / (u.mas / u.year) + def d_delay_astrometry_d_PMDEC(mjd_timestamps, ssb_to_earth_vector, ra_radians, dec_radians, posepoch_mjd): """Calculate the derivative wrt PMDEC @@ -225,21 +236,26 @@ def d_delay_astrometry_d_PMDEC(mjd_timestamps, ssb_to_earth_vector, ra_radians, """ import astropy.units as u import astropy.constants as const - - earth_ra, earth_dec = np.arctan2(ssb_to_earth_vector[:,1], ssb_to_earth_vector[:,0]), np.arcsin(ssb_to_earth_vector[:,2] / np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1))) - + + earth_ra, earth_dec = np.arctan2(ssb_to_earth_vector[:, 1], ssb_to_earth_vector[:, 0]), np.arcsin( + ssb_to_earth_vector[:, 2] / np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) + ) + te = (mjd_timestamps - posepoch_mjd) * u.day - geom = np.cos(earth_dec) * np.sin(dec_radians) * np.cos(ra_radians - earth_ra) - np.cos(dec_radians) * np.sin(earth_dec) + geom = np.cos(earth_dec) * np.sin(dec_radians) * np.cos(ra_radians - earth_ra) - np.cos(dec_radians) * np.sin( + earth_dec + ) deriv = np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) * geom * te / (const.c * u.radian) dd_dpmdec = deriv * u.mas / u.year return dd_dpmdec.decompose(u.si.bases) / (u.mas / u.year) + def create_astrometry_timing_model(toas, raj, decj, posepoch): """ Create an astrometry-only timing model for pulsar (raj, decj) - + :param toas: The TOAs of the pulsar [sec] :param raj: Right ascension of the pulsar [rad] :param decj: Declination of the pulsar [rad] @@ -251,21 +267,22 @@ def create_astrometry_timing_model(toas, raj, decj, posepoch): toas_mjds = toas * 86400.0 posepoch_mjd = posepoch * 86400.0 - parameter_names = ['RAJ', 'DECJ', 'PMRA', 'PMDEC', 'PX'] + parameter_names = ["RAJ", "DECJ", "PMRA", "PMDEC", "PX"] designmatrix = np.zeros((len(toas_mjds), 5)) ssb_to_earth = ssb_to_earth_vector(toas_mjds) ssb_to_pulsar = ssb_to_pulsar_vector(raj, decj, 1.0) ssb_to_pulsar_norm = ssb_to_pulsar / np.linalg.norm(ssb_to_pulsar) - designmatrix[:,0] = d_delay_astrometry_d_RAJ(ssb_to_earth, raj, decj).value - designmatrix[:,1] = d_delay_astrometry_d_DECJ(ssb_to_earth, raj, decj).value - designmatrix[:,2] = d_delay_astrometry_d_PMRA(toas_mjds, ssb_to_earth, raj, posepoch_mjd).value - designmatrix[:,3] = d_delay_astrometry_d_PMDEC(toas_mjds, ssb_to_earth, raj, decj, posepoch_mjd).value - designmatrix[:,4] = d_delay_astrometry_d_PX(ssb_to_earth, ssb_to_pulsar_norm).value + designmatrix[:, 0] = d_delay_astrometry_d_RAJ(ssb_to_earth, raj, decj).value + designmatrix[:, 1] = d_delay_astrometry_d_DECJ(ssb_to_earth, raj, decj).value + designmatrix[:, 2] = d_delay_astrometry_d_PMRA(toas_mjds, ssb_to_earth, raj, posepoch_mjd).value + designmatrix[:, 3] = d_delay_astrometry_d_PMDEC(toas_mjds, ssb_to_earth, raj, decj, posepoch_mjd).value + designmatrix[:, 4] = d_delay_astrometry_d_PX(ssb_to_earth, ssb_to_pulsar_norm).value return designmatrix, parameter_names + def create_spindown_timing_model(toas, order=2): """ Create a spindown-only timing model with some order (default 2 = QSD) @@ -276,16 +293,17 @@ def create_spindown_timing_model(toas, order=2): :returns: design matrix, parameter names """ - avetoas = (toas-np.mean(toas)) / np.mean(toas) - designmatrix = np.vstack([avetoas**ii for ii in range(1+order)]).T - parameter_names = ['Offset'] + ['F{ii}'.format(ii) for ii in range(order)] + avetoas = (toas - np.mean(toas)) / np.mean(toas) + designmatrix = np.vstack([avetoas**ii for ii in range(1 + order)]).T + parameter_names = ["Offset"] + ["F{ii}".format(ii) for ii in range(order)] return designmatrix, parameter_names + def create_astrometry_spin_timing_model(toas, raj, decj, posepoch, spindown_order=2): """ Create a spindown + astrometry timing model - + :param toas: The TOAs of the pulsar [sec] :param raj: Right ascension of the pulsar [rad] :param decj: Declination of the pulsar [rad] @@ -303,6 +321,7 @@ def create_astrometry_spin_timing_model(toas, raj, decj, posepoch, spindown_orde return designmatrix, parameter_names + class ConditionalGP: def __init__(self, pta, phiinv_method="cliques"): """This class allows the computation of conditional means and @@ -576,12 +595,10 @@ def create_stabletimingdesignmatrix(designmat, fastDesign=True): Mm = designmat.copy() if fastDesign: - norm = np.sqrt(np.sum(Mm**2, axis=0)) Mm /= norm else: - u, s, v = np.linalg.svd(Mm) Mm = u[:, : len(s)] @@ -594,7 +611,6 @@ def create_stabletimingdesignmatrix(designmat, fastDesign=True): def make_ecc_interpolant(): - """ Make interpolation function from eccentricity file to determine number of harmonics to use for a given @@ -611,7 +627,6 @@ def make_ecc_interpolant(): def get_edot(F, mc, e): - """ Compute eccentricity derivative from Taylor et al. (2016) @@ -1321,7 +1336,6 @@ def createfourierdesignmatrix_physicalephem( "jup_orb_elements", "sat_orb_elements", ]: - ppar = locals()[parname] if ppar: if parname not in ["jup_orb_elements", "sat_orb_elements"]: @@ -1361,7 +1375,6 @@ def physical_ephem_delay( sat_orbit=None, equatorial=True, ): - # convert toas to MJD mjd = toas / 86400 From 782a706eec0d8aed45deaa7fbf990f553efe092e Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Mon, 27 Nov 2023 10:15:14 +0100 Subject: [PATCH 06/14] Reformatting and adding optional import check of astropy --- enterprise/signals/utils.py | 147 ++++++++++++++++++++---------------- 1 file changed, 84 insertions(+), 63 deletions(-) diff --git a/enterprise/signals/utils.py b/enterprise/signals/utils.py index d2478891..600857de 100644 --- a/enterprise/signals/utils.py +++ b/enterprise/signals/utils.py @@ -86,8 +86,12 @@ def ssb_to_earth_vector(mjd_timestamps): :returns: 2D array with the position of the Earth wrt the SSB """ - from astropy.time import Time - from astropy.coordinates import get_body_barycentric, solar_system_ephemeris, ICRS, GeocentricTrueEcliptic + try: + from astropy.time import Time + from astropy.coordinates import get_body_barycentric, solar_system_ephemeris, ICRS, GeocentricTrueEcliptic + except ImportError: # pragma: no cover + log.error("Astropy required for native astrometry timing models") + raise # Set solar system ephemeris to 'builtin' for offline calculations solar_system_ephemeris.set("builtin") @@ -99,9 +103,7 @@ def ssb_to_earth_vector(mjd_timestamps): earth_positions_icrs = get_body_barycentric("earth", times) # Extract x, y, and z coordinates - ssb_obs = np.array([earth_positions_icrs.x.value, earth_positions_icrs.y.value, earth_positions_icrs.z.value]).T - - return ssb_obs + return np.array([earth_positions_icrs.x.value, earth_positions_icrs.y.value, earth_positions_icrs.z.value]).T def ssb_to_pulsar_vector(ra_radians, dec_radians, distance_parsecs): @@ -115,8 +117,12 @@ def ssb_to_pulsar_vector(ra_radians, dec_radians, distance_parsecs): :returns: 3D vector of the pulsar position in ICRS coordinates """ - from astropy.coordinates import SkyCoord - import astropy.units as u + try: + from astropy.coordinates import SkyCoord + import astropy.units as u + except ImportError: # pragma: no cover + log.error("Astropy required for native astrometry timing models") + raise # Create a SkyCoord object with the given RA, DEC, and distance pulsar_coord = SkyCoord( @@ -124,33 +130,35 @@ def ssb_to_pulsar_vector(ra_radians, dec_radians, distance_parsecs): ) # Convert to Cartesian coordinates (x, y, z) - pulsar_vector = pulsar_coord.cartesian.xyz.value - - return pulsar_vector + return pulsar_coord.cartesian.xyz.value -def d_delay_astrometry_d_PX(ssb_to_earth_vector, ssb_to_pulsar_vector): +def d_delay_d_PX(ssb_to_earth_vector, ssb_to_pulsar_vector): """ Calculate the derivative wrt PX :param ssb_to_earth_vector: Position of Earth wrt SSB per timestamp :param ssb_to_pulsar_vector: Direction of Pulsar wrt SSB - :returns: d_delay / d_PX + :returns: d_delay / d_PX [sec / mas] """ - import astropy.units as u - import astropy.constants as const + try: + import astropy.units as u + import astropy.constants as ac + except ImportError: # pragma: no cover + log.error("Astropy required for native astrometry timing models") + raise - ssb_obs_r = np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) - in_psr_obs = np.sum(ssb_to_earth_vector * ssb_to_pulsar_vector, axis=1) + ssb_earh_r = np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) + in_prod = np.sum(ssb_to_earth_vector * ssb_to_pulsar_vector, axis=1) - px_r = np.sqrt(ssb_obs_r**2 - in_psr_obs**2) - dd_dpx = 0.5 * (px_r**2 / (u.AU * const.c)) * (u.mas / u.radian) + px_r = np.sqrt(ssb_earh_r**2 - in_prod**2) + dd_dpx = 0.5 * (px_r**2 / (u.AU * ac.c)) * (u.mas / u.radian) return dd_dpx.decompose(u.si.bases) / u.mas -def d_delay_astrometry_d_RAJ(ssb_to_earth_vector, ra_radians, dec_radians): +def d_delay_d_RAJ(ssb_to_earth_vector, ra_radians, dec_radians): """ Calculate the derivative wrt RAJ @@ -160,20 +168,23 @@ def d_delay_astrometry_d_RAJ(ssb_to_earth_vector, ra_radians, dec_radians): :returns: d_delay / d_RAJ """ - import astropy.units as u - import astropy.constants as const + try: + import astropy.units as u + import astropy.constants as ac + except ImportError: # pragma: no cover + log.error("Astropy required for native astrometry timing models") + raise - earth_ra, earth_dec = np.arctan2(ssb_to_earth_vector[:, 1], ssb_to_earth_vector[:, 0]), np.arcsin( - ssb_to_earth_vector[:, 2] / np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) - ) + earth_ra = np.arctan2(ssb_to_earth_vector[:, 1], ssb_to_earth_vector[:, 0]) + earth_dec = np.arcsin(ssb_to_earth_vector[:, 2] / np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1))) - geom = np.cos(earth_dec) * np.cos(dec_radians) * np.sin(ra_radians - earth_ra) - dd_draj = np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) * geom / (const.c * u.radian) + geometric = np.cos(earth_dec) * np.cos(dec_radians) * np.sin(ra_radians - earth_ra) + dd_draj = np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) * geometric / (ac.c * u.radian) return dd_draj.decompose(u.si.bases) -def d_delay_astrometry_d_DECJ(ssb_to_earth_vector, ra_radians, dec_radians): +def d_delay_d_DECJ(ssb_to_earth_vector, ra_radians, dec_radians): """Calculate the derivative wrt DECJ :param ssb_to_earth_vector: Position of Earth wrt SSB per timestamp @@ -182,22 +193,25 @@ def d_delay_astrometry_d_DECJ(ssb_to_earth_vector, ra_radians, dec_radians): :returns: d_delay / d_DECJ """ - import astropy.units as u - import astropy.constants as const + try: + import astropy.units as u + import astropy.constants as ac + except ImportError: # pragma: no cover + log.error("Astropy required for native astrometry timing models") + raise - earth_ra, earth_dec = np.arctan2(ssb_to_earth_vector[:, 1], ssb_to_earth_vector[:, 0]), np.arcsin( - ssb_to_earth_vector[:, 2] / np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) - ) + earth_ra = np.arctan2(ssb_to_earth_vector[:, 1], ssb_to_earth_vector[:, 0]) + earth_dec = np.arcsin(ssb_to_earth_vector[:, 2] / np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1))) - geom = np.cos(earth_dec) * np.sin(dec_radians) * np.cos(ra_radians - earth_ra) - np.sin(earth_dec) * np.cos( + geometric = np.cos(earth_dec) * np.sin(dec_radians) * np.cos(ra_radians - earth_ra) - np.sin(earth_dec) * np.cos( dec_radians ) - dd_ddecj = np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) * geom / (const.c * u.radian) + dd_ddecj = np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) * geometric / (ac.c * u.radian) return dd_ddecj.decompose(u.si.bases) -def d_delay_astrometry_d_PMRA(mjd_timestamps, ssb_to_earth_vector, ra_radians, posepoch_mjd): +def d_delay_d_PMRA(mjd_timestamps, ssb_to_earth_vector, ra_radians, posepoch_mjd): """Calculate the derivative wrt PMRA :param mjd_timestamps: Timestamps [MJD] @@ -207,23 +221,27 @@ def d_delay_astrometry_d_PMRA(mjd_timestamps, ssb_to_earth_vector, ra_radians, p :returns: d_delay / d_PMRA """ - import astropy.units as u - import astropy.constants as const + try: + import astropy.units as u + import astropy.constants as ac + except ImportError: # pragma: no cover + log.error("Astropy required for native astrometry timing models") + raise earth_ra = np.arctan2(ssb_to_earth_vector[:, 1], ssb_to_earth_vector[:, 0]) - te = (mjd_timestamps - posepoch_mjd) * u.day - geom = np.cos(np.arcsin(ssb_to_earth_vector[:, 2] / np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)))) * np.sin( - ra_radians - earth_ra - ) + time_earth = (mjd_timestamps - posepoch_mjd) * u.day + geometric = np.cos( + np.arcsin(ssb_to_earth_vector[:, 2] / np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1))) + ) * np.sin(ra_radians - earth_ra) - deriv = np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) * geom * te / (const.c * u.radian) - dd_dpmra = deriv * u.mas / u.year + ddelay_pmra = np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) * geometric * time_earth / (ac.c * u.radian) + ddelay_dpmra_u = ddelay_pmra * u.mas / u.year - return dd_dpmra.decompose(u.si.bases) / (u.mas / u.year) + return ddelay_dpmra_u.decompose(u.si.bases) / (u.mas / u.year) -def d_delay_astrometry_d_PMDEC(mjd_timestamps, ssb_to_earth_vector, ra_radians, dec_radians, posepoch_mjd): +def d_delay_d_PMDEC(mjd_timestamps, ssb_to_earth_vector, ra_radians, dec_radians, posepoch_mjd): """Calculate the derivative wrt PMDEC :param mjd_timestamps: Timestamps [MJD] @@ -234,22 +252,25 @@ def d_delay_astrometry_d_PMDEC(mjd_timestamps, ssb_to_earth_vector, ra_radians, :returns: d_delay / d_PMDEC """ - import astropy.units as u - import astropy.constants as const + try: + import astropy.units as u + import astropy.constants as ac + except ImportError: # pragma: no cover + log.error("Astropy required for native astrometry timing models") + raise - earth_ra, earth_dec = np.arctan2(ssb_to_earth_vector[:, 1], ssb_to_earth_vector[:, 0]), np.arcsin( - ssb_to_earth_vector[:, 2] / np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) - ) + earth_ra = np.arctan2(ssb_to_earth_vector[:, 1], ssb_to_earth_vector[:, 0]) + earth_dec = np.arcsin(ssb_to_earth_vector[:, 2] / np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1))) - te = (mjd_timestamps - posepoch_mjd) * u.day - geom = np.cos(earth_dec) * np.sin(dec_radians) * np.cos(ra_radians - earth_ra) - np.cos(dec_radians) * np.sin( + time_earth = (mjd_timestamps - posepoch_mjd) * u.day + geometric = np.cos(earth_dec) * np.sin(dec_radians) * np.cos(ra_radians - earth_ra) - np.cos(dec_radians) * np.sin( earth_dec ) - deriv = np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) * geom * te / (const.c * u.radian) - dd_dpmdec = deriv * u.mas / u.year + ddelay_dpmdec = np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) * geometric * time_earth / (ac.c * u.radian) + ddelay_dpmdec_u = ddelay_dpmdec * u.mas / u.year - return dd_dpmdec.decompose(u.si.bases) / (u.mas / u.year) + return ddelay_dpmdec_u.decompose(u.si.bases) / (u.mas / u.year) def create_astrometry_timing_model(toas, raj, decj, posepoch): @@ -264,21 +285,21 @@ def create_astrometry_timing_model(toas, raj, decj, posepoch): :returns: design matrix, parameter names """ - toas_mjds = toas * 86400.0 - posepoch_mjd = posepoch * 86400.0 + toas_mjds = (toas * u.sec).to(u.day) + posepoch_mjd = (posepoch * u.sec).to(u.day) parameter_names = ["RAJ", "DECJ", "PMRA", "PMDEC", "PX"] designmatrix = np.zeros((len(toas_mjds), 5)) - ssb_to_earth = ssb_to_earth_vector(toas_mjds) + ssb_to_earth = ssb_to_earth_vector(toas_mjds.value) ssb_to_pulsar = ssb_to_pulsar_vector(raj, decj, 1.0) ssb_to_pulsar_norm = ssb_to_pulsar / np.linalg.norm(ssb_to_pulsar) - designmatrix[:, 0] = d_delay_astrometry_d_RAJ(ssb_to_earth, raj, decj).value - designmatrix[:, 1] = d_delay_astrometry_d_DECJ(ssb_to_earth, raj, decj).value - designmatrix[:, 2] = d_delay_astrometry_d_PMRA(toas_mjds, ssb_to_earth, raj, posepoch_mjd).value - designmatrix[:, 3] = d_delay_astrometry_d_PMDEC(toas_mjds, ssb_to_earth, raj, decj, posepoch_mjd).value - designmatrix[:, 4] = d_delay_astrometry_d_PX(ssb_to_earth, ssb_to_pulsar_norm).value + designmatrix[:, 0] = d_delay_d_RAJ(ssb_to_earth, raj, decj).value + designmatrix[:, 1] = d_delay_d_DECJ(ssb_to_earth, raj, decj).value + designmatrix[:, 2] = d_delay_d_PMRA(toas_mjds, ssb_to_earth, raj, posepoch_mjd.value).value + designmatrix[:, 3] = d_delay_d_PMDEC(toas_mjds, ssb_to_earth, raj, decj, posepoch_mjd.value).value + designmatrix[:, 4] = d_delay_d_PX(ssb_to_earth, ssb_to_pulsar_norm).value return designmatrix, parameter_names From 727d3bc21f261bd8db8c4e9309f0caa74d00e3ed Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Tue, 28 Nov 2023 11:54:16 +0100 Subject: [PATCH 07/14] FIXED: units of astrometry derivatives --- enterprise/signals/utils.py | 133 +++++++++++++++++++----------------- 1 file changed, 70 insertions(+), 63 deletions(-) diff --git a/enterprise/signals/utils.py b/enterprise/signals/utils.py index 600857de..55586ed2 100644 --- a/enterprise/signals/utils.py +++ b/enterprise/signals/utils.py @@ -88,6 +88,7 @@ def ssb_to_earth_vector(mjd_timestamps): try: from astropy.time import Time + import astropy.units as u from astropy.coordinates import get_body_barycentric, solar_system_ephemeris, ICRS, GeocentricTrueEcliptic except ImportError: # pragma: no cover log.error("Astropy required for native astrometry timing models") @@ -103,7 +104,7 @@ def ssb_to_earth_vector(mjd_timestamps): earth_positions_icrs = get_body_barycentric("earth", times) # Extract x, y, and z coordinates - return np.array([earth_positions_icrs.x.value, earth_positions_icrs.y.value, earth_positions_icrs.z.value]).T + return np.array([earth_positions_icrs.x.value, earth_positions_icrs.y.value, earth_positions_icrs.z.value]).T * u.AU def ssb_to_pulsar_vector(ra_radians, dec_radians, distance_parsecs): @@ -126,47 +127,22 @@ def ssb_to_pulsar_vector(ra_radians, dec_radians, distance_parsecs): # Create a SkyCoord object with the given RA, DEC, and distance pulsar_coord = SkyCoord( - ra=ra_radians * u.radian, dec=dec_radians * u.radian, distance=distance_parsecs * u.parsec, frame="icrs" + ra=ra_radians, dec=dec_radians, distance=distance_parsecs, frame="icrs" ) # Convert to Cartesian coordinates (x, y, z) return pulsar_coord.cartesian.xyz.value -def d_delay_d_PX(ssb_to_earth_vector, ssb_to_pulsar_vector): - """ - Calculate the derivative wrt PX - - :param ssb_to_earth_vector: Position of Earth wrt SSB per timestamp - :param ssb_to_pulsar_vector: Direction of Pulsar wrt SSB - - :returns: d_delay / d_PX [sec / mas] - """ - try: - import astropy.units as u - import astropy.constants as ac - except ImportError: # pragma: no cover - log.error("Astropy required for native astrometry timing models") - raise - - ssb_earh_r = np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) - in_prod = np.sum(ssb_to_earth_vector * ssb_to_pulsar_vector, axis=1) - - px_r = np.sqrt(ssb_earh_r**2 - in_prod**2) - dd_dpx = 0.5 * (px_r**2 / (u.AU * ac.c)) * (u.mas / u.radian) - - return dd_dpx.decompose(u.si.bases) / u.mas - - -def d_delay_d_RAJ(ssb_to_earth_vector, ra_radians, dec_radians): +def d_delay_d_RAJ(ssb_to_earth_v, ra_radians, dec_radians): """ Calculate the derivative wrt RAJ - :param ssb_to_earth_vector: Position of Earth wrt SSB per timestamp + :param ssb_to_earth_v: Position of Earth wrt SSB per timestamp :param ra_radians: Right ascension of the Pulsar [rad] :param dec_radians: Declination of the Pulsar [rad] - :returns: d_delay / d_RAJ + :returns: d_delay / d_RAJ [sec / rad] """ try: import astropy.units as u @@ -175,23 +151,23 @@ def d_delay_d_RAJ(ssb_to_earth_vector, ra_radians, dec_radians): log.error("Astropy required for native astrometry timing models") raise - earth_ra = np.arctan2(ssb_to_earth_vector[:, 1], ssb_to_earth_vector[:, 0]) - earth_dec = np.arcsin(ssb_to_earth_vector[:, 2] / np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1))) + earth_ra = np.arctan2(ssb_to_earth_v[:, 1], ssb_to_earth_v[:, 0]) + earth_dec = np.arcsin(ssb_to_earth_v[:, 2] / np.sqrt(np.sum(ssb_to_earth_v**2, axis=1))) geometric = np.cos(earth_dec) * np.cos(dec_radians) * np.sin(ra_radians - earth_ra) - dd_draj = np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) * geometric / (ac.c * u.radian) + dd_draj = np.sqrt(np.sum(ssb_to_earth_v**2, axis=1)) * geometric / (ac.c * u.radian) - return dd_draj.decompose(u.si.bases) + return dd_draj.to(u.sec / u.rad) -def d_delay_d_DECJ(ssb_to_earth_vector, ra_radians, dec_radians): +def d_delay_d_DECJ(ssb_to_earth_v, ra_radians, dec_radians): """Calculate the derivative wrt DECJ - :param ssb_to_earth_vector: Position of Earth wrt SSB per timestamp + :param ssb_to_earth_v: Position of Earth wrt SSB per timestamp :param ra_radians: Right ascension of the Pulsar :param dec_radians: Declination of the Pulsar - :returns: d_delay / d_DECJ + :returns: d_delay / d_DECJ [sec / rad] """ try: import astropy.units as u @@ -200,22 +176,22 @@ def d_delay_d_DECJ(ssb_to_earth_vector, ra_radians, dec_radians): log.error("Astropy required for native astrometry timing models") raise - earth_ra = np.arctan2(ssb_to_earth_vector[:, 1], ssb_to_earth_vector[:, 0]) - earth_dec = np.arcsin(ssb_to_earth_vector[:, 2] / np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1))) + earth_ra = np.arctan2(ssb_to_earth_v[:, 1], ssb_to_earth_v[:, 0]) + earth_dec = np.arcsin(ssb_to_earth_v[:, 2] / np.sqrt(np.sum(ssb_to_earth_v**2, axis=1))) geometric = np.cos(earth_dec) * np.sin(dec_radians) * np.cos(ra_radians - earth_ra) - np.sin(earth_dec) * np.cos( dec_radians ) - dd_ddecj = np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) * geometric / (ac.c * u.radian) + dd_ddecj = np.sqrt(np.sum(ssb_to_earth_v**2, axis=1)) * geometric / (ac.c * u.radian) - return dd_ddecj.decompose(u.si.bases) + return dd_ddecj.to(u.sec / u.rad) -def d_delay_d_PMRA(mjd_timestamps, ssb_to_earth_vector, ra_radians, posepoch_mjd): +def d_delay_d_PMRA(mjd_timestamps, ssb_to_earth_v, ra_radians, posepoch_mjd): """Calculate the derivative wrt PMRA :param mjd_timestamps: Timestamps [MJD] - :param ssb_to_earth_vector: Position of Earth wrt SSB per timestamp [rad] + :param ssb_to_earth_v: Position of Earth wrt SSB per timestamp [rad] :param ra_radians: Right ascension of the Pulsar [rad] :param posepoch_mjd: Position epoch [MJD] @@ -228,24 +204,24 @@ def d_delay_d_PMRA(mjd_timestamps, ssb_to_earth_vector, ra_radians, posepoch_mjd log.error("Astropy required for native astrometry timing models") raise - earth_ra = np.arctan2(ssb_to_earth_vector[:, 1], ssb_to_earth_vector[:, 0]) + earth_ra = np.arctan2(ssb_to_earth_v[:, 1], ssb_to_earth_v[:, 0]) - time_earth = (mjd_timestamps - posepoch_mjd) * u.day + time_earth = (mjd_timestamps - posepoch_mjd) geometric = np.cos( - np.arcsin(ssb_to_earth_vector[:, 2] / np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1))) + np.arcsin(ssb_to_earth_v[:, 2] / np.sqrt(np.sum(ssb_to_earth_v**2, axis=1))) ) * np.sin(ra_radians - earth_ra) - ddelay_pmra = np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) * geometric * time_earth / (ac.c * u.radian) - ddelay_dpmra_u = ddelay_pmra * u.mas / u.year + ddelay_pmra = np.sqrt(np.sum(ssb_to_earth_v**2, axis=1)) * geometric * time_earth / (ac.c * u.radian) + #ddelay_dpmra_u = ddelay_pmra * u.mas / u.year - return ddelay_dpmra_u.decompose(u.si.bases) / (u.mas / u.year) + return ddelay_dpmra.decompose(u.si.bases) #/ (u.mas / u.year) -def d_delay_d_PMDEC(mjd_timestamps, ssb_to_earth_vector, ra_radians, dec_radians, posepoch_mjd): +def d_delay_d_PMDEC(mjd_timestamps, ssb_to_earth_v, ra_radians, dec_radians, posepoch_mjd): """Calculate the derivative wrt PMDEC :param mjd_timestamps: Timestamps [MJD] - :param ssb_to_earth_vector: Position of Earth wrt SSB per timestamp + :param ssb_to_earth_v: Position of Earth wrt SSB per timestamp :param ra_radians: Right ascension of the Pulsar [rad] :param dec_radians: Declination of the Pulsar [rad] :param posepoch_mjd: Position epoch [MJD] @@ -259,18 +235,41 @@ def d_delay_d_PMDEC(mjd_timestamps, ssb_to_earth_vector, ra_radians, dec_radians log.error("Astropy required for native astrometry timing models") raise - earth_ra = np.arctan2(ssb_to_earth_vector[:, 1], ssb_to_earth_vector[:, 0]) - earth_dec = np.arcsin(ssb_to_earth_vector[:, 2] / np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1))) + earth_ra = np.arctan2(ssb_to_earth_v[:, 1], ssb_to_earth_v[:, 0]) + earth_dec = np.arcsin(ssb_to_earth_v[:, 2] / np.sqrt(np.sum(ssb_to_earth_v**2, axis=1))) - time_earth = (mjd_timestamps - posepoch_mjd) * u.day + time_earth = (mjd_timestamps - posepoch_mjd) geometric = np.cos(earth_dec) * np.sin(dec_radians) * np.cos(ra_radians - earth_ra) - np.cos(dec_radians) * np.sin( earth_dec ) - ddelay_dpmdec = np.sqrt(np.sum(ssb_to_earth_vector**2, axis=1)) * geometric * time_earth / (ac.c * u.radian) - ddelay_dpmdec_u = ddelay_dpmdec * u.mas / u.year + ddelay_dpmdec = np.sqrt(np.sum(ssb_to_earth_v**2, axis=1)) * geometric * time_earth / (ac.c * u.radian) + return ddelay_dpmdec.decompose(u.si.bases) + + +def d_delay_d_PX(ssb_to_earth_v, ssb_to_pulsar_v): + """ + Calculate the derivative wrt PX + + :param ssb_to_earth_v: Position of Earth wrt SSB per timestamp + :param ssb_to_pulsar_v: Direction of Pulsar wrt SSB + + :returns: d_delay / d_PX [sec / mas] + """ + try: + import astropy.units as u + import astropy.constants as ac + except ImportError: # pragma: no cover + log.error("Astropy required for native astrometry timing models") + raise + + ssb_earh_r = np.sqrt(np.sum(ssb_to_earth_v**2, axis=1)) + in_prod = np.sum(ssb_to_earth_v * ssb_to_pulsar_v, axis=1) - return ddelay_dpmdec_u.decompose(u.si.bases) / (u.mas / u.year) + px_radius = np.sqrt(ssb_earh_r**2 - in_prod**2) + dd_dpx = 0.5 * (px_radius**2 / (u.AU * ac.c)) / u.radian + + return dd_dpx.to(u.second / u.mas) def create_astrometry_timing_model(toas, raj, decj, posepoch): @@ -284,21 +283,29 @@ def create_astrometry_timing_model(toas, raj, decj, posepoch): :returns: design matrix, parameter names """ + try: + import astropy.units as u + except ImportError: # pragma: no cover + log.error("Astropy required for native astrometry timing models") + raise + + raj = raj * u.rad + decj = decj * u.rad - toas_mjds = (toas * u.sec).to(u.day) - posepoch_mjd = (posepoch * u.sec).to(u.day) + toas_mjds = (toas * u.second).to(u.day) + posepoch_mjd = (posepoch * u.second).to(u.day) parameter_names = ["RAJ", "DECJ", "PMRA", "PMDEC", "PX"] designmatrix = np.zeros((len(toas_mjds), 5)) - ssb_to_earth = ssb_to_earth_vector(toas_mjds.value) - ssb_to_pulsar = ssb_to_pulsar_vector(raj, decj, 1.0) + ssb_to_earth = ssb_to_earth_vector(toas_mjds) + ssb_to_pulsar = ssb_to_pulsar_vector(raj, decj, 1.0 * u.parsec) ssb_to_pulsar_norm = ssb_to_pulsar / np.linalg.norm(ssb_to_pulsar) designmatrix[:, 0] = d_delay_d_RAJ(ssb_to_earth, raj, decj).value designmatrix[:, 1] = d_delay_d_DECJ(ssb_to_earth, raj, decj).value - designmatrix[:, 2] = d_delay_d_PMRA(toas_mjds, ssb_to_earth, raj, posepoch_mjd.value).value - designmatrix[:, 3] = d_delay_d_PMDEC(toas_mjds, ssb_to_earth, raj, decj, posepoch_mjd.value).value + designmatrix[:, 2] = d_delay_d_PMRA(toas_mjds, ssb_to_earth, raj, posepoch_mjd).value + designmatrix[:, 3] = d_delay_d_PMDEC(toas_mjds, ssb_to_earth, raj, decj, posepoch_mjd).value designmatrix[:, 4] = d_delay_d_PX(ssb_to_earth, ssb_to_pulsar_norm).value return designmatrix, parameter_names From 2396b0dfb99c6a2cd3c9d2c29932adacb44ff52d Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Tue, 28 Nov 2023 13:22:21 +0100 Subject: [PATCH 08/14] Added Astrometry derivatives unit tests --- enterprise/pulsar.py | 6 ++--- enterprise/signals/utils.py | 25 ++++++++++----------- tests/test_utils.py | 44 +++++++++++++++++++++++++++++++++++++ 3 files changed, 59 insertions(+), 16 deletions(-) diff --git a/enterprise/pulsar.py b/enterprise/pulsar.py index bfc1c1b2..8fe28e31 100644 --- a/enterprise/pulsar.py +++ b/enterprise/pulsar.py @@ -10,7 +10,7 @@ import astropy.constants as const import astropy.units as u import numpy as np -from ephem import Ecliptic, Equatorial +from ephem import Ecliptic, Equatorial, J2000 import enterprise from enterprise.signals import utils @@ -626,8 +626,8 @@ def __init__( self.name = utils.get_psrname_from_pos(elong=elong, elat=elat, raj=raj, decj=decj) if elong is not None and elat is not None: - ec = ephem.Ecliptic(elong * np.pi / 180.0, elat * np.pi / 180.0) - eq = ephem.Equatorial(ec, epoch=ephem.J2000) + ec = Ecliptic(elong * np.pi / 180.0, elat * np.pi / 180.0) + eq = Equatorial(ec, epoch=J2000) raj, decj = np.double(eq.ra), np.double(eq.dec) self._raj, self._decj = raj, decj diff --git a/enterprise/signals/utils.py b/enterprise/signals/utils.py index 55586ed2..4b4b66e6 100644 --- a/enterprise/signals/utils.py +++ b/enterprise/signals/utils.py @@ -126,9 +126,7 @@ def ssb_to_pulsar_vector(ra_radians, dec_radians, distance_parsecs): raise # Create a SkyCoord object with the given RA, DEC, and distance - pulsar_coord = SkyCoord( - ra=ra_radians, dec=dec_radians, distance=distance_parsecs, frame="icrs" - ) + pulsar_coord = SkyCoord(ra=ra_radians, dec=dec_radians, distance=distance_parsecs, frame="icrs") # Convert to Cartesian coordinates (x, y, z) return pulsar_coord.cartesian.xyz.value @@ -157,7 +155,7 @@ def d_delay_d_RAJ(ssb_to_earth_v, ra_radians, dec_radians): geometric = np.cos(earth_dec) * np.cos(dec_radians) * np.sin(ra_radians - earth_ra) dd_draj = np.sqrt(np.sum(ssb_to_earth_v**2, axis=1)) * geometric / (ac.c * u.radian) - return dd_draj.to(u.sec / u.rad) + return dd_draj.to(u.second / u.rad) def d_delay_d_DECJ(ssb_to_earth_v, ra_radians, dec_radians): @@ -184,7 +182,7 @@ def d_delay_d_DECJ(ssb_to_earth_v, ra_radians, dec_radians): ) dd_ddecj = np.sqrt(np.sum(ssb_to_earth_v**2, axis=1)) * geometric / (ac.c * u.radian) - return dd_ddecj.to(u.sec / u.rad) + return dd_ddecj.to(u.second / u.rad) def d_delay_d_PMRA(mjd_timestamps, ssb_to_earth_v, ra_radians, posepoch_mjd): @@ -206,15 +204,15 @@ def d_delay_d_PMRA(mjd_timestamps, ssb_to_earth_v, ra_radians, posepoch_mjd): earth_ra = np.arctan2(ssb_to_earth_v[:, 1], ssb_to_earth_v[:, 0]) - time_earth = (mjd_timestamps - posepoch_mjd) - geometric = np.cos( - np.arcsin(ssb_to_earth_v[:, 2] / np.sqrt(np.sum(ssb_to_earth_v**2, axis=1))) - ) * np.sin(ra_radians - earth_ra) + time_earth = mjd_timestamps - posepoch_mjd + geometric = np.cos(np.arcsin(ssb_to_earth_v[:, 2] / np.sqrt(np.sum(ssb_to_earth_v**2, axis=1)))) * np.sin( + ra_radians - earth_ra + ) - ddelay_pmra = np.sqrt(np.sum(ssb_to_earth_v**2, axis=1)) * geometric * time_earth / (ac.c * u.radian) - #ddelay_dpmra_u = ddelay_pmra * u.mas / u.year + ddelay_dpmra = np.sqrt(np.sum(ssb_to_earth_v**2, axis=1)) * geometric * time_earth / (ac.c * u.radian) - return ddelay_dpmra.decompose(u.si.bases) #/ (u.mas / u.year) + # TODO: these units are not correct + return ddelay_dpmra.decompose(u.si.bases) def d_delay_d_PMDEC(mjd_timestamps, ssb_to_earth_v, ra_radians, dec_radians, posepoch_mjd): @@ -238,11 +236,12 @@ def d_delay_d_PMDEC(mjd_timestamps, ssb_to_earth_v, ra_radians, dec_radians, pos earth_ra = np.arctan2(ssb_to_earth_v[:, 1], ssb_to_earth_v[:, 0]) earth_dec = np.arcsin(ssb_to_earth_v[:, 2] / np.sqrt(np.sum(ssb_to_earth_v**2, axis=1))) - time_earth = (mjd_timestamps - posepoch_mjd) + time_earth = mjd_timestamps - posepoch_mjd geometric = np.cos(earth_dec) * np.sin(dec_radians) * np.cos(ra_radians - earth_ra) - np.cos(dec_radians) * np.sin( earth_dec ) + # TODO: these units are not correct ddelay_dpmdec = np.sqrt(np.sum(ssb_to_earth_v**2, axis=1)) * geometric * time_earth / (ac.c * u.radian) return ddelay_dpmdec.decompose(u.si.bases) diff --git a/tests/test_utils.py b/tests/test_utils.py index a07af86e..22bf46d2 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -217,3 +217,47 @@ def test_orf(self): vals = [(hd, hd_exp), (dp, dp_exp), (mp, mp_exp), (anis_orf, anis_orf_exp)] for key, val in zip(keys, vals): assert val[0] == val[1], msg.format(key) + + +class TestAstrometry(unittest.TestCase): + @classmethod + def setUpClass(cls): + """Setup the Pulsar object.""" + + # initialize Pulsar class that uses Equatorial Coordinates + cls.psr = Pulsar( + datadir + "/1713.Sep.T2.par", datadir + "/1713.Sep.T2.tim", timing_package="tempo2", drop_t2pulsar=False + ) + + cls.Mmat = cls.psr.t2pulsar.designmatrix(fixunits=False, fixsigns=True, incoffset=True) + cls.posepoch = cls.psr.t2pulsar["POSEPOCH"].val * 86400.0 + cls.dm, cls.dmp = utils.create_astrometry_timing_model(cls.psr.toas, cls.psr._raj, cls.psr._decj, cls.posepoch) + + def test_ddelay_dastrometry(self): + """Test the derivatives of the astrometry parameters""" + + incorrect_units = ["PMRA", "PMDEC"] + + for pp, pname in enumerate(self.dmp): + dmc_self = self.dm[:, pp] + + t2dm = ("Offset",) + self.psr.t2pulsar.pars(which="fit") + pind = t2dm.index(pname) + dmc_t2 = self.Mmat[self.psr._isort, pind] + a = dmc_t2 + b = dmc_self + + if pname not in incorrect_units: + rel_diff = (np.abs(a) - np.abs(b)) / (np.abs(a) + np.abs(b)) + + msg = f"ddelay_d{pname} is not consistent with Tempo2" + assert np.allclose(rel_diff, 0.0, atol=0.05), msg + + else: + # Proper motion doesn't have the right units + conv = np.mean(np.abs(a) / np.abs(b)) + + rel_diff = (np.abs(a) - conv * np.abs(b)) / (np.abs(a) + conv * np.abs(b)) + + msg = f"ddelay_d{pname} is not consistent with Tempo2" + assert np.allclose(rel_diff, 0.0, atol=0.05), msg From 50825c1dd493c5c394b364c7f6c10f2694473475 Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Tue, 28 Nov 2023 13:50:20 +0100 Subject: [PATCH 09/14] Fixed flake8 suggestions --- enterprise/signals/utils.py | 25 ++++++++++++------------- 1 file changed, 12 insertions(+), 13 deletions(-) diff --git a/enterprise/signals/utils.py b/enterprise/signals/utils.py index 4b4b66e6..15e12682 100644 --- a/enterprise/signals/utils.py +++ b/enterprise/signals/utils.py @@ -14,7 +14,7 @@ from scipy.integrate import odeint from scipy.interpolate import interp1d from sksparse.cholmod import cholesky -from ephem import Ecliptic, Equatorial +from ephem import Ecliptic, Equatorial, J2000 import enterprise from enterprise import constants as const @@ -68,8 +68,8 @@ def get_psrname_from_pos(elong=None, elat=None, raj=None, decj=None): """ if elong is not None and elat is not None: - ec = ephem.Ecliptic(elong * np.pi / 180.0, elat * np.pi / 180.0) - eq = ephem.Equatorial(ec, epoch=ephem.J2000) + ec = Ecliptic(elong * np.pi / 180.0, elat * np.pi / 180.0) + eq = Equatorial(ec, epoch=J2000) raj, decj = float(eq.ra), float(eq.dec) elif raj is None or decj is None: raise ValueError("Need to provide either raj/decj or elong/elat!") @@ -89,9 +89,9 @@ def ssb_to_earth_vector(mjd_timestamps): try: from astropy.time import Time import astropy.units as u - from astropy.coordinates import get_body_barycentric, solar_system_ephemeris, ICRS, GeocentricTrueEcliptic + from astropy.coordinates import get_body_barycentric, solar_system_ephemeris except ImportError: # pragma: no cover - log.error("Astropy required for native astrometry timing models") + logger.error("Astropy required for native astrometry timing models") raise # Set solar system ephemeris to 'builtin' for offline calculations @@ -120,9 +120,8 @@ def ssb_to_pulsar_vector(ra_radians, dec_radians, distance_parsecs): try: from astropy.coordinates import SkyCoord - import astropy.units as u except ImportError: # pragma: no cover - log.error("Astropy required for native astrometry timing models") + logger.error("Astropy required for native astrometry timing models") raise # Create a SkyCoord object with the given RA, DEC, and distance @@ -146,7 +145,7 @@ def d_delay_d_RAJ(ssb_to_earth_v, ra_radians, dec_radians): import astropy.units as u import astropy.constants as ac except ImportError: # pragma: no cover - log.error("Astropy required for native astrometry timing models") + logger.error("Astropy required for native astrometry timing models") raise earth_ra = np.arctan2(ssb_to_earth_v[:, 1], ssb_to_earth_v[:, 0]) @@ -171,7 +170,7 @@ def d_delay_d_DECJ(ssb_to_earth_v, ra_radians, dec_radians): import astropy.units as u import astropy.constants as ac except ImportError: # pragma: no cover - log.error("Astropy required for native astrometry timing models") + logger.error("Astropy required for native astrometry timing models") raise earth_ra = np.arctan2(ssb_to_earth_v[:, 1], ssb_to_earth_v[:, 0]) @@ -199,7 +198,7 @@ def d_delay_d_PMRA(mjd_timestamps, ssb_to_earth_v, ra_radians, posepoch_mjd): import astropy.units as u import astropy.constants as ac except ImportError: # pragma: no cover - log.error("Astropy required for native astrometry timing models") + logger.error("Astropy required for native astrometry timing models") raise earth_ra = np.arctan2(ssb_to_earth_v[:, 1], ssb_to_earth_v[:, 0]) @@ -230,7 +229,7 @@ def d_delay_d_PMDEC(mjd_timestamps, ssb_to_earth_v, ra_radians, dec_radians, pos import astropy.units as u import astropy.constants as ac except ImportError: # pragma: no cover - log.error("Astropy required for native astrometry timing models") + logger.error("Astropy required for native astrometry timing models") raise earth_ra = np.arctan2(ssb_to_earth_v[:, 1], ssb_to_earth_v[:, 0]) @@ -259,7 +258,7 @@ def d_delay_d_PX(ssb_to_earth_v, ssb_to_pulsar_v): import astropy.units as u import astropy.constants as ac except ImportError: # pragma: no cover - log.error("Astropy required for native astrometry timing models") + logger.error("Astropy required for native astrometry timing models") raise ssb_earh_r = np.sqrt(np.sum(ssb_to_earth_v**2, axis=1)) @@ -285,7 +284,7 @@ def create_astrometry_timing_model(toas, raj, decj, posepoch): try: import astropy.units as u except ImportError: # pragma: no cover - log.error("Astropy required for native astrometry timing models") + logger.error("Astropy required for native astrometry timing models") raise raj = raj * u.rad From af0264be10de52609a78ea1b0163b6af1f6bfe17 Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Tue, 28 Nov 2023 14:53:37 +0100 Subject: [PATCH 10/14] Added MockPulsar unittests --- enterprise/pulsar.py | 5 ++-- enterprise/signals/utils.py | 2 +- tests/test_pulsar.py | 59 ++++++++++++++++++++++++++++++++++++- 3 files changed, 61 insertions(+), 5 deletions(-) diff --git a/enterprise/pulsar.py b/enterprise/pulsar.py index 8fe28e31..f93710c7 100644 --- a/enterprise/pulsar.py +++ b/enterprise/pulsar.py @@ -619,6 +619,7 @@ def __init__( residuals=None, toaerrs=1e-6, sort=True, + flags={}, telescope="GBT", spindown_order=2, inc_astrometry=True, @@ -637,7 +638,7 @@ def __init__( self._toas = np.double(obs_times_mjd) * 86400 self._stoas = np.double(obs_times_mjd) * 86400 - self._residuals = residuals if residuals else np.zeros_like(obs_times_mjd) + self._residuals = residuals if residuals is not None else np.zeros_like(obs_times_mjd) self._toaerrs = np.ones_like(obs_times_mjd) * toaerrs self._posepoch = np.mean(self._toas) @@ -653,8 +654,6 @@ def __init__( # set parameters self.setpars = [fp for fp in self.fitpars] - flags = {} - # new-style storage of flags as a numpy record array (previously, psr._flags = flags) self._flags = np.zeros(len(self._toas), dtype=[(key, val.dtype) for key, val in flags.items()]) for key, val in flags.items(): diff --git a/enterprise/signals/utils.py b/enterprise/signals/utils.py index 15e12682..d1aa38f4 100644 --- a/enterprise/signals/utils.py +++ b/enterprise/signals/utils.py @@ -321,7 +321,7 @@ def create_spindown_timing_model(toas, order=2): avetoas = (toas - np.mean(toas)) / np.mean(toas) designmatrix = np.vstack([avetoas**ii for ii in range(1 + order)]).T - parameter_names = ["Offset"] + ["F{ii}".format(ii) for ii in range(order)] + parameter_names = ["Offset"] + ["F{ii}".format(ii=ii) for ii in range(order)] return designmatrix, parameter_names diff --git a/tests/test_pulsar.py b/tests/test_pulsar.py index 6ab617e3..4a8d1dda 100644 --- a/tests/test_pulsar.py +++ b/tests/test_pulsar.py @@ -18,7 +18,7 @@ import numpy as np -from enterprise.pulsar import Pulsar +from enterprise.pulsar import Pulsar, MockPulsar from tests.enterprise_test_data import datadir import pint.models.timing_model @@ -230,3 +230,60 @@ def test_no_planet(self): msg += "`planet` flag is not True in `toas` or further Pint " msg += "development to add additional planets is needed." self.assertTrue(msg in context.exception) + + +class TestPulsarMock(TestPulsar): + @classmethod + def setUpClass(cls): + """Setup the Pulsar object.""" + + toas = np.linspace(53000.0, 58000.0, 4005) + flags = {"f": np.array(["nosystem"] * 4005), "fe": np.array(["nofrontend"] * 4005)} + decj, raj = (0.16848694562363042, 4.9533700839400492) + + cls.psr = MockPulsar( + obs_times_mjd=toas, + raj=raj, + decj=decj, + ssbfreqs=1440.0 * np.ones_like(toas), + residuals=np.zeros_like(toas), + toaerrs=1e-6 * np.ones_like(toas), + sort=True, + flags=flags, + telescope="GBT", + spindown_order=2, + inc_astrometry=True, + ) + + def test_deflate_inflate(self): + pass + + def test_tearDownClass(self): + pass + + def test_dm(self): + pass + + def test_design_matrix(self): + """Check design matrix shape.""" + + msg = "Design matrix shape incorrect." + assert self.psr.Mmat.shape == (4005, 8), msg + + def test_planetssb(self): + """Place holder for filter_data tests.""" + assert hasattr(self.psr, "_planetssb") + + def test_sunssb(self): + """Place holder for filter_data tests.""" + assert hasattr(self.psr, "_sunssb") + + def test_wrong_input(self): + pass + + def test_value_error(self): + pass + + def test_to_pickle(self): + """Place holder for to_pickle tests.""" + pass From 94bb7df0b175fad12aa636aaa63eacb0f3e79e9b Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Tue, 28 Nov 2023 22:30:02 +0100 Subject: [PATCH 11/14] Added unit tests for MockPulsar and related additions --- tests/test_pulsar.py | 31 +++++++++++++++++++++++++++++-- tests/test_utils.py | 23 +++++++++++++++++++++++ 2 files changed, 52 insertions(+), 2 deletions(-) diff --git a/tests/test_pulsar.py b/tests/test_pulsar.py index 4a8d1dda..47b91bc5 100644 --- a/tests/test_pulsar.py +++ b/tests/test_pulsar.py @@ -23,6 +23,7 @@ import pint.models.timing_model from pint.models import get_model_and_toas +import ephem class TestPulsar(unittest.TestCase): @@ -240,11 +241,14 @@ def setUpClass(cls): toas = np.linspace(53000.0, 58000.0, 4005) flags = {"f": np.array(["nosystem"] * 4005), "fe": np.array(["nofrontend"] * 4005)} decj, raj = (0.16848694562363042, 4.9533700839400492) + eq = ephem.Equatorial(raj, decj, epoch=ephem.J2000) + ec = ephem.Ecliptic(eq) + elong, elat = ec.lon * 180 / np.pi, ec.lat * 180 / np.pi cls.psr = MockPulsar( obs_times_mjd=toas, - raj=raj, - decj=decj, + elong=elong, + elat=elat, ssbfreqs=1440.0 * np.ones_like(toas), residuals=np.zeros_like(toas), toaerrs=1e-6 * np.ones_like(toas), @@ -255,6 +259,21 @@ def setUpClass(cls): inc_astrometry=True, ) + cls.psr_spin = MockPulsar( + obs_times_mjd=toas, + elong=elong, + elat=elat, + ssbfreqs=1440.0 * np.ones_like(toas), + residuals=np.zeros_like(toas), + toaerrs=1e-6 * np.ones_like(toas), + sort=True, + flags=flags, + telescope="GBT", + spindown_order=2, + inc_astrometry=False, + ) + cls.psr_spin.set_residuals(np.ones_like(toas)) + def test_deflate_inflate(self): pass @@ -270,6 +289,8 @@ def test_design_matrix(self): msg = "Design matrix shape incorrect." assert self.psr.Mmat.shape == (4005, 8), msg + assert self.psr_spin.Mmat.shape == (4005, 3), msg + def test_planetssb(self): """Place holder for filter_data tests.""" assert hasattr(self.psr, "_planetssb") @@ -287,3 +308,9 @@ def test_value_error(self): def test_to_pickle(self): """Place holder for to_pickle tests.""" pass + + def test_set_residuals(self): + """Check Residual shape.""" + + msg = "Residuals from set_residuals incorrect" + assert np.all(self.psr_spin.residuals == np.ones(4005)), msg diff --git a/tests/test_utils.py b/tests/test_utils.py index 22bf46d2..307d4e56 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -18,6 +18,8 @@ from enterprise.signals import utils from tests.enterprise_test_data import datadir +import ephem + class TestUtils(unittest.TestCase): @classmethod @@ -224,6 +226,7 @@ class TestAstrometry(unittest.TestCase): def setUpClass(cls): """Setup the Pulsar object.""" + # TODO: use a different pulsar so it's faster? # initialize Pulsar class that uses Equatorial Coordinates cls.psr = Pulsar( datadir + "/1713.Sep.T2.par", datadir + "/1713.Sep.T2.tim", timing_package="tempo2", drop_t2pulsar=False @@ -261,3 +264,23 @@ def test_ddelay_dastrometry(self): msg = f"ddelay_d{pname} is not consistent with Tempo2" assert np.allclose(rel_diff, 0.0, atol=0.05), msg + + def test_get_psrname_from_pos(self): + """Test the functionality to derive pulsar names""" + + # Pulsar B1855+09 (= J1857+09..) + decj, raj = (0.16848694562363042, 4.9533700839400492) + eq = ephem.Equatorial(raj, decj, epoch=ephem.J2000) + ec = ephem.Ecliptic(eq) + elong, elat = ec.lon * 180 / np.pi, ec.lat * 180 / np.pi + + msg = "Name from elong/elat not consistent with real pulsar name" + psrname = utils.get_psrname_from_pos(elong=elong, elat=elat, raj=None, decj=None) + assert psrname == "J1855+0939", msg + + msg = "Name from raj/decj not consistent with real pulsar name" + psrname = utils.get_psrname_from_pos(elong=None, elat=None, raj=raj, decj=decj) + assert psrname == "J1855+0939", msg + + with self.assertRaises(ValueError) as context: + psrname = utils.get_psrname_from_pos(elong=None, elat=None, raj=None, decj=None) From c090d036a534d46441a220c3b50f36c7923e3bcb Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Wed, 29 Nov 2023 06:53:08 +0100 Subject: [PATCH 12/14] Removed context variable for linting --- tests/test_utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_utils.py b/tests/test_utils.py index 307d4e56..d38584c7 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -282,5 +282,5 @@ def test_get_psrname_from_pos(self): psrname = utils.get_psrname_from_pos(elong=None, elat=None, raj=raj, decj=decj) assert psrname == "J1855+0939", msg - with self.assertRaises(ValueError) as context: + with self.assertRaises(ValueError): psrname = utils.get_psrname_from_pos(elong=None, elat=None, raj=None, decj=None) From dbf5779eb53c5880228f4c6995aaf316cd6fae63 Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Wed, 29 Nov 2023 09:52:19 +0100 Subject: [PATCH 13/14] Turn off astrometry parameters if there's no astropy, and issue a warning --- enterprise/pulsar.py | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/enterprise/pulsar.py b/enterprise/pulsar.py index f93710c7..57350b1c 100644 --- a/enterprise/pulsar.py +++ b/enterprise/pulsar.py @@ -603,10 +603,7 @@ def destroy(psr): # pragma: py-lt-38 class MockPulsar(BasePulsar): """Class to allow mock pulsars to be used with Enterprise""" - # TODO: allow units? - # test utils.get_psrname_from_raj_decj - # utils.get_psrname_from_pos - # utils.create_spindown_timing_model + _noastropy_warning_issued = False # Class variable def __init__( self, @@ -624,6 +621,16 @@ def __init__( spindown_order=2, inc_astrometry=True, ): + if inc_astrometry and not hasattr(const, "c"): # pragma: no cover + # We requested astromery parameters, but there's no astropy + if not MockPulsar._noastropy_warning_issued: + logger.warning( + "Astropy not installed but requested astrometry timing model. Switching off astrometry parameters in all instances of MockPulsar." + ) + MockPulsar._noastropy_warning_issued = True + + inc_astrometry = False + self.name = utils.get_psrname_from_pos(elong=elong, elat=elat, raj=raj, decj=decj) if elong is not None and elat is not None: From 143bb6b864ab25c0d086f29d442a2a93900f1afb Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Wed, 29 Nov 2023 10:00:59 +0100 Subject: [PATCH 14/14] Linting: too long a line in pulsar.py --- enterprise/pulsar.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/enterprise/pulsar.py b/enterprise/pulsar.py index 57350b1c..71bb437a 100644 --- a/enterprise/pulsar.py +++ b/enterprise/pulsar.py @@ -624,9 +624,11 @@ def __init__( if inc_astrometry and not hasattr(const, "c"): # pragma: no cover # We requested astromery parameters, but there's no astropy if not MockPulsar._noastropy_warning_issued: - logger.warning( - "Astropy not installed but requested astrometry timing model. Switching off astrometry parameters in all instances of MockPulsar." - ) + msg = "WARNING: Astropy not installed but user requested " + msg += "astrometry timing model in MockPulsar. " + msg += "Switching off astrometry in all instances of MockPulsar." + logger.warning(msg) + MockPulsar._noastropy_warning_issued = True inc_astrometry = False