From 37457340e43ba71878478b9f8fdb2409124e33e1 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Fri, 29 Sep 2023 09:48:47 -0500 Subject: [PATCH 01/21] faster implementations of astrometric ssb calculations --- CHANGELOG-unreleased.md | 2 + src/pint/models/astrometry.py | 181 ++++++++++++++++++++++++++++++++-- tests/test_astrometry.py | 80 +++++++++++++++ 3 files changed, 257 insertions(+), 6 deletions(-) diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index bd2cb7039..21722454a 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -12,6 +12,8 @@ the released changes. - `WAVE` parameters can be added to a `Wave` model with `add_wave_component()` in `wave.py` - Moved design matrix normalization code from `pint.fitter` to the new `pint.utils.normalize_designmatrix()` function. - Made `Residuals` independent of `GLSFitter` (GLS chi2 is now computed using the new function `Residuals._calc_gls_chi2()`). +- `ssb_to_psb_ICRS` implementation is now a lot faster (uses erfa functions directly) +- `ssb_to_psb_ECL` implementation within `AstrometryEcliptic` model is now a lot faster (uses erfa functions directly) ### Added - CHI2, CHI2R, TRES, DMRES now in postfit par files - Added `WaveX` model as a `DelayComponent` with Fourier amplitudes as fitted parameters diff --git a/src/pint/models/astrometry.py b/src/pint/models/astrometry.py index 1d737849c..8867ffe6b 100644 --- a/src/pint/models/astrometry.py +++ b/src/pint/models/astrometry.py @@ -11,7 +11,7 @@ from loguru import logger as log -from erfa import ErfaWarning +from erfa import ErfaWarning, pmsafe from pint import ls from pint.models.parameter import ( AngleParameter, @@ -61,9 +61,59 @@ def ssb_to_psb_xyz_ICRS(self, epoch=None): """Returns unit vector(s) from SSB to pulsar system barycenter under ICRS. If epochs (MJD) are given, proper motion is included in the calculation. + + Parameters + ---------- + epoch : float or astropy.time.Time, optional + + Returns + ------- + np.ndarray : + (len(epoch), 3) array of unit vectors """ # TODO: would it be better for this to return a 6-vector (pos, vel)? - return self.coords_as_ICRS(epoch=epoch).cartesian.xyz.transpose() + + # this is somewhat slow, since it repeatedly created difference SkyCoord Objects + # return self.coords_as_ICRS(epoch=epoch).cartesian.xyz.transpose() + + # Instead look at what https://docs.astropy.org/en/stable/_modules/astropy/coordinates/sky_coordinate.html#SkyCoord.apply_space_motion + # does, which is to use https://github.com/liberfa/erfa/blob/master/src/starpm.c + # and then just use the relevant pieces of that + if epoch is None: + return self.coords_as_ICRS(epoch=epoch).cartesian.xyz.transpose() + epoch = ( + epoch + if isinstance(epoch, Time) + else Time(epoch, scale="tdb", format="pulsar_mjd") + ) + # in the general case we don't know what system the coordinates will be in + # so explicitly transform to ICRS + cc = self.get_psr_coords() + icrsrep = cc.icrs.represent_as( + coords.SphericalRepresentation, + coords.SphericalDifferential, + ) + icrsvel = icrsrep.differentials["s"] + with warnings.catch_warnings(): + warnings.simplefilter("ignore", ErfaWarning) + starpm = pmsafe( + icrsrep.lon.radian, + icrsrep.lat.radian, + icrsvel.d_lon.to_value(u.radian / u.yr), + icrsvel.d_lat.to_value(u.radian / u.yr), + 0.0, + 0.0, + self.POSEPOCH.quantity.jd1, + self.POSEPOCH.quantity.jd2, + epoch.jd1, + epoch.jd2, + ) + # ra,dec now in radians + ra, dec = starpm[0], starpm[1] + x = np.cos(ra) * np.cos(dec) + y = np.sin(ra) * np.cos(dec) + z = np.sin(dec) + return u.Quantity([x, y, z]).T def ssb_to_psb_xyz_ECL(self, epoch=None, ecl=None): """Returns unit vector(s) from SSB to pulsar system barycenter under Ecliptic coordinates. @@ -75,6 +125,11 @@ def ssb_to_psb_xyz_ECL(self, epoch=None, ecl=None): epoch : float, optional ecl : str, optional Obliquity (IERS2010 by default) + + Returns + ------- + np.ndarray : + (len(epoch), 3) array of unit vectors """ # TODO: would it be better for this to return a 6-vector (pos, vel)? return self.coords_as_ECL(epoch=epoch, ecl=ecl).cartesian.xyz.transpose() @@ -346,7 +401,9 @@ def get_psr_coords(self, epoch=None): frame=coords.ICRS, ) newepoch = ( - epoch if isinstance(epoch, Time) else Time(epoch, scale="tdb", format="mjd") + epoch + if isinstance(epoch, Time) + else Time(epoch, scale="tdb", format="pulsar_mjd") ) position_now = add_dummy_distance(self.get_psr_coords()) with warnings.catch_warnings(): @@ -390,6 +447,60 @@ def get_params_as_ICRS(self): "PMDEC": self.PMDEC.quantity, } + def ssb_to_psb_xyz_ICRS(self, epoch=None): + """Returns unit vector(s) from SSB to pulsar system barycenter under ICRS. + + If epochs (MJD) are given, proper motion is included in the calculation. + + Parameters + ---------- + epoch : float or astropy.time.Time, optional + + Returns + ------- + np.ndarray : + (len(epoch), 3) array of unit vectors + """ + # TODO: would it be better for this to return a 6-vector (pos, vel)? + + # this was somewhat slow, since it repeatedly created difference SkyCoord Objects + # return self.coords_as_ICRS(epoch=epoch).cartesian.xyz.transpose() + + # Instead look at what https://docs.astropy.org/en/stable/_modules/astropy/coordinates/sky_coordinate.html#SkyCoord.apply_space_motion + # does, which is to use https://github.com/liberfa/erfa/blob/master/src/starpm.c + # and then just use the relevant pieces of that + if epoch is None: + return self.coords_as_ICRS(epoch=epoch).cartesian.xyz.transpose() + epoch = ( + epoch + if isinstance(epoch, Time) + else Time(epoch, scale="tdb", format="pulsar_mjd") + ) + # compared to the general case above we can assume that the coordinates are ICRS + # so just access those components + with warnings.catch_warnings(): + warnings.simplefilter("ignore", ErfaWarning) + # note that pmsafe wants mu_alpha not mu_alpha * cos(delta) + starpm = pmsafe( + self.RAJ.quantity.to_value(u.radian), + self.DECJ.quantity.to_value(u.radian), + self.PMRA.quantity.to_value(u.radian / u.yr) + / np.cos(self.DECJ.quantity).value, + self.PMDEC.quantity.to_value(u.radian / u.yr), + 0.0, + 0.0, + self.POSEPOCH.quantity.jd1, + self.POSEPOCH.quantity.jd2, + epoch.jd1, + epoch.jd2, + ) + # ra,dec now in radians + ra, dec = starpm[0], starpm[1] + x = np.cos(ra) * np.cos(dec) + y = np.sin(ra) * np.cos(dec) + z = np.sin(dec) + return u.Quantity([x, y, z]).T + def d_delay_astrometry_d_RAJ(self, toas, param="", acc_delay=None): """Calculate the derivative wrt RAJ @@ -486,7 +597,7 @@ def change_posepoch(self, new_epoch): if isinstance(new_epoch, Time): new_epoch = Time(new_epoch, scale="tdb", precision=9) else: - new_epoch = Time(new_epoch, scale="tdb", format="mjd", precision=9) + new_epoch = Time(new_epoch, scale="tdb", format="pulsar_mjd", precision=9) if self.POSEPOCH.value is None: raise ValueError("POSEPOCH is not currently set.") @@ -725,7 +836,9 @@ def get_psr_coords(self, epoch=None): ) # Compute for each time because there is proper motion newepoch = ( - epoch if isinstance(epoch, Time) else Time(epoch, scale="tdb", format="mjd") + epoch + if isinstance(epoch, Time) + else Time(epoch, scale="tdb", format="pulsar_mjd") ) position_now = add_dummy_distance(self.get_psr_coords()) return remove_dummy_distance( @@ -754,6 +867,62 @@ def coords_as_ECL(self, epoch=None, ecl=None): pos_ecl = pos_ecl.transform_to(PulsarEcliptic(ecl=ecl)) return pos_ecl + def ssb_to_psb_xyz_ECL(self, epoch=None, ecl=None): + """Returns unit vector(s) from SSB to pulsar system barycenter under ECL. + + If epochs (MJD) are given, proper motion is included in the calculation. + + Parameters + ---------- + epoch : float or astropy.time.Time, optional + ecl : str, optional + Obliquity (IERS2010 by default) + + Returns + ------- + np.ndarray : + (len(epoch), 3) array of unit vectors + """ + # TODO: would it be better for this to return a 6-vector (pos, vel)? + + # this was somewhat slow, since it repeatedly created difference SkyCoord Objects + # return self.coords_as_ICRS(epoch=epoch).cartesian.xyz.transpose() + + # Instead look at what https://docs.astropy.org/en/stable/_modules/astropy/coordinates/sky_coordinate.html#SkyCoord.apply_space_motion + # does, which is to use https://github.com/liberfa/erfa/blob/master/src/starpm.c + # and then just use the relevant pieces of that + if epoch is None: + return self.coords_as_ECL(epoch=epoch, ecl=ecl).cartesian.xyz.transpose() + epoch = ( + epoch + if isinstance(epoch, Time) + else Time(epoch, scale="tdb", format="pulsar_mjd") + ) + # compared to the general case above we can assume that the coordinates are ECL + # so just access those components + with warnings.catch_warnings(): + warnings.simplefilter("ignore", ErfaWarning) + # note that pmsafe wants mu_lon not mu_lon * cos(lat) + starpm = pmsafe( + self.ELONG.quantity.to_value(u.radian), + self.ELAT.quantity.to_value(u.radian), + self.PMELONG.quantity.to_value(u.radian / u.yr) + / np.cos(self.ELAT.quantity).value, + self.PMELAT.quantity.to_value(u.radian / u.yr), + 0.0, + 0.0, + self.POSEPOCH.quantity.jd1, + self.POSEPOCH.quantity.jd2, + epoch.jd1, + epoch.jd2, + ) + # lon,lat now in radians + lon, lat = starpm[0], starpm[1] + x = np.cos(lon) * np.cos(lat) + y = np.sin(lon) * np.cos(lat) + z = np.sin(lat) + return u.Quantity([x, y, z]).T + def get_d_delay_quantities_ecliptical(self, toas): """Calculate values needed for many d_delay_d_param functions.""" # TODO: Move all these calculations in a separate class for elegance @@ -904,7 +1073,7 @@ def change_posepoch(self, new_epoch): if isinstance(new_epoch, Time): new_epoch = Time(new_epoch, scale="tdb", precision=9) else: - new_epoch = Time(new_epoch, scale="tdb", format="mjd", precision=9) + new_epoch = Time(new_epoch, scale="tdb", format="pulsar_mjd", precision=9) if self.POSEPOCH.value is None: raise ValueError("POSEPOCH is not currently set.") diff --git a/tests/test_astrometry.py b/tests/test_astrometry.py index 152cd7c27..241875fa5 100644 --- a/tests/test_astrometry.py +++ b/tests/test_astrometry.py @@ -81,3 +81,83 @@ def test_pm_acquires_posepoch_equatorial(): def test_pm_one_set_other_not(): m = get_model(StringIO("\n".join([par_basic_equatorial, "PMRA 7"]))) assert m.POSEPOCH.quantity == m.PEPOCH.quantity + + +parICRS = """PSR 1748-2021E +RAJ 17:48:52.75 1 +DECJ -20:21:29.0 1 +PMRA 1000 +PMDEC -500 +F0 61.485476554 1 +F1 -1.181D-15 1 +PEPOCH 53750.000000 +POSEPOCH 53750.000000 +DM 223.9 1 +SOLARN0 0.00 +EPHEM DE421 +CLK UTC(NIST) +UNITS TDB +TIMEEPH FB90 +T2CMETHOD TEMPO +CORRECT_TROPOSPHERE N +PLANET_SHAPIRO N +DILATEFREQ N +TZRMJD 53801.38605118223 +TZRFRQ 1949.609 +TZRSITE 1 +""" + + +def test_ssb_accuracy_ICRS(): + m = get_model(StringIO(parICRS)) + t0 = m.POSEPOCH.quantity + t0.format = "pulsar_mjd" + t = t0 + np.linspace(0, 20) * u.yr + xyzastropy = m.coords_as_ICRS(epoch=t).cartesian.xyz.transpose() + xyznew = m.ssb_to_psb_xyz_ICRS(epoch=t) + assert np.allclose(xyznew, xyzastropy) + + +parECL = """PSR 1748-2021E +ELONG 300.0 +ELAT -55.0 +PMELONG 1000 +PMELAT -500 +F0 61.485476554 1 +F1 -1.181D-15 1 +PEPOCH 53750.000000 +POSEPOCH 53750.000000 +DM 223.9 1 +SOLARN0 0.00 +EPHEM DE421 +CLK UTC(NIST) +UNITS TDB +TIMEEPH FB90 +T2CMETHOD TEMPO +CORRECT_TROPOSPHERE N +PLANET_SHAPIRO N +DILATEFREQ N +TZRMJD 53801.38605118223 +TZRFRQ 1949.609 +TZRSITE 1 +""" + + +def test_ssb_accuracy_ICRS_ECLmodel(): + m = get_model(StringIO(parECL)) + t0 = m.POSEPOCH.quantity + t0.format = "pulsar_mjd" + t = t0 + np.linspace(0, 20) * u.yr + xyzastropy = m.coords_as_ICRS(epoch=t).cartesian.xyz.transpose() + xyznew = m.ssb_to_psb_xyz_ICRS(epoch=t) + assert np.allclose(xyznew, xyzastropy) + + +def test_ssb_accuracy_ECL_ECLmodel(): + m = get_model(StringIO(parECL)) + t0 = m.POSEPOCH.quantity + t0.format = "pulsar_mjd" + t = t0 + np.linspace(0, 20) * u.yr + xyzastropy = m.coords_as_ECL(epoch=t).cartesian.xyz.transpose() + xyznew = m.ssb_to_psb_xyz_ECL(epoch=t) + assert np.allclose(xyznew, xyzastropy) From 58ffbc24cda83ade56581e2f2fba559dc7382639 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Fri, 29 Sep 2023 10:24:58 -0500 Subject: [PATCH 02/21] checks for obliquity in ecliptic calculations --- src/pint/models/astrometry.py | 22 +++++++++++++++++----- 1 file changed, 17 insertions(+), 5 deletions(-) diff --git a/src/pint/models/astrometry.py b/src/pint/models/astrometry.py index 8867ffe6b..fad389cf1 100644 --- a/src/pint/models/astrometry.py +++ b/src/pint/models/astrometry.py @@ -891,6 +891,11 @@ def ssb_to_psb_xyz_ECL(self, epoch=None, ecl=None): # Instead look at what https://docs.astropy.org/en/stable/_modules/astropy/coordinates/sky_coordinate.html#SkyCoord.apply_space_motion # does, which is to use https://github.com/liberfa/erfa/blob/master/src/starpm.c # and then just use the relevant pieces of that + + # but we need to check that the obliquity is the same + if ecl is not None and ecl != self.ECL.quantity: + return super().ssb_to_psb_xyz_ECL(epoch=epoch, ecl=ecl) + if epoch is None: return self.coords_as_ECL(epoch=epoch, ecl=ecl).cartesian.xyz.transpose() epoch = ( @@ -900,15 +905,22 @@ def ssb_to_psb_xyz_ECL(self, epoch=None, ecl=None): ) # compared to the general case above we can assume that the coordinates are ECL # so just access those components + lon = self.ELONG.quantity.to_value(u.radian) + lat = self.ELAT.quantity.to_value(u.radian) + pm_lon = ( + self.PMELONG.quantity.to_value(u.radian / u.yr) + / np.cos(self.ELAT.quantity).value + ) + pm_lat = self.PMELAT.quantity.to_value(u.radian / u.yr) + with warnings.catch_warnings(): warnings.simplefilter("ignore", ErfaWarning) # note that pmsafe wants mu_lon not mu_lon * cos(lat) starpm = pmsafe( - self.ELONG.quantity.to_value(u.radian), - self.ELAT.quantity.to_value(u.radian), - self.PMELONG.quantity.to_value(u.radian / u.yr) - / np.cos(self.ELAT.quantity).value, - self.PMELAT.quantity.to_value(u.radian / u.yr), + lon, + lat, + pm_lon, + pm_lat, 0.0, 0.0, self.POSEPOCH.quantity.jd1, From c1f4b18e15bb23c22de7e24f8293275a8565f8f3 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Fri, 29 Sep 2023 10:56:31 -0500 Subject: [PATCH 03/21] revert pulsar_mjd -> mjd --- src/pint/models/astrometry.py | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/src/pint/models/astrometry.py b/src/pint/models/astrometry.py index fad389cf1..00d252993 100644 --- a/src/pint/models/astrometry.py +++ b/src/pint/models/astrometry.py @@ -401,9 +401,7 @@ def get_psr_coords(self, epoch=None): frame=coords.ICRS, ) newepoch = ( - epoch - if isinstance(epoch, Time) - else Time(epoch, scale="tdb", format="pulsar_mjd") + epoch if isinstance(epoch, Time) else Time(epoch, scale="tdb", format="mjd") ) position_now = add_dummy_distance(self.get_psr_coords()) with warnings.catch_warnings(): @@ -597,7 +595,7 @@ def change_posepoch(self, new_epoch): if isinstance(new_epoch, Time): new_epoch = Time(new_epoch, scale="tdb", precision=9) else: - new_epoch = Time(new_epoch, scale="tdb", format="pulsar_mjd", precision=9) + new_epoch = Time(new_epoch, scale="tdb", format="mjd", precision=9) if self.POSEPOCH.value is None: raise ValueError("POSEPOCH is not currently set.") @@ -836,9 +834,7 @@ def get_psr_coords(self, epoch=None): ) # Compute for each time because there is proper motion newepoch = ( - epoch - if isinstance(epoch, Time) - else Time(epoch, scale="tdb", format="pulsar_mjd") + epoch if isinstance(epoch, Time) else Time(epoch, scale="tdb", format="mjd") ) position_now = add_dummy_distance(self.get_psr_coords()) return remove_dummy_distance( @@ -1085,7 +1081,7 @@ def change_posepoch(self, new_epoch): if isinstance(new_epoch, Time): new_epoch = Time(new_epoch, scale="tdb", precision=9) else: - new_epoch = Time(new_epoch, scale="tdb", format="pulsar_mjd", precision=9) + new_epoch = Time(new_epoch, scale="tdb", format="mjd", precision=9) if self.POSEPOCH.value is None: raise ValueError("POSEPOCH is not currently set.") From 137bdb49e44504e813826e044bb481408c1bfbba Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Fri, 29 Sep 2023 11:00:55 -0500 Subject: [PATCH 04/21] revert pulsar_mjd -> mjd --- src/pint/models/astrometry.py | 12 +++--------- 1 file changed, 3 insertions(+), 9 deletions(-) diff --git a/src/pint/models/astrometry.py b/src/pint/models/astrometry.py index 00d252993..3fd8d1c49 100644 --- a/src/pint/models/astrometry.py +++ b/src/pint/models/astrometry.py @@ -82,9 +82,7 @@ def ssb_to_psb_xyz_ICRS(self, epoch=None): if epoch is None: return self.coords_as_ICRS(epoch=epoch).cartesian.xyz.transpose() epoch = ( - epoch - if isinstance(epoch, Time) - else Time(epoch, scale="tdb", format="pulsar_mjd") + epoch if isinstance(epoch, Time) else Time(epoch, scale="tdb", format="mjd") ) # in the general case we don't know what system the coordinates will be in # so explicitly transform to ICRS @@ -470,9 +468,7 @@ def ssb_to_psb_xyz_ICRS(self, epoch=None): if epoch is None: return self.coords_as_ICRS(epoch=epoch).cartesian.xyz.transpose() epoch = ( - epoch - if isinstance(epoch, Time) - else Time(epoch, scale="tdb", format="pulsar_mjd") + epoch if isinstance(epoch, Time) else Time(epoch, scale="tdb", format="mjd") ) # compared to the general case above we can assume that the coordinates are ICRS # so just access those components @@ -895,9 +891,7 @@ def ssb_to_psb_xyz_ECL(self, epoch=None, ecl=None): if epoch is None: return self.coords_as_ECL(epoch=epoch, ecl=ecl).cartesian.xyz.transpose() epoch = ( - epoch - if isinstance(epoch, Time) - else Time(epoch, scale="tdb", format="pulsar_mjd") + epoch if isinstance(epoch, Time) else Time(epoch, scale="tdb", format="mjd") ) # compared to the general case above we can assume that the coordinates are ECL # so just access those components From c81591a026a337e33c2623b1e2f6cebe16cf9e2d Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Fri, 29 Sep 2023 11:46:41 -0500 Subject: [PATCH 05/21] fix when no POSEPOCH defined --- src/pint/models/astrometry.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/pint/models/astrometry.py b/src/pint/models/astrometry.py index 3fd8d1c49..55371fc03 100644 --- a/src/pint/models/astrometry.py +++ b/src/pint/models/astrometry.py @@ -79,7 +79,7 @@ def ssb_to_psb_xyz_ICRS(self, epoch=None): # Instead look at what https://docs.astropy.org/en/stable/_modules/astropy/coordinates/sky_coordinate.html#SkyCoord.apply_space_motion # does, which is to use https://github.com/liberfa/erfa/blob/master/src/starpm.c # and then just use the relevant pieces of that - if epoch is None: + if epoch is None or self.POSEPOCH.quantity is None: return self.coords_as_ICRS(epoch=epoch).cartesian.xyz.transpose() epoch = ( epoch if isinstance(epoch, Time) else Time(epoch, scale="tdb", format="mjd") @@ -465,7 +465,7 @@ def ssb_to_psb_xyz_ICRS(self, epoch=None): # Instead look at what https://docs.astropy.org/en/stable/_modules/astropy/coordinates/sky_coordinate.html#SkyCoord.apply_space_motion # does, which is to use https://github.com/liberfa/erfa/blob/master/src/starpm.c # and then just use the relevant pieces of that - if epoch is None: + if epoch is None or self.POSEPOCH is None: return self.coords_as_ICRS(epoch=epoch).cartesian.xyz.transpose() epoch = ( epoch if isinstance(epoch, Time) else Time(epoch, scale="tdb", format="mjd") @@ -888,7 +888,7 @@ def ssb_to_psb_xyz_ECL(self, epoch=None, ecl=None): if ecl is not None and ecl != self.ECL.quantity: return super().ssb_to_psb_xyz_ECL(epoch=epoch, ecl=ecl) - if epoch is None: + if epoch is None or self.POSEPOCH is None: return self.coords_as_ECL(epoch=epoch, ecl=ecl).cartesian.xyz.transpose() epoch = ( epoch if isinstance(epoch, Time) else Time(epoch, scale="tdb", format="mjd") From 83c98594b3e569d9ee93c430ccb60435e4ca49b6 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Fri, 29 Sep 2023 12:58:05 -0500 Subject: [PATCH 06/21] fix when no POSEPOCH defined --- src/pint/models/astrometry.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pint/models/astrometry.py b/src/pint/models/astrometry.py index 55371fc03..2a76c7bf3 100644 --- a/src/pint/models/astrometry.py +++ b/src/pint/models/astrometry.py @@ -465,7 +465,7 @@ def ssb_to_psb_xyz_ICRS(self, epoch=None): # Instead look at what https://docs.astropy.org/en/stable/_modules/astropy/coordinates/sky_coordinate.html#SkyCoord.apply_space_motion # does, which is to use https://github.com/liberfa/erfa/blob/master/src/starpm.c # and then just use the relevant pieces of that - if epoch is None or self.POSEPOCH is None: + if epoch is None or self.POSEPOCH.quantity is None: return self.coords_as_ICRS(epoch=epoch).cartesian.xyz.transpose() epoch = ( epoch if isinstance(epoch, Time) else Time(epoch, scale="tdb", format="mjd") @@ -888,7 +888,7 @@ def ssb_to_psb_xyz_ECL(self, epoch=None, ecl=None): if ecl is not None and ecl != self.ECL.quantity: return super().ssb_to_psb_xyz_ECL(epoch=epoch, ecl=ecl) - if epoch is None or self.POSEPOCH is None: + if epoch is None or self.POSEPOCH.quantity is None: return self.coords_as_ECL(epoch=epoch, ecl=ecl).cartesian.xyz.transpose() epoch = ( epoch if isinstance(epoch, Time) else Time(epoch, scale="tdb", format="mjd") From 4a2d1387d864fdc489351a654f85c4c6cacfa94e Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Fri, 29 Sep 2023 16:47:33 -0500 Subject: [PATCH 07/21] result value to value not quantity --- src/pint/models/timing_model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index 657319fbf..113ff7fa1 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -1885,7 +1885,7 @@ def d_phase_d_param_num(self, toas, param, step=1e-2): res_f = -phase_f[:, 0] + phase_f[:, 1] result = (res_i + res_f) / (2.0 * h * unit) # shift value back to the original value - par.quantity = ori_value + par.value = ori_value return result def d_delay_d_param_num(self, toas, param, step=1e-2): From 2adcdd793aa136b11b2be39c6806dcbb161c6ce9 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Fri, 29 Sep 2023 16:48:07 -0500 Subject: [PATCH 08/21] remove second pytest import --- tests/test_B1855.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/test_B1855.py b/tests/test_B1855.py index e1248e755..605b6953a 100644 --- a/tests/test_B1855.py +++ b/tests/test_B1855.py @@ -1,7 +1,6 @@ """Various tests to assess the performance of the B1855+09.""" import logging import os -import pytest import astropy.units as u import numpy as np From 48f8e732eb0809b88cff91677534718e286adea6 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Fri, 29 Sep 2023 16:48:37 -0500 Subject: [PATCH 09/21] fix some comments, use jd1/jd2 without going through Time conversion when possible --- src/pint/models/astrometry.py | 51 ++++++++++++++++++++--------------- 1 file changed, 30 insertions(+), 21 deletions(-) diff --git a/src/pint/models/astrometry.py b/src/pint/models/astrometry.py index 2a76c7bf3..7219e61c0 100644 --- a/src/pint/models/astrometry.py +++ b/src/pint/models/astrometry.py @@ -73,7 +73,7 @@ def ssb_to_psb_xyz_ICRS(self, epoch=None): """ # TODO: would it be better for this to return a 6-vector (pos, vel)? - # this is somewhat slow, since it repeatedly created difference SkyCoord Objects + # this is somewhat slow, since it repeatedly created different SkyCoord Objects # return self.coords_as_ICRS(epoch=epoch).cartesian.xyz.transpose() # Instead look at what https://docs.astropy.org/en/stable/_modules/astropy/coordinates/sky_coordinate.html#SkyCoord.apply_space_motion @@ -81,10 +81,12 @@ def ssb_to_psb_xyz_ICRS(self, epoch=None): # and then just use the relevant pieces of that if epoch is None or self.POSEPOCH.quantity is None: return self.coords_as_ICRS(epoch=epoch).cartesian.xyz.transpose() - epoch = ( - epoch if isinstance(epoch, Time) else Time(epoch, scale="tdb", format="mjd") - ) - # in the general case we don't know what system the coordinates will be in + if isinstance(epoch, Time): + jd1 = epoch.jd1 + jd2 = epoch.jd2 + else: + jd1 = 2400000.5 + jd2 = epoch # in the general case we don't know what system the coordinates will be in # so explicitly transform to ICRS cc = self.get_psr_coords() icrsrep = cc.icrs.represent_as( @@ -103,8 +105,8 @@ def ssb_to_psb_xyz_ICRS(self, epoch=None): 0.0, self.POSEPOCH.quantity.jd1, self.POSEPOCH.quantity.jd2, - epoch.jd1, - epoch.jd2, + jd1, + jd2, ) # ra,dec now in radians ra, dec = starpm[0], starpm[1] @@ -459,17 +461,21 @@ def ssb_to_psb_xyz_ICRS(self, epoch=None): """ # TODO: would it be better for this to return a 6-vector (pos, vel)? - # this was somewhat slow, since it repeatedly created difference SkyCoord Objects + # this was somewhat slow, since it repeatedly created different SkyCoord Objects # return self.coords_as_ICRS(epoch=epoch).cartesian.xyz.transpose() # Instead look at what https://docs.astropy.org/en/stable/_modules/astropy/coordinates/sky_coordinate.html#SkyCoord.apply_space_motion # does, which is to use https://github.com/liberfa/erfa/blob/master/src/starpm.c # and then just use the relevant pieces of that - if epoch is None or self.POSEPOCH.quantity is None: + if epoch is None or (self.PMRA.quantity == 0 and self.PMDEC.quantity == 0): return self.coords_as_ICRS(epoch=epoch).cartesian.xyz.transpose() - epoch = ( - epoch if isinstance(epoch, Time) else Time(epoch, scale="tdb", format="mjd") - ) + if isinstance(epoch, Time): + jd1 = epoch.jd1 + jd2 = epoch.jd2 + else: + # assume MJD + jd1 = 2400000.5 + jd2 = epoch # compared to the general case above we can assume that the coordinates are ICRS # so just access those components with warnings.catch_warnings(): @@ -485,8 +491,8 @@ def ssb_to_psb_xyz_ICRS(self, epoch=None): 0.0, self.POSEPOCH.quantity.jd1, self.POSEPOCH.quantity.jd2, - epoch.jd1, - epoch.jd2, + jd1, + jd2, ) # ra,dec now in radians ra, dec = starpm[0], starpm[1] @@ -877,7 +883,7 @@ def ssb_to_psb_xyz_ECL(self, epoch=None, ecl=None): """ # TODO: would it be better for this to return a 6-vector (pos, vel)? - # this was somewhat slow, since it repeatedly created difference SkyCoord Objects + # this was somewhat slow, since it repeatedly created different SkyCoord Objects # return self.coords_as_ICRS(epoch=epoch).cartesian.xyz.transpose() # Instead look at what https://docs.astropy.org/en/stable/_modules/astropy/coordinates/sky_coordinate.html#SkyCoord.apply_space_motion @@ -888,11 +894,14 @@ def ssb_to_psb_xyz_ECL(self, epoch=None, ecl=None): if ecl is not None and ecl != self.ECL.quantity: return super().ssb_to_psb_xyz_ECL(epoch=epoch, ecl=ecl) - if epoch is None or self.POSEPOCH.quantity is None: + if epoch is None or (self.PMELONG.value == 0 and self.PMELAT.value == 0): return self.coords_as_ECL(epoch=epoch, ecl=ecl).cartesian.xyz.transpose() - epoch = ( - epoch if isinstance(epoch, Time) else Time(epoch, scale="tdb", format="mjd") - ) + if isinstance(epoch, Time): + jd1 = epoch.jd1 + jd2 = epoch.jd2 + else: + jd1 = 2400000.5 + jd2 = epoch # compared to the general case above we can assume that the coordinates are ECL # so just access those components lon = self.ELONG.quantity.to_value(u.radian) @@ -915,8 +924,8 @@ def ssb_to_psb_xyz_ECL(self, epoch=None, ecl=None): 0.0, self.POSEPOCH.quantity.jd1, self.POSEPOCH.quantity.jd2, - epoch.jd1, - epoch.jd2, + jd1, + jd2, ) # lon,lat now in radians lon, lat = starpm[0], starpm[1] From d887f067c6271da803772d14482d60e72da92c09 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Tue, 3 Oct 2023 10:28:30 -0500 Subject: [PATCH 10/21] test version with multiple implementations --- src/pint/models/astrometry.py | 39 +++++++++++++++++++++++++++++++---- 1 file changed, 35 insertions(+), 4 deletions(-) diff --git a/src/pint/models/astrometry.py b/src/pint/models/astrometry.py index 7219e61c0..f7b42ee13 100644 --- a/src/pint/models/astrometry.py +++ b/src/pint/models/astrometry.py @@ -173,7 +173,7 @@ def sun_angle(self, toas, heliocenter=True, also_distance=False): def barycentric_radio_freq(self, toas): raise NotImplementedError - def solar_system_geometric_delay(self, toas, acc_delay=None): + def solar_system_geometric_delay(self, toas, acc_delay=None, method=1): """Returns geometric delay (in sec) due to position of site in solar system. This includes Roemer delay and parallax. @@ -185,7 +185,9 @@ def solar_system_geometric_delay(self, toas, acc_delay=None): # c selects the non-barycentric TOAs that need actual calculation c = np.logical_and.reduce(tbl["ssb_obs_pos"] != 0, axis=1) if np.any(c): - L_hat = self.ssb_to_psb_xyz_ICRS(epoch=tbl["tdbld"][c].astype(np.float64)) + L_hat = self.ssb_to_psb_xyz_ICRS( + epoch=tbl["tdbld"][c].astype(np.float64), method=method + ) re_dot_L = np.sum(tbl["ssb_obs_pos"][c] * L_hat, axis=1) delay[c] = -re_dot_L.to(ls).value if self.PX.value != 0.0: @@ -445,7 +447,7 @@ def get_params_as_ICRS(self): "PMDEC": self.PMDEC.quantity, } - def ssb_to_psb_xyz_ICRS(self, epoch=None): + def ssb_to_psb_xyz_ICRS(self, epoch=None, method=1): """Returns unit vector(s) from SSB to pulsar system barycenter under ICRS. If epochs (MJD) are given, proper motion is included in the calculation. @@ -467,8 +469,37 @@ def ssb_to_psb_xyz_ICRS(self, epoch=None): # Instead look at what https://docs.astropy.org/en/stable/_modules/astropy/coordinates/sky_coordinate.html#SkyCoord.apply_space_motion # does, which is to use https://github.com/liberfa/erfa/blob/master/src/starpm.c # and then just use the relevant pieces of that - if epoch is None or (self.PMRA.quantity == 0 and self.PMDEC.quantity == 0): + if ( + epoch is None + or (self.PMRA.quantity == 0 and self.PMDEC.quantity == 0) + or method == 0 + ): return self.coords_as_ICRS(epoch=epoch).cartesian.xyz.transpose() + elif method == 1: + ra0, dec0 = self.RAJ.quantity, self.DECJ.quantity + pmra, pmdec = self.PMRA.quantity, self.PMDEC.quantity + posepoch = self.POSEPOCH.quantity + + if epoch is None or (pmra == 0 and pmdec == 0): + ra, dec = ra0, dec0 + else: + epoch = ( + epoch + if isinstance(epoch, Time) + else Time(epoch, scale="tdb", format="mjd") + ) + dt = epoch - posepoch + ra = ra0 + pmra * dt / np.cos(dec0) + dec = dec0 + pmdec * dt + + ra = ra.to_value(u.rad) + dec = dec.to_value(u.rad) + + cos_dec = np.cos(dec) + return u.Quantity( + [cos_dec * np.cos(ra), cos_dec * np.sin(ra), np.sin(dec)] + ).T + if isinstance(epoch, Time): jd1 = epoch.jd1 jd2 = epoch.jd2 From e9ee4abb6168496dbfad1b2d6e249c22cae79ac5 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Tue, 3 Oct 2023 15:04:30 -0500 Subject: [PATCH 11/21] include parallax in position calculation --- src/pint/models/astrometry.py | 67 +++++++++++------------------------ 1 file changed, 20 insertions(+), 47 deletions(-) diff --git a/src/pint/models/astrometry.py b/src/pint/models/astrometry.py index f7b42ee13..7cb55482c 100644 --- a/src/pint/models/astrometry.py +++ b/src/pint/models/astrometry.py @@ -11,7 +11,7 @@ from loguru import logger as log -from erfa import ErfaWarning, pmsafe +from erfa import ErfaWarning, pmsafe, starpm from pint import ls from pint.models.parameter import ( AngleParameter, @@ -23,7 +23,6 @@ from pint.pulsar_ecliptic import OBL, PulsarEcliptic from pint.utils import add_dummy_distance, remove_dummy_distance - astropy_version = sys.modules["astropy"].__version__ mas_yr = u.mas / u.yr @@ -33,6 +32,8 @@ "Astrometry", ] +pmsafe = starpm + class Astrometry(DelayComponent): """Common tools for astrometric calculations.""" @@ -81,12 +82,14 @@ def ssb_to_psb_xyz_ICRS(self, epoch=None): # and then just use the relevant pieces of that if epoch is None or self.POSEPOCH.quantity is None: return self.coords_as_ICRS(epoch=epoch).cartesian.xyz.transpose() + if isinstance(epoch, Time): jd1 = epoch.jd1 jd2 = epoch.jd2 else: jd1 = 2400000.5 - jd2 = epoch # in the general case we don't know what system the coordinates will be in + jd2 = epoch + # in the general case we don't know what system the coordinates will be in # so explicitly transform to ICRS cc = self.get_psr_coords() icrsrep = cc.icrs.represent_as( @@ -96,12 +99,12 @@ def ssb_to_psb_xyz_ICRS(self, epoch=None): icrsvel = icrsrep.differentials["s"] with warnings.catch_warnings(): warnings.simplefilter("ignore", ErfaWarning) - starpm = pmsafe( + starpmout = pmsafe( icrsrep.lon.radian, icrsrep.lat.radian, icrsvel.d_lon.to_value(u.radian / u.yr), icrsvel.d_lat.to_value(u.radian / u.yr), - 0.0, + cc.PX.quantity.to_value(u.arcsec), 0.0, self.POSEPOCH.quantity.jd1, self.POSEPOCH.quantity.jd2, @@ -109,7 +112,7 @@ def ssb_to_psb_xyz_ICRS(self, epoch=None): jd2, ) # ra,dec now in radians - ra, dec = starpm[0], starpm[1] + ra, dec = starpmout[0], starpmout[1] x = np.cos(ra) * np.cos(dec) y = np.sin(ra) * np.cos(dec) z = np.sin(dec) @@ -173,7 +176,7 @@ def sun_angle(self, toas, heliocenter=True, also_distance=False): def barycentric_radio_freq(self, toas): raise NotImplementedError - def solar_system_geometric_delay(self, toas, acc_delay=None, method=1): + def solar_system_geometric_delay(self, toas, acc_delay=None): """Returns geometric delay (in sec) due to position of site in solar system. This includes Roemer delay and parallax. @@ -185,9 +188,7 @@ def solar_system_geometric_delay(self, toas, acc_delay=None, method=1): # c selects the non-barycentric TOAs that need actual calculation c = np.logical_and.reduce(tbl["ssb_obs_pos"] != 0, axis=1) if np.any(c): - L_hat = self.ssb_to_psb_xyz_ICRS( - epoch=tbl["tdbld"][c].astype(np.float64), method=method - ) + L_hat = self.ssb_to_psb_xyz_ICRS(epoch=tbl["tdbld"][c].astype(np.float64)) re_dot_L = np.sum(tbl["ssb_obs_pos"][c] * L_hat, axis=1) delay[c] = -re_dot_L.to(ls).value if self.PX.value != 0.0: @@ -447,7 +448,7 @@ def get_params_as_ICRS(self): "PMDEC": self.PMDEC.quantity, } - def ssb_to_psb_xyz_ICRS(self, epoch=None, method=1): + def ssb_to_psb_xyz_ICRS(self, epoch=None): """Returns unit vector(s) from SSB to pulsar system barycenter under ICRS. If epochs (MJD) are given, proper motion is included in the calculation. @@ -469,36 +470,8 @@ def ssb_to_psb_xyz_ICRS(self, epoch=None, method=1): # Instead look at what https://docs.astropy.org/en/stable/_modules/astropy/coordinates/sky_coordinate.html#SkyCoord.apply_space_motion # does, which is to use https://github.com/liberfa/erfa/blob/master/src/starpm.c # and then just use the relevant pieces of that - if ( - epoch is None - or (self.PMRA.quantity == 0 and self.PMDEC.quantity == 0) - or method == 0 - ): + if epoch is None or (self.PMRA.quantity == 0 and self.PMDEC.quantity == 0): return self.coords_as_ICRS(epoch=epoch).cartesian.xyz.transpose() - elif method == 1: - ra0, dec0 = self.RAJ.quantity, self.DECJ.quantity - pmra, pmdec = self.PMRA.quantity, self.PMDEC.quantity - posepoch = self.POSEPOCH.quantity - - if epoch is None or (pmra == 0 and pmdec == 0): - ra, dec = ra0, dec0 - else: - epoch = ( - epoch - if isinstance(epoch, Time) - else Time(epoch, scale="tdb", format="mjd") - ) - dt = epoch - posepoch - ra = ra0 + pmra * dt / np.cos(dec0) - dec = dec0 + pmdec * dt - - ra = ra.to_value(u.rad) - dec = dec.to_value(u.rad) - - cos_dec = np.cos(dec) - return u.Quantity( - [cos_dec * np.cos(ra), cos_dec * np.sin(ra), np.sin(dec)] - ).T if isinstance(epoch, Time): jd1 = epoch.jd1 @@ -511,14 +484,14 @@ def ssb_to_psb_xyz_ICRS(self, epoch=None, method=1): # so just access those components with warnings.catch_warnings(): warnings.simplefilter("ignore", ErfaWarning) - # note that pmsafe wants mu_alpha not mu_alpha * cos(delta) - starpm = pmsafe( + # note that starpm wants mu_alpha not mu_alpha * cos(delta) + starpmout = pmsafe( self.RAJ.quantity.to_value(u.radian), self.DECJ.quantity.to_value(u.radian), self.PMRA.quantity.to_value(u.radian / u.yr) / np.cos(self.DECJ.quantity).value, self.PMDEC.quantity.to_value(u.radian / u.yr), - 0.0, + self.PX.quantity.to_value(u.arcsec), 0.0, self.POSEPOCH.quantity.jd1, self.POSEPOCH.quantity.jd2, @@ -526,7 +499,7 @@ def ssb_to_psb_xyz_ICRS(self, epoch=None, method=1): jd2, ) # ra,dec now in radians - ra, dec = starpm[0], starpm[1] + ra, dec = starpmout[0], starpmout[1] x = np.cos(ra) * np.cos(dec) y = np.sin(ra) * np.cos(dec) z = np.sin(dec) @@ -946,12 +919,12 @@ def ssb_to_psb_xyz_ECL(self, epoch=None, ecl=None): with warnings.catch_warnings(): warnings.simplefilter("ignore", ErfaWarning) # note that pmsafe wants mu_lon not mu_lon * cos(lat) - starpm = pmsafe( + starpmout = pmsafe( lon, lat, pm_lon, pm_lat, - 0.0, + self.PX.quantity.to_value(u.arcsec), 0.0, self.POSEPOCH.quantity.jd1, self.POSEPOCH.quantity.jd2, @@ -959,7 +932,7 @@ def ssb_to_psb_xyz_ECL(self, epoch=None, ecl=None): jd2, ) # lon,lat now in radians - lon, lat = starpm[0], starpm[1] + lon, lat = starpmout[0], starpmout[1] x = np.cos(lon) * np.cos(lat) y = np.sin(lon) * np.cos(lat) z = np.sin(lat) From 8bd2640bd7efbe09b14f1ce1a5bc3de605b7fcce Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Tue, 3 Oct 2023 15:09:03 -0500 Subject: [PATCH 12/21] fix location of PX --- src/pint/models/astrometry.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pint/models/astrometry.py b/src/pint/models/astrometry.py index 7cb55482c..2b45f71c8 100644 --- a/src/pint/models/astrometry.py +++ b/src/pint/models/astrometry.py @@ -104,7 +104,7 @@ def ssb_to_psb_xyz_ICRS(self, epoch=None): icrsrep.lat.radian, icrsvel.d_lon.to_value(u.radian / u.yr), icrsvel.d_lat.to_value(u.radian / u.yr), - cc.PX.quantity.to_value(u.arcsec), + self.PX.quantity.to_value(u.arcsec), 0.0, self.POSEPOCH.quantity.jd1, self.POSEPOCH.quantity.jd2, From 8ed1bff1ed9a55b68f4cab7de15a724981c90e93 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Wed, 4 Oct 2023 08:42:55 -0500 Subject: [PATCH 13/21] use pmsafe --- src/pint/models/astrometry.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/pint/models/astrometry.py b/src/pint/models/astrometry.py index 2b45f71c8..e08597426 100644 --- a/src/pint/models/astrometry.py +++ b/src/pint/models/astrometry.py @@ -11,7 +11,7 @@ from loguru import logger as log -from erfa import ErfaWarning, pmsafe, starpm +from erfa import ErfaWarning, pmsafe from pint import ls from pint.models.parameter import ( AngleParameter, @@ -32,8 +32,6 @@ "Astrometry", ] -pmsafe = starpm - class Astrometry(DelayComponent): """Common tools for astrometric calculations.""" From 8c5fb9f025ef7c1f6c3b87d1d1e5b9be0fbc4bef Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Thu, 5 Oct 2023 08:53:45 -0500 Subject: [PATCH 14/21] changed numerical derivative steps for PMRA, PMDEC --- tests/test_derivative_utils.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/tests/test_derivative_utils.py b/tests/test_derivative_utils.py index ca35028ae..650caa37c 100644 --- a/tests/test_derivative_utils.py +++ b/tests/test_derivative_utils.py @@ -68,10 +68,8 @@ def get_derivative_params(model): h = 2e-1 elif p in ["FD2"]: h = 3e-2 - elif p in ["F1", "M2", "PMELONG", "PMELAT"]: + elif p in ["F1", "M2", "PMELONG", "PMELAT", "PMRA", "PMDEC"]: h = 2 - elif p in ["PMDEC"]: - h = 3e-2 elif p in ["PBDOT"]: h = 200 elif p in ["FB1"]: From 983c96f4f4263c82d4f164fbd6b39f0a79682159 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Thu, 5 Oct 2023 08:54:00 -0500 Subject: [PATCH 15/21] simplified, cleaned up --- src/pint/models/astrometry.py | 49 +++++------------------------------ 1 file changed, 6 insertions(+), 43 deletions(-) diff --git a/src/pint/models/astrometry.py b/src/pint/models/astrometry.py index e08597426..2d739a13d 100644 --- a/src/pint/models/astrometry.py +++ b/src/pint/models/astrometry.py @@ -73,48 +73,8 @@ def ssb_to_psb_xyz_ICRS(self, epoch=None): # TODO: would it be better for this to return a 6-vector (pos, vel)? # this is somewhat slow, since it repeatedly created different SkyCoord Objects - # return self.coords_as_ICRS(epoch=epoch).cartesian.xyz.transpose() - - # Instead look at what https://docs.astropy.org/en/stable/_modules/astropy/coordinates/sky_coordinate.html#SkyCoord.apply_space_motion - # does, which is to use https://github.com/liberfa/erfa/blob/master/src/starpm.c - # and then just use the relevant pieces of that - if epoch is None or self.POSEPOCH.quantity is None: - return self.coords_as_ICRS(epoch=epoch).cartesian.xyz.transpose() - - if isinstance(epoch, Time): - jd1 = epoch.jd1 - jd2 = epoch.jd2 - else: - jd1 = 2400000.5 - jd2 = epoch - # in the general case we don't know what system the coordinates will be in - # so explicitly transform to ICRS - cc = self.get_psr_coords() - icrsrep = cc.icrs.represent_as( - coords.SphericalRepresentation, - coords.SphericalDifferential, - ) - icrsvel = icrsrep.differentials["s"] - with warnings.catch_warnings(): - warnings.simplefilter("ignore", ErfaWarning) - starpmout = pmsafe( - icrsrep.lon.radian, - icrsrep.lat.radian, - icrsvel.d_lon.to_value(u.radian / u.yr), - icrsvel.d_lat.to_value(u.radian / u.yr), - self.PX.quantity.to_value(u.arcsec), - 0.0, - self.POSEPOCH.quantity.jd1, - self.POSEPOCH.quantity.jd2, - jd1, - jd2, - ) - # ra,dec now in radians - ra, dec = starpmout[0], starpmout[1] - x = np.cos(ra) * np.cos(dec) - y = np.sin(ra) * np.cos(dec) - z = np.sin(dec) - return u.Quantity([x, y, z]).T + # but for consistency only change the method in the subclasses below + return self.coords_as_ICRS(epoch=epoch).cartesian.xyz.transpose() def ssb_to_psb_xyz_ECL(self, epoch=None, ecl=None): """Returns unit vector(s) from SSB to pulsar system barycenter under Ecliptic coordinates. @@ -427,7 +387,7 @@ def coords_as_ECL(self, epoch=None, ecl=None): If `ecl` is left unspecified, the global default IERS2010 will be used. """ if ecl is None: - log.info("ECL not specified; using IERS2010.") + log.debug("ECL not specified; using IERS2010.") ecl = "IERS2010" pos_icrs = self.get_psr_coords(epoch=epoch) @@ -896,6 +856,9 @@ def ssb_to_psb_xyz_ECL(self, epoch=None, ecl=None): if ecl is not None and ecl != self.ECL.quantity: return super().ssb_to_psb_xyz_ECL(epoch=epoch, ecl=ecl) + if ecl is None: + log.debug("ECL not specified; using IERS2010.") + ecl = "IERS2010" if epoch is None or (self.PMELONG.value == 0 and self.PMELAT.value == 0): return self.coords_as_ECL(epoch=epoch, ecl=ecl).cartesian.xyz.transpose() if isinstance(epoch, Time): From b5636fb5a2738fb49588ab3f83438af1bd4aa8a3 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Thu, 5 Oct 2023 12:24:27 -0500 Subject: [PATCH 16/21] more explicit tolerance on test --- tests/test_modelconversions.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/tests/test_modelconversions.py b/tests/test_modelconversions.py index 99b328e63..ffcccdb44 100644 --- a/tests/test_modelconversions.py +++ b/tests/test_modelconversions.py @@ -155,7 +155,10 @@ def test_ICRS_to_ECL_uncertainties(): m2 = fit_ECL.model.as_ICRS() for p in ("RAJ", "DECJ", "PMRA", "PMDEC"): - assert np.isclose(m1.__getitem__(p).value, m2.__getitem__(p).value) + tol = 1e-12 if p in ["RAJ", "DECJ"] else 1e-6 + assert np.isclose( + m1.__getitem__(p).value, m2.__getitem__(p).value, atol=tol + ), f"Paramter {p} with values {m1.__getitem__(p).value}, {m2.__getitem__(p).value} was too far apart (difference={m1.__getitem__(p).value-m2.__getitem__(p).value})" # do a generous test here since the uncertainties could change assert np.isclose( m1.__getitem__(p).uncertainty, m2.__getitem__(p).uncertainty, rtol=0.5 From 0a00cc04aa224c81e365bffa71c2f05052aadd7d Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Thu, 5 Oct 2023 14:56:25 -0500 Subject: [PATCH 17/21] more explicit tolerance on test --- tests/test_modelconversions.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_modelconversions.py b/tests/test_modelconversions.py index ffcccdb44..34a65f0d1 100644 --- a/tests/test_modelconversions.py +++ b/tests/test_modelconversions.py @@ -155,7 +155,7 @@ def test_ICRS_to_ECL_uncertainties(): m2 = fit_ECL.model.as_ICRS() for p in ("RAJ", "DECJ", "PMRA", "PMDEC"): - tol = 1e-12 if p in ["RAJ", "DECJ"] else 1e-6 + tol = 1e-12 if p in ["RAJ", "DECJ"] else 1e-3 assert np.isclose( m1.__getitem__(p).value, m2.__getitem__(p).value, atol=tol ), f"Paramter {p} with values {m1.__getitem__(p).value}, {m2.__getitem__(p).value} was too far apart (difference={m1.__getitem__(p).value-m2.__getitem__(p).value})" From 7b3cb57265ec5a10563bea10a9e03e893160a73f Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Fri, 6 Oct 2023 09:03:59 -0500 Subject: [PATCH 18/21] broader testing --- tests/test_astrometry.py | 52 +++++++++++++++++++++++++++++++++++----- 1 file changed, 46 insertions(+), 6 deletions(-) diff --git a/tests/test_astrometry.py b/tests/test_astrometry.py index 241875fa5..73ff00ecb 100644 --- a/tests/test_astrometry.py +++ b/tests/test_astrometry.py @@ -108,8 +108,22 @@ def test_pm_one_set_other_not(): """ -def test_ssb_accuracy_ICRS(): - m = get_model(StringIO(parICRS)) +ntests = 200 + + +@pytest.mark.parametrize( + ("ra,dec,pmra,pmdec"), + list( + zip( + np.random.uniform(0, 360, ntests), + np.random.uniform(-90, 90, ntests), + np.random.normal(0, scale=10, size=ntests), + np.random.normal(0, scale=10, size=ntests), + ) + ), +) +def test_ssb_accuracy_ICRS(ra, dec, pmra, pmdec): + m = get_model(StringIO(parICRS), RAJ=ra, DECJ=dec, PMRA=pmra, PMDEC=pmdec) t0 = m.POSEPOCH.quantity t0.format = "pulsar_mjd" t = t0 + np.linspace(0, 20) * u.yr @@ -143,8 +157,21 @@ def test_ssb_accuracy_ICRS(): """ -def test_ssb_accuracy_ICRS_ECLmodel(): - m = get_model(StringIO(parECL)) +@pytest.mark.parametrize( + ("elong,elat,pmelong,pmelat"), + list( + zip( + np.random.uniform(0, 360, ntests), + np.random.uniform(-90, 90, ntests), + np.random.normal(0, scale=10, size=ntests), + np.random.normal(0, scale=10, size=ntests), + ) + ), +) +def test_ssb_accuracy_ICRS_ECLmodel(elong, elat, pmelong, pmelat): + m = get_model( + StringIO(parECL), ELONG=elong, ELAT=elat, PMELONG=pmelong, PMELAT=pmelat + ) t0 = m.POSEPOCH.quantity t0.format = "pulsar_mjd" t = t0 + np.linspace(0, 20) * u.yr @@ -153,8 +180,21 @@ def test_ssb_accuracy_ICRS_ECLmodel(): assert np.allclose(xyznew, xyzastropy) -def test_ssb_accuracy_ECL_ECLmodel(): - m = get_model(StringIO(parECL)) +@pytest.mark.parametrize( + ("elong,elat,pmelong,pmelat"), + list( + zip( + np.random.uniform(0, 360, ntests), + np.random.uniform(-90, 90, ntests), + np.random.normal(0, scale=10, size=ntests), + np.random.normal(0, scale=10, size=ntests), + ) + ), +) +def test_ssb_accuracy_ECL_ECLmodel(elong, elat, pmelong, pmelat): + m = get_model( + StringIO(parECL), ELONG=elong, ELAT=elat, PMELONG=pmelong, PMELAT=pmelat + ) t0 = m.POSEPOCH.quantity t0.format = "pulsar_mjd" t = t0 + np.linspace(0, 20) * u.yr From 5939d9b4e877850fb35977981e48b6992ea3530b Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Fri, 13 Oct 2023 10:27:11 -0500 Subject: [PATCH 19/21] remove duplicate parallax --- src/pint/models/astrometry.py | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/src/pint/models/astrometry.py b/src/pint/models/astrometry.py index 2d739a13d..bc5aeabf7 100644 --- a/src/pint/models/astrometry.py +++ b/src/pint/models/astrometry.py @@ -149,16 +149,6 @@ def solar_system_geometric_delay(self, toas, acc_delay=None): L_hat = self.ssb_to_psb_xyz_ICRS(epoch=tbl["tdbld"][c].astype(np.float64)) re_dot_L = np.sum(tbl["ssb_obs_pos"][c] * L_hat, axis=1) delay[c] = -re_dot_L.to(ls).value - if self.PX.value != 0.0: - L = (1.0 / self.PX.value) * u.kpc - # TODO: np.sum currently loses units in some cases... - re_sqr = ( - np.sum(tbl["ssb_obs_pos"][c] ** 2, axis=1) - * tbl["ssb_obs_pos"].unit ** 2 - ) - delay[c] += ( - (0.5 * (re_sqr / L) * (1.0 - re_dot_L**2 / re_sqr)).to(ls).value - ) return delay * u.second def get_d_delay_quantities(self, toas): From fe769b2a70ce1be1c41ce4019c6b29a11b245831 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Fri, 13 Oct 2023 10:38:22 -0500 Subject: [PATCH 20/21] put parallax back in only 1 place? --- src/pint/models/astrometry.py | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/src/pint/models/astrometry.py b/src/pint/models/astrometry.py index bc5aeabf7..7a37e1203 100644 --- a/src/pint/models/astrometry.py +++ b/src/pint/models/astrometry.py @@ -149,6 +149,16 @@ def solar_system_geometric_delay(self, toas, acc_delay=None): L_hat = self.ssb_to_psb_xyz_ICRS(epoch=tbl["tdbld"][c].astype(np.float64)) re_dot_L = np.sum(tbl["ssb_obs_pos"][c] * L_hat, axis=1) delay[c] = -re_dot_L.to(ls).value + if self.PX.value != 0.0: + L = (1.0 / self.PX.value) * u.kpc + # TODO: np.sum currently loses units in some cases... + re_sqr = ( + np.sum(tbl["ssb_obs_pos"][c] ** 2, axis=1) + * tbl["ssb_obs_pos"].unit ** 2 + ) + delay[c] += ( + (0.5 * (re_sqr / L) * (1.0 - re_dot_L**2 / re_sqr)).to(ls).value + ) return delay * u.second def get_d_delay_quantities(self, toas): @@ -439,7 +449,8 @@ def ssb_to_psb_xyz_ICRS(self, epoch=None): self.PMRA.quantity.to_value(u.radian / u.yr) / np.cos(self.DECJ.quantity).value, self.PMDEC.quantity.to_value(u.radian / u.yr), - self.PX.quantity.to_value(u.arcsec), + # self.PX.quantity.to_value(u.arcsec), + 0.0, 0.0, self.POSEPOCH.quantity.jd1, self.POSEPOCH.quantity.jd2, @@ -875,7 +886,8 @@ def ssb_to_psb_xyz_ECL(self, epoch=None, ecl=None): lat, pm_lon, pm_lat, - self.PX.quantity.to_value(u.arcsec), + # self.PX.quantity.to_value(u.arcsec), + 0.0, 0.0, self.POSEPOCH.quantity.jd1, self.POSEPOCH.quantity.jd2, From deb3ec4c1114e950ca47754a1ee0e60c556b52d2 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Fri, 13 Oct 2023 12:40:11 -0500 Subject: [PATCH 21/21] parallax back in PMSAFE and explicitly in delay --- src/pint/models/astrometry.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/pint/models/astrometry.py b/src/pint/models/astrometry.py index 7a37e1203..2d739a13d 100644 --- a/src/pint/models/astrometry.py +++ b/src/pint/models/astrometry.py @@ -449,8 +449,7 @@ def ssb_to_psb_xyz_ICRS(self, epoch=None): self.PMRA.quantity.to_value(u.radian / u.yr) / np.cos(self.DECJ.quantity).value, self.PMDEC.quantity.to_value(u.radian / u.yr), - # self.PX.quantity.to_value(u.arcsec), - 0.0, + self.PX.quantity.to_value(u.arcsec), 0.0, self.POSEPOCH.quantity.jd1, self.POSEPOCH.quantity.jd2, @@ -886,8 +885,7 @@ def ssb_to_psb_xyz_ECL(self, epoch=None, ecl=None): lat, pm_lon, pm_lat, - # self.PX.quantity.to_value(u.arcsec), - 0.0, + self.PX.quantity.to_value(u.arcsec), 0.0, self.POSEPOCH.quantity.jd1, self.POSEPOCH.quantity.jd2,