Merge branch 'master' into clock-search

CHANGELOG.md

@@ -20,6 +20,7 @@ and this project, at least loosely, adheres to [Semantic Versioning](https://sem
 - Observatories using TEMPO-format clock files now refer to the clock file (for example time_ao.dat) rather than via time.dat
 - Observatories that don't use clock corrections (for example CHIME uses GPS time directly) just use an empty list of clock correction files rather than files containing only zeros
 - Updated Jodrell clock corrections to include post-1997
+- DDK model will now use ICRS or ECL coordinates depending on what the input model is
 
 ## [0.8.6 == 0.8.7] 2022-05-10
 ### Added

src/pint/models/astrometry.py

@@ -63,13 +63,19 @@ def ssb_to_psb_xyz_ICRS(self, epoch=None):
         # TODO: would it be better for this to return a 6-vector (pos, vel)?
         return self.coords_as_ICRS(epoch=epoch).cartesian.xyz.transpose()
 
-    def ssb_to_psb_xyz_ECL(self, epoch=None):
+    def ssb_to_psb_xyz_ECL(self, epoch=None, ecl=None):
         """Returns unit vector(s) from SSB to pulsar system barycenter under Ecliptic coordinates.
 
         If epochs (MJD) are given, proper motion is included in the calculation.
+
+        Parameters
+        ----------
+        epoch : float, optional
+        ecl : str, optional
+            Obliquity (IERS2010 by default)
         """
         # TODO: would it be better for this to return a 6-vector (pos, vel)?
-        return self.coords_as_ECL(epoch=epoch).cartesian.xyz.transpose()
+        return self.coords_as_ECL(epoch=epoch, ecl=ecl).cartesian.xyz.transpose()
 
     def sun_angle(self, toas, heliocenter=True, also_distance=False):
         """Compute the pulsar-observatory-Sun angle.
@@ -1011,8 +1017,8 @@ def as_ICRS(self, epoch=None):
             frame=PulsarEcliptic,
         )
         c_ICRS = c.transform_to(coords.ICRS)
-        m_eq.RAJ.uncertainty = c_ICRS.pm_ra_cosdec * dt / np.cos(c_ICRS.dec)
-        m_eq.DECJ.uncertainty = c_ICRS.pm_dec * dt
+        m_eq.RAJ.uncertainty = np.abs(c_ICRS.pm_ra_cosdec * dt / np.cos(c_ICRS.dec))
+        m_eq.DECJ.uncertainty = np.abs(c_ICRS.pm_dec * dt)
         # use fake proper motions to convert uncertainties on proper motion
         # assume that PMELONG uncertainty includes cos(DEC)
         c = coords.SkyCoord(
@@ -1025,8 +1031,8 @@ def as_ICRS(self, epoch=None):
             frame=PulsarEcliptic,
         )
         c_ICRS = c.transform_to(coords.ICRS)
-        m_eq.PMRA.uncertainty = c_ICRS.pm_ra_cosdec
-        m_eq.PMDEC.uncertainty = c_ICRS.pm_dec
+        m_eq.PMRA.uncertainty = np.abs(c_ICRS.pm_ra_cosdec)
+        m_eq.PMDEC.uncertainty = np.abs(c_ICRS.pm_dec)
         # freeze comparable parameters
         m_eq.RAJ.frozen = self.ELONG.frozen
         m_eq.DECJ.frozen = self.ELAT.frozen

src/pint/models/binary_ddk.py

@@ -2,6 +2,8 @@
 import warnings
 import numpy as np
 from astropy import units as u
+from loguru import logger as log
+
 from pint.models.binary_dd import BinaryDD
 from pint.models.parameter import boolParameter, floatParameter
 from pint.models.stand_alone_psr_binaries.DDK_model import DDKmodel
@@ -16,25 +18,36 @@ class BinaryDDK(BinaryDD):
     of the system from Earth, the finite size of the system, and the
     interaction of these with the proper motion.
 
+    From Kopeikin (1995) this includes :math:`\Delta_{\pi M}` (Equation 17), the mixed annual-orbital parallax term, which changes :math:`a_1` and :math:`\omega`
+    (:meth:`~pint.models.stand_alone_psr_binaries.DDK_model.DDKmodel.delta_a1_parallax` and :meth:`~pint.models.stand_alone_psr_binaries.DDK_model.DDKmodel.delta_omega_parallax`).
+
+    It does not include :math:`\Delta_{\pi P}`, the pure pulsar orbital parallax term (Equation 14).
+
+    From Kopeikin (1996) this includes apparent changes in :math:`\omega`, :math:`a_1`, and :math:`i` due to the proper motion
+    (:meth:`~pint.models.stand_alone_psr_binaries.DDK_model.DDKmodel.delta_omega_proper_motion`, :meth:`~pint.models.stand_alone_psr_binaries.DDK_model.DDKmodel.delta_a1_proper_motion`,
+    :meth:`~pint.models.stand_alone_psr_binaries.DDK_model.DDKmodel.delta_kin_proper_motion`) (Equations 8, 9, 10).
+
     The actual calculations for this are done in
     :class:`pint.models.stand_alone_psr_binaries.DDK_model.DDKmodel`.
 
     It supports all the parameters defined in :class:`pint.models.pulsar_binary.PulsarBinary`
-    and :class:`pint.models.pulsar_binary.BinaryDDK` plus:
+    and :class:`pint.models.binary_dd.BinaryDD` plus:
 
-        - KIN - inclination angle (deg)
-        - KOM - the longitude of the ascending node, Kopeikin (1995) Eq 9. OMEGA (deg)
-        - K96 - flag for Kopeikin binary model proper motion correction
+       KIN
+            the inclination angle: :math:`i`
+       KOM
+            the longitude of the ascending node, Kopeikin (1995) Eq 9: :math:`\Omega`
+       K96
+            flag for Kopeikin binary model proper motion correction
 
     It also removes:
 
-        - SINI - use KIN instead
+       SINI
+            use ``KIN`` instead
 
     Note
     ----
-    This model defines KOM with reference to celestial north regardless of the astrometry
-    model. This is incompatible with tempo2, which defines KOM with reference to ecliptic
-    north when using ecliptic coordinates.
+    This model defines KOM with reference to north, either equatorial or ecliptic depending on how the model is defined.
 
     Parameters supported:
 
@@ -43,7 +56,12 @@ class BinaryDDK(BinaryDD):
 
     References
     ----------
-    KOPEIKIN. 1995, 1996
+    - Kopeikin (1995), ApJ, 439, L5 [1]_
+    - Kopeikin (1996), ApJ, 467, L93 [2]_
+
+    .. [1] https://ui.adsabs.harvard.edu/abs/1995ApJ...439L...5K/abstract
+    .. [2] https://ui.adsabs.harvard.edu/abs/1996ApJ...467L..93K/abstract
+
     """
 
     register = True
@@ -76,42 +94,54 @@ def __init__(
             )
         )
         self.remove_param("SINI")
