Replace decimal year conversions with julian year conversions.

orbitize/gaia.py

@@ -7,7 +7,7 @@
     from astroquery.gaia import Gaia
 from astropy import units as u
 import astropy.io.fits as fits
-import astropy.time as time
+from astropy.time import Time
 from astropy.io.ascii import read
 from astropy.coordinates import get_body_barycentric_posvel
 import numpy.linalg
@@ -61,12 +61,12 @@ def __init__(self, gaia_num, hiplogprob, dr="dr2", query=True, gaia_data=None):
         self.dr = dr
 
         if self.dr == "edr3":
-            self.gaia_epoch = 2016.0
+            self.gaia_epoch = Time("J2016", format='jyear_str')
         elif self.dr == "dr2":
-            self.gaia_epoch = 2015.5
+            self.gaia_epoch = Time("J2015.5", format='jyear_str')
         else:
             raise ValueError("`dr` must be either `dr2` or `edr3`")
-        self.hipparcos_epoch = 1991.25
+        self.hipparcos_epoch = Time("J1991.25", format='jyear_str')
 
         if query:
             query = """SELECT
@@ -134,12 +134,12 @@ def compute_lnlike(self, raoff_model, deoff_model, samples, param_idx):
         """
 
         alpha_H0 = samples[param_idx["alpha0"]]  # [deg]
-        pm_ra = samples[param_idx["pm_ra"]]  # [mas/yr]
-        delta_alpha_from_pm = pm_ra * (self.gaia_epoch - self.hipparcos_epoch)  # [mas]
+        pm_ra = samples[param_idx["pm_ra"]]  # [mas/jyr]
+        delta_alpha_from_pm = pm_ra * (self.gaia_epoch.jyear - self.hipparcos_epoch.jyear)  # [mas]
 
         delta_H0 = samples[param_idx["delta0"]]  # [deg]
-        pm_dec = samples[param_idx["pm_dec"]]  # [mas/yr]
-        delta_delta_from_pm = pm_dec * (self.gaia_epoch - self.hipparcos_epoch)  # [mas]
+        pm_dec = samples[param_idx["pm_dec"]]  # [mas/jyr]
+        delta_delta_from_pm = pm_dec * (self.gaia_epoch.jyear - self.hipparcos_epoch.jyear)  # [mas]
 
         # difference in position due to orbital motion between Hipparcos & Gaia epochs
         alpha_diff_orbit = raoff_model[1, :] - raoff_model[0, :]  # [mas]
@@ -162,6 +162,7 @@ def compute_lnlike(self, raoff_model, deoff_model, samples, param_idx):
         # technically this is an angle so we should wrap it, but the precision
         # of Hipparcos and Gaia is so good that we'll never have to.
         alpha_resid = alpha_model - alpha_data
