From df4fb7ff51b3234d1f341359adf14978d29ae271 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Wed, 17 Dec 2025 16:26:28 +0530 Subject: [PATCH 1/5] vlbi_coord_rotation --- src/pint/models/astrometry.py | 57 +++++++++++++++++++++++++++++++++-- 1 file changed, 55 insertions(+), 2 deletions(-) diff --git a/src/pint/models/astrometry.py b/src/pint/models/astrometry.py index c08a6fd49..2a2ff558f 100644 --- a/src/pint/models/astrometry.py +++ b/src/pint/models/astrometry.py @@ -65,6 +65,36 @@ def __init__(self): ) ) + self.add_param( + floatParameter( + name="VLBIAX", + units="mas", + value=None, + description="X component of the offset between VLBI and pulsar timing coordinate systems", + tcb2tdb_scale_factor=1, + ) + ) + + self.add_param( + floatParameter( + name="VLBIAY", + units="mas", + value=None, + description="Y component of the offset between VLBI and pulsar timing coordinate systems", + tcb2tdb_scale_factor=1, + ) + ) + + self.add_param( + floatParameter( + name="VLBIAZ", + units="mas", + value=None, + description="Z component of the offset between VLBI and pulsar timing coordinate systems", + tcb2tdb_scale_factor=1, + ) + ) + self.delay_funcs_component += [self.solar_system_geometric_delay] self.register_deriv_funcs(self.d_delay_astrometry_d_PX, "PX") @@ -268,6 +298,25 @@ def as_ECL(self, epoch=None, ecl="IERS2010"): def as_ICRS(self, epoch=None, ecl="IERS2010"): raise NotImplementedError + def vlbi_coord_rotation(self) -> Optional[np.ndarray]: + if ( + self.VLBIAX.quantity is not None + and self.VLBIAY.quantity is not None + and self.VLBIAZ.quantity is not None + ): + Ax = self.VLBIAX.quantity.to_value(u.rad) + Ay = self.VLBIAY.quantity.to_value(u.rad) + Az = self.VLBIAZ.quantity.to_value(u.rad) + return np.array( + [ + [1, Az, -Ay], + [-Az, 1, Ax], + [Ay, -Ax, 1], + ] + ) + else: + return None + class AstrometryEquatorial(Astrometry): """Astrometry in equatorial coordinates. @@ -525,7 +574,9 @@ def ssb_to_psb_xyz_ICRS(self, epoch: Optional[time_like] = None) -> u.Quantity: ) # ra,dec now in radians ra, dec = starpmout[0], starpmout[1] - return self.xyz_from_radec(ra, dec) + Omega = self.vlbi_coord_rotation() + Khat = self.xyz_from_radec(ra, dec) + return Omega @ Khat if Omega is not None else Khat def xyz_from_radec(self, ra, dec): x = np.cos(ra) * np.cos(dec) @@ -1023,7 +1074,9 @@ def ssb_to_psb_xyz_ECL( ) # lon,lat now in radians lon, lat = starpmout[0], starpmout[1] - return self.xyz_from_latlong(lon, lat) + Khat = self.xyz_from_latlong(lon, lat) + Omega = self.vlbi_coord_rotation() + return Omega @ Khat if Omega is not None else Khat def xyz_from_latlong(self, lon, lat): x = np.cos(lon) * np.cos(lat) From 502d3351b332056beeffd65986968a21c0e31535 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Thu, 18 Dec 2025 10:20:27 +0530 Subject: [PATCH 2/5] test_vlbi_rotation --- tests/test_vlbi_rotation.py | 77 +++++++++++++++++++++++++++++++++++++ 1 file changed, 77 insertions(+) create mode 100644 tests/test_vlbi_rotation.py diff --git a/tests/test_vlbi_rotation.py b/tests/test_vlbi_rotation.py new file mode 100644 index 000000000..4e32d8937 --- /dev/null +++ b/tests/test_vlbi_rotation.py @@ -0,0 +1,77 @@ +import io +from pint.models import get_model +from pint.simulation import make_fake_toas_uniform +from pint.residuals import Residuals +from pint.fitter import Fitter +import pytest +import copy +import numpy as np + +pars = [ + """ + RAJ 05:00:00 1 + DECJ 15:00:00 1 + POSEPOCH 55000 + VLBIAX -0.15 + VLBIAY 0.05 + VLBIAZ 0.21 + F0 100 1 + F1 -1e-15 1 + PEPOCH 55000 + TZRMJD 55000 + TZRFRQ 1400 + TZRSITE gbt + UNITS TDB + EPHEM DE440 + """, + """ + ELAT 30.0 1 + ELONG 75.0 1 + POSEPOCH 55000 + VLBIAX -0.15 + VLBIAY 0.05 + VLBIAZ 0.21 + F0 100 1 + F1 -1e-15 1 + PEPOCH 55000 + TZRMJD 55000 + TZRFRQ 1400 + TZRSITE gbt + UNITS TDB + EPHEM DE440 + """, +] + + +@pytest.mark.parametrize("par", pars) +def test_vlbi_rotation(par): + m0 = get_model(io.StringIO(par)) + + assert all(param in m0 for param in ["VLBIAX", "VLBIAY", "VLBIAZ"]) + + t = make_fake_toas_uniform( + startMJD=54000, + endMJD=56000, + ntoas=200, + model=m0, + add_noise=True, + ) + + assert Residuals(t, m0).reduced_chi2 < 1.5 + + ftr = Fitter.auto(t, m0) + ftr.fit_toas() + assert ftr.resids.reduced_chi2 < 1.5 + + # ------------ + + m1 = copy.deepcopy(m0) + m1["VLBIAX"].value = m1["VLBIAY"].value = m1["VLBIAZ"].value = 0 + + ftr1 = Fitter.auto(t, m1) + ftr1.fit_toas() + assert ftr1.resids.reduced_chi2 < 1.5 + + assert np.allclose( + ftr1.model.solar_system_geometric_delay(t), m0.solar_system_geometric_delay(t) + ) From 545d73ce2d8e9caaceb67cdb51dc213b3caf6116 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Thu, 18 Dec 2025 10:21:45 +0530 Subject: [PATCH 3/5] CHANGELOG-unreleased --- CHANGELOG-unreleased.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index c3dbcd584..94c610df9 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -14,5 +14,6 @@ the released changes. - KS test to check if the whitened residuals are unit-normal distributed - Warning about setting of TZRMJD from TOAs - Method to zero out mean residual based on TZRMJD +- Use VLBI astrometric measurements along with coordinate offset in the timing model ### Fixed ### Removed From add71d149fadcc185b9a98c7f6ea6d9a93306955 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Thu, 18 Dec 2025 10:26:25 +0530 Subject: [PATCH 4/5] docs --- src/pint/models/astrometry.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/pint/models/astrometry.py b/src/pint/models/astrometry.py index 2a2ff558f..a35d87de1 100644 --- a/src/pint/models/astrometry.py +++ b/src/pint/models/astrometry.py @@ -299,6 +299,8 @@ def as_ICRS(self, epoch=None, ecl="IERS2010"): raise NotImplementedError def vlbi_coord_rotation(self) -> Optional[np.ndarray]: + """Returns the coordinate rotation matrix between VLBI and pulsar + timing coordinate systems if the rotation parameters are given.""" if ( self.VLBIAX.quantity is not None and self.VLBIAY.quantity is not None From 1cab23a63ca5b23d4cb90ae2664870d8cbf8e001 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Fri, 19 Dec 2025 09:55:39 +0530 Subject: [PATCH 5/5] reference --- src/pint/models/astrometry.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/src/pint/models/astrometry.py b/src/pint/models/astrometry.py index a35d87de1..7a4cd06e1 100644 --- a/src/pint/models/astrometry.py +++ b/src/pint/models/astrometry.py @@ -300,7 +300,11 @@ def as_ICRS(self, epoch=None, ecl="IERS2010"): def vlbi_coord_rotation(self) -> Optional[np.ndarray]: """Returns the coordinate rotation matrix between VLBI and pulsar - timing coordinate systems if the rotation parameters are given.""" + timing coordinate systems if the rotation parameters are given. + + Reference: + Madison+ 2013, The Astrophysical Journal 777 104 (Equation 9) + """ if ( self.VLBIAX.quantity is not None and self.VLBIAY.quantity is not None @@ -576,6 +580,8 @@ def ssb_to_psb_xyz_ICRS(self, epoch: Optional[time_like] = None) -> u.Quantity: ) # ra,dec now in radians ra, dec = starpmout[0], starpmout[1] + + # Reference: Madison+ 2023, The Astrophysical Journal 777 104 (Equations 9, 10) Omega = self.vlbi_coord_rotation() Khat = self.xyz_from_radec(ra, dec) return Omega @ Khat if Omega is not None else Khat @@ -1077,6 +1083,8 @@ def ssb_to_psb_xyz_ECL( # lon,lat now in radians lon, lat = starpmout[0], starpmout[1] Khat = self.xyz_from_latlong(lon, lat) + + # Reference: Madison+ 2023, The Astrophysical Journal 777 104 (Equations 9, 10) Omega = self.vlbi_coord_rotation() return Omega @ Khat if Omega is not None else Khat