-        self.internal_params += ["PMRA_DDK", "PMDEC_DDK"]
+        self.internal_params += ["PMLONG_DDK", "PMLAT_DDK"]
 
     @property
-    def PMRA_DDK(self):
-        params = self._parent.get_params_as_ICRS()
-        par_obj = floatParameter(
-            name="PMRA",
-            units="mas/year",
-            value=params["PMRA"],
-            description="Proper motion in RA",
-        )
-        return par_obj
+    def PMLONG_DDK(self):
+        """Proper motion in longitude (RA or ecliptic longitude)"""
+        if "AstrometryEquatorial" in self._parent.components:
+            return self._parent.PMRA
+        elif "AstrometryEcliptic" in self._parent.components:
+            return self._parent.PMELONG
+        else:
+            raise TimingModelError(
+                "No valid AstrometryEcliptic or AstrometryEquatorial component found"
+            )
 
     @property
-    def PMDEC_DDK(self):
-        params = self._parent.get_params_as_ICRS()
-        par_obj = floatParameter(
-            name="PMDEC",
-            units="mas/year",
-            value=params["PMDEC"],
-            description="Proper motion in DEC",
-        )
-        return par_obj
+    def PMLAT_DDK(self):
+        """Proper motion in latitude (Dec or ecliptic latitude)"""
+        if "AstrometryEquatorial" in self._parent.components:
+            return self._parent.PMDEC
+        elif "AstrometryEcliptic" in self._parent.components:
+            return self._parent.PMELAT
+        else:
+            raise TimingModelError(
+                "No valid AstrometryEcliptic or AstrometryEquatorial component found"
+            )
 
     def validate(self):
         """Validate parameters."""
         super().validate()
-        if "PMRA" not in self._parent.params or "PMDEC" not in self._parent.params:
-            # Check ecliptic coordinates proper motion.
+        if "AstrometryEquatorial" in self._parent.components:
+            log.debug("Validating DDK model in ICRS coordinates")
+            if "PMRA" not in self._parent.params or "PMDEC" not in self._parent.params:
+                raise MissingParameter(
+                    "DDK", "DDK model needs proper motion parameters."
+                )
+        elif "AstrometryEcliptic" in self._parent.components:
+            log.debug("Validating DDK model in ECL coordinates")
             if (
                 "PMELONG" not in self._parent.params
                 or "PMELAT" not in self._parent.params
             ):
                 raise MissingParameter(
                     "DDK", "DDK model needs proper motion parameters."
                 )