+
         alpha_chi2 = (alpha_resid / alpha_unc) ** 2
 
         delta_model = self.hiplogprob.delta0 + self.mas2deg * (  # [deg]
@@ -296,9 +297,9 @@ def __init__(self, hip_id, hiplogprob, gost_filepath, hgca_filepath=None):
         self.gaia_epoch_dec = entry["epoch_dec_gaia"][0]
         # read in the GOST file to get the estimated Gaia epochs and scan angles
         gost_dat = read(gost_filepath, converters={"*": [int, float, bytes]})
-        self.gaia_epoch = time.Time(
+        self.gaia_epoch = Time(
             gost_dat["ObservationTimeAtGaia[UTC]"]
-        ).decimalyear  # in decimal year
+        )  # in julian year
         gaia_scan_theta = np.array(gost_dat["scanAngle[rad]"])
         gaia_scan_phi = gaia_scan_theta + np.pi / 2
         self.gaia_cos_phi = np.cos(gaia_scan_phi)
@@ -443,9 +444,9 @@ def _linear_pm_fit(
         """
         # Sovle y = A * x
         # construct A matrix
-        A_pmra = cos_phi * (epochs - epoch_ref_ra) / errs
+        A_pmra = cos_phi * (epochs.jyear - epoch_ref_ra) / errs
         A_raoff = cos_phi / errs
-        A_pmdec = sin_phi * (epochs - epoch_ref_dec) / errs
+        A_pmdec = sin_phi * (epochs.jyear - epoch_ref_dec) / errs
         A_decoff = sin_phi / errs
         A_matrix = np.vstack((A_raoff, A_decoff, A_pmra, A_pmdec)).T
 

orbitize/hipparcos.py

@@ -1,5 +1,4 @@
 import numpy as np
-from astropy.io import ascii
 import pandas as pd
 import emcee
 from scipy.stats import norm
@@ -17,29 +16,28 @@ class PMPlx_Motion(object):
     parallax and proper motion model (NO orbital motion is added in this class).
 
     Args:
-        times_mjd (np.array of float): times (in mjd) at which we have absolute astrometric
+        epochs_mjd (np.array of float): times (in mjd) at which we have absolute astrometric
             measurements
         alpha0 (float): measured RA position (in degrees) of the object at alphadec0_epoch (see below).
         delta0 (float): measured Dec position (in degrees) of the object at alphadec0_epoch (see below).
         alphadec0_epoch (float): a (fixed) reference time. For stars with Hipparcos data, this
-            should generally be 1991.25, but you can define it however you want. Absolute
+            should generally be J1991.25, but you can define it however you want. Absolute
             astrometric data (passed in via an orbitize! data table) should be defined
             as offsets from the reported position of the object at this epoch (with propagated
             uncertainties). For example, if you have two absolute astrometric measurements
-            taken with GRAVITY, as well as a Hipparcos-derived position (at epoch 1991.25),
-            alphadec0_epoch should be 1991.25, and you should pass in absolute astrometry
+            taken with GRAVITY, as well as a Hipparcos-derived position (at epoch J1991.25),
+            alphadec0_epoch should be J1991.25, and you should pass in absolute astrometry
             in terms of mas *offset* from the Hipparcos catalog position, with propagated
             errors of your measurement and the Hipparcos measurement.
     """
 
-    def __init__(self, epochs_mjd, alpha0, delta0, alphadec0_epoch=1991.25):
+    def __init__(self, epochs_mjd, alpha0, delta0, alphadec0_epoch=Time("J1991.25", format="jyear_str").jyear):
         self.epochs_mjd = epochs_mjd
         self.alphadec0_epoch = alphadec0_epoch
         self.alpha0 = alpha0
         self.delta0 = delta0
 
         epochs = Time(epochs_mjd, format="mjd")
-        self.epochs = epochs.decimalyear
 
         # compute Earth XYZ position in barycentric coordinates
         bary_pos, _ = get_body_barycentric_posvel("earth", epochs)
@@ -61,7 +59,7 @@ def compute_astrometric_model(self, samples, param_idx, epochs=None):
                 indices in an array of fitting parameters (generally
                 set to System.basis.param_idx).
             epochs: if None, use self.epochs for astrometric predictions. Otherwise,
-                use this array passed in [in decimalyear].
+                use this array passed in [in mjd].
 
         Returns:
             tuple of:
@@ -78,14 +76,14 @@ def compute_astrometric_model(self, samples, param_idx, epochs=None):
         delta_H0 = samples[param_idx["delta0"]]
 
         if epochs is None:
-            epochs = self.epochs
+            epochs = Time(self.epochs_mjd, format='mjd')
             X = self.X
             Y = self.Y
             Z = self.Z
         else:
             # compute Earth XYZ position in barycentric coordinates
             bary_pos, _ = get_body_barycentric_posvel(
-                "earth", Time(epochs, format="decimalyear")
+                "earth", Time(epochs, format="mjd")
             )
             X = bary_pos.x.value  # [au]
             Y = bary_pos.y.value  # [au]
@@ -105,7 +103,7 @@ def compute_astrometric_model(self, samples, param_idx, epochs=None):
                     X[i] * np.sin(np.radians(self.alpha0))
                     - Y[i] * np.cos(np.radians(self.alpha0))
                 )
-                + (epochs[i] - self.alphadec0_epoch) * pm_ra
+                + (epochs[i].jyear - self.alphadec0_epoch) * pm_ra
             )
             delta_C_array[i] = (
                 delta_H0
@@ -119,7 +117,7 @@ def compute_astrometric_model(self, samples, param_idx, epochs=None):
                     * np.sin(np.radians(self.delta0))
                     - Z[i] * np.cos(np.radians(self.delta0))
                 )
