Use VLBI astrometric measurements along with coordinate system offset in the timing model

CHANGELOG-unreleased.md

@@ -15,5 +15,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

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,31 @@ 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]:
+        """Returns the coordinate rotation matrix between VLBI and pulsar
+        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
+            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 +580,11 @@ 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)
+
+        # 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
 
     def xyz_from_radec(self, ra, dec):
         x = np.cos(ra) * np.cos(dec)
@@ -1023,7 +1082,11 @@ 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)
+
+        # 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
 
     def xyz_from_latlong(self, lon, lat):
         x = np.cos(lon) * np.cos(lat)

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)
+    )