+        else:
+            raise TimingModelError(
+                "No valid AstrometryEcliptic or AstrometryEquatorial component found"
+            )
 
         if hasattr(self._parent, "PX"):
             if self._parent.PX.value <= 0.0 or self._parent.PX.value is None:
@@ -123,16 +153,14 @@ def validate(self):
 
         if "A1DOT" in self.params and self.A1DOT.value != 0:
             warnings.warn("Using A1DOT with a DDK model is not advised.")
-        # Should we warn if the model is using ecliptic coordinates?
-        # Should we support KOM_PINT that works this way and KOM that works the way tempo2 does?
 
     def alternative_solutions(self):
         """Alternative Kopeikin solutions (potential local minima)
 
         There are 4 potential local minima for a DDK model where a1dot is the same
         These are given by where Eqn. 8 in Kopeikin (1996) is equal to the best-fit value.
 
-        We first define the symmetry point where a1dot is zero:
+        We first define the symmetry point where a1dot is zero (in equatorial coordinates):
 
         :math:`KOM_0 = \\tan^{-1} (\mu_{\delta} / \mu_{\\alpha})`
 
@@ -157,7 +185,9 @@ def alternative_solutions(self):
         y0 = self.KOM.quantity
         solutions = [(x0, y0)]
         # where Eqn. 8 in Kopeikin (1996) that is equal to 0
-        KOM_zero = np.arctan2(self.PMDEC_DDK.quantity, self.PMRA_DDK.quantity).to(u.deg)
+        KOM_zero = np.arctan2(self.PMLAT_DDK.quantity, self.PMLONG_DDK.quantity).to(
+            u.deg
+        )
         # second one in the same banana
         solutions += [(x0, (2 * (KOM_zero) - y0 - 180 * u.deg) % (360 * u.deg))]
         # and the other banana

src/pint/models/pulsar_binary.py

@@ -7,13 +7,15 @@
 import astropy.units as u
 import numpy as np
 from astropy.time import Time
+from astropy.coordinates import SkyCoord
 from loguru import logger as log
 
 from pint import ls
 from pint.models.parameter import MJDParameter, floatParameter, prefixParameter
 from pint.models.stand_alone_psr_binaries import binary_orbits as bo
 from pint.models.timing_model import DelayComponent, MissingParameter, UnknownParameter
 from pint.utils import taylor_horner_deriv
+from pint.pulsar_ecliptic import PulsarEcliptic
 
 
 class PulsarBinary(DelayComponent):
@@ -269,6 +271,11 @@ def update_binary_object(self, toas, acc_delay=None):
             the standard barycentering. The stand alone binary receives the
             input TOAs - acc_delay.
 
+        Notes
+        -----
+        The values for ``obs_pos`` (the observatory position wrt the Solar System Barycenter) and ``psr_pos``
+        (the pulsar position wrt the Solar System Barycenter) are both computed in the same reference frame, ICRS or ECL depending on the model.
+
         Warns
         -----
         If passing 'None' to 'toa' argument, the stand alone binary model will use
@@ -291,10 +298,25 @@ def update_binary_object(self, toas, acc_delay=None):
             else:
                 self.barycentric_time = tbl["tdbld"] * u.day - acc_delay
             updates["barycentric_toa"] = self.barycentric_time
-            updates["obs_pos"] = tbl["ssb_obs_pos"].quantity
-            updates["psr_pos"] = self._parent.ssb_to_psb_xyz_ICRS(
-                epoch=tbl["tdbld"].astype(np.float64)
-            )
+            if "AstrometryEquatorial" in self._parent.components:
+                # it's already in ICRS
+                updates["obs_pos"] = tbl["ssb_obs_pos"].quantity
+                updates["psr_pos"] = self._parent.ssb_to_psb_xyz_ICRS(
+                    epoch=tbl["tdbld"].astype(np.float64)
+                )
+            elif "AstrometryEcliptic" in self._parent.components:
+                # convert from ICRS to ECL
+                obs_pos = SkyCoord(
+                    tbl["ssb_obs_pos"].quantity,
+                    representation_type="cartesian",
+                    frame="icrs",
+                )
+                updates["obs_pos"] = obs_pos.transform_to(
+                    PulsarEcliptic(ecl=self._parent.ECL.value)
+                ).cartesian.xyz.transpose()
+                updates["psr_pos"] = self._parent.ssb_to_psb_xyz_ECL(
+                    epoch=tbl["tdbld"].astype(np.float64), ecl=self._parent.ECL.value
+                )
         for par in self.binary_instance.binary_params:
             if par in self.binary_instance.param_aliases.keys():
                 alias = self.binary_instance.param_aliases[par]