-                + (epochs[i] - self.alphadec0_epoch) * pm_dec
+                + (epochs[i].jyear - self.alphadec0_epoch) * pm_dec
             )
         return alpha_C_st_array, delta_C_array
 
@@ -137,9 +135,9 @@ class HipparcosLogProb(object):
     are described here for completeness. See Nielsen+ 2020 for more detail.
 
     - alpha0: RA offset from the reported Hipparcos position at a particular
-        epoch (usually 1991.25) [mas]
+        epoch (usually J1991.25) [mas]
     - delta0: Dec offset from the reported Hipparcos position at a particular
-        epoch (usually 1991.25) [mas]
+        epoch (usually J1991.25) [mas]
     - pm_ra: RA proper motion [mas/yr]
     - pm_dec: Dec proper motion [mas/yr]
     - plx: parallax [mas]
@@ -157,7 +155,7 @@ class HipparcosLogProb(object):
             zeros in the prefix if number is <100,000. (i.e. 27321 should be
             passed in as '027321').
         num_secondary_bodies (int): number of companions in the system
-        alphadec0_epoch (float): epoch (in decimal year) that the fitting
+        alphadec0_epoch (float): epoch (in Julian decimal year) that the fitting
             parameters alpha0 and delta0 are defined relative to (see above).
         renormalize_errors (bool): if True, normalize the scan errors to get
             chisq_red = 1, following Nielsen+ 2020 (eq 10). In general, this
@@ -173,7 +171,7 @@ def __init__(
         path_to_iad_file,
         hip_num,
         num_secondary_bodies,
-        alphadec0_epoch=1991.25,
+        alphadec0_epoch=Time("J1991.25", format="jyear_str").jyear,
         renormalize_errors=False,
     ):
         self.path_to_iad_file = path_to_iad_file
@@ -284,7 +282,11 @@ def __init__(
 
         n_lines = len(iad)
 
-        times = iad[1] + 1991.25
+
+        self.epochs = Time(
+            iad[1] * 365.25 + Time("J1991.25",format="jyear_str").jd, format='jd'
+        )
+
         self.cos_phi = iad[3]  # scan direction
         self.sin_phi = iad[4]
         self.R = iad[5]  # abscissa residual [mas]
@@ -295,18 +297,18 @@ def __init__(
 
         if n_lines - len(good_scans) > 0:
             print("{} Hipparcos scans rejected.".format(n_lines - len(good_scans)))
-        times = times[good_scans]
+        self.epochs = self.epochs[good_scans]
         self.cos_phi = self.cos_phi[good_scans]
         self.sin_phi = self.sin_phi[good_scans]
         self.R = self.R[good_scans]
         self.eps = self.eps[good_scans]
 
+        self.epochs_mjd = self.epochs.mjd
+
         # if the star has a type 1 (stochastic) solution, we need to undo the addition of a jitter term in quadrature
         self.eps = np.sqrt(self.eps**2 - self.var**2)
 
-        epochs = Time(times, format="decimalyear")
-        self.epochs = epochs.decimalyear
-        self.epochs_mjd = epochs.mjd
+
 
         self.hipparcos_plxpm_predictor = PMPlx_Motion(
             self.epochs_mjd,
@@ -316,7 +318,7 @@ def __init__(
         )
 
         if self.renormalize_errors:
-            D = len(epochs) - 6
+            D = len(self.epochs) - 6
             G = f2
 
             f = (G * np.sqrt(2 / (9 * D)) + 1 - (2 / (9 * D))) ** (3 / 2)
@@ -333,7 +335,7 @@ def __init__(
                 self.hipparcos_plxpm_predictor.X * np.sin(np.radians(self.alpha0))
                 - self.hipparcos_plxpm_predictor.Y * np.cos(np.radians(self.alpha0))
             )
-            + (self.epochs - 1991.25) * self.pm_ra0
+            + (self.epochs.jyear - self.alphadec0_epoch) * self.pm_ra0
         )
 
         changein_delta = (
@@ -347,7 +349,7 @@ def __init__(
                 * np.sin(np.radians(self.delta0))
                 - self.hipparcos_plxpm_predictor.Z * np.cos(np.radians(self.delta0))
             )
-            + (self.epochs - 1991.25) * self.pm_dec0
+            + (self.epochs.jyear - self.alphadec0_epoch) * self.pm_dec0
         )
 
         # compute abcissa point (Nielsen+ Eq 3)

orbitize/sampler.py

@@ -126,12 +126,9 @@ def _logl(self, params):
             )
 
         if self.system.gaia is not None:
-            gaiahip_epochs = Time(
-                np.append(
-                    self.system.gaia.hipparcos_epoch, self.system.gaia.gaia_epoch
-                ),
-                format="decimalyear",
-            ).mjd
+            gaiahip_epochs = np.append(
+                self.system.gaia.hipparcos_epoch.mjd, self.system.gaia.gaia_epoch.mjd
+            )
 
             # compute Ra/Dec predictions at the Gaia epoch
             raoff_model, deoff_model, _ = self.system.compute_all_orbits(

tests/test_gaia.py

@@ -156,11 +156,11 @@ def test_orbit_calculation():
     myGaia.ra = myHip.alpha0 + (
         myGaia.mas2deg
         * pm_a
-        * (myGaia.gaia_epoch - myGaia.hipparcos_epoch)
+        * (myGaia.gaia_epoch.jyear - myGaia.hipparcos_epoch.jyear)
         / np.cos(np.radians(myHip.delta0))
     )
     myGaia.dec = myHip.delta0 + (
-        myGaia.mas2deg * pm_d * (myGaia.gaia_epoch - myGaia.hipparcos_epoch)
+        myGaia.mas2deg * pm_d * (myGaia.gaia_epoch.jyear - myGaia.hipparcos_epoch.jyear)
     )
     test_samples = [sma, ecc, inc, aop, pan, tau, plx, m0, m1, a0, d0, pm_a, pm_d]
 
@@ -194,7 +194,7 @@ def test_orbit_calculation():
     myGaia.dec = myHip.delta0 + 1
 
     sma = 2 * (myGaia.dec - myHip.delta0) * deg2arcsec * (plx * mas2arcsec)  # [au]
-    per = 2 * (myGaia.gaia_epoch - myGaia.hipparcos_epoch)  # [yr]
+    per = 2 * (myGaia.gaia_epoch.jyear - myGaia.hipparcos_epoch.jyear)  # [yr]
     mtot = sma**3 / per**2
 
     test_samples[param_idx["sma1"]] = sma
@@ -203,7 +203,7 @@ def test_orbit_calculation():
 
     # passes through peri (+sma decl for e=0 orbits) at Hipparcos epoch
     # -> @ Gaia epoch, primary should be at +sma decl
-    tau = basis.tp_to_tau(myGaia.hipparcos_epoch, 58849, per)
+    tau = basis.tp_to_tau(myGaia.hipparcos_epoch.jyear, 58849, per)
     test_samples[param_idx["tau1"]] = tau
 
     # choose sma and mass so that Hipparcos/Gaia difference is only due to orbit.
@@ -308,8 +308,8 @@ def test_nointernet():
 
 if __name__ == "__main__":
     test_nointernet()
-    # test_dr2_edr3()
-    # test_system_setup()
-    # test_valueerror()
-    # test_orbit_calculation()
+    test_dr2_edr3()
+    test_system_setup()
+    test_valueerror()
+    test_orbit_calculation()
     test_hgca()