Version 1.3.0

.github/workflows/tests.yaml

@@ -24,7 +24,8 @@ jobs:
       - name: Set fortran compiler on macOS
         if: runner.os == 'macOS'
         run: |
-          echo "PATH=$(brew --prefix gfortran)/bin:$PATH" >> $GITHUB_ENV
+          brew install gfortran
+          export PATH="$(brew --prefix gfortran)/bin:$PATH" >> $GITHUB_ENV
 
       - name: Install package
         run: |

.gitignore

@@ -6,6 +6,7 @@
 doc
 src/pride/displacements/old.py
 dev
+validation
 
 # Compiled extension modules
 *.so

src/pride/astro/__init__.py

@@ -0,0 +1,17 @@
+from .lorentz_transformations import transform_position_from_gcrf_to_bcrf
+from .ephemerides import (
+    get_icrf_position_vector,
+    get_icrf_state_vector,
+    get_gcrf_position_vector,
+    get_gcrf_state_vector,
+    get_body_gravitational_parameter,
+)
+
+__all__ = [
+    "transform_position_from_gcrf_to_bcrf",
+    "get_icrf_state_vector",
+    "get_icrf_position_vector",
+    "get_gcrf_state_vector",
+    "get_gcrf_position_vector",
+    "get_body_gravitational_parameter",
+]

src/pride/astro/ephemerides.py

@@ -0,0 +1,83 @@
+import spiceypy as spice
+import numpy as np
+from ..logger import log
+
+
+def get_icrf_state_vector(target: str, epochs: np.ndarray) -> np.ndarray:
+    """Get cartesian state vector of a target in ICRF (SI units)
+
+    Uses spice.spkerz to get the state of a target in J2000, with respect to SSB, and with aberrations set to NONE. The output is transformed to SI units and casted to an (N, 6) numpy array.
+
+    :param target: Target name known by spice
+    :param epochs: Array of epochs in ET (TDB seconds past J2000)
+    :return state_vector: Aberrated, cartesian state vector of the target with respect to SSB in J2000 (ICRF) frame.
+    """
+
+    # Get state and LT from spice
+    cstate_tuple, _ = spice.spkezr(target, epochs, "J2000", "NONE", "SSB")
+
+    # Convert to meters and numpy array
+    return np.array(cstate_tuple, dtype=np.float64) * 1e3
+
+
+def get_icrf_position_vector(target: str, epochs: np.ndarray) -> np.ndarray:
+    """Get cartesian position vector of a target in ICRF (SI units)
+
+    Uses spice.spkpos to get the position of a target in J2000, with respect to SSB, and with aberrations set to NONE. The output is transformed to SI units and casted to an (N, 3) numpy array.
+
+    :param target: Target name known by spice
+    :param epochs: Array of epochs in ET (TDB seconds past J2000)
+    :return position_vector: Aberrated, cartesian position vector of the target with respect to SSB in J2000 (ICRF) frame.
+    """
+
+    # Get position and LT from spice
+    cpos_tuple, _ = spice.spkpos(target, epochs, "J2000", "NONE", "SSB")
+
+    # Convert to meters and numpy array
+    return np.array(cpos_tuple, dtype=np.float64) * 1e3
+
+
+def get_gcrf_position_vector(target: str, epochs: np.ndarray) -> np.ndarray:
+    """Get cartesian position vector of a target in GCRF (SI units)
+
+    Uses spice.spkpos to get the position of a target in J2000, with respect to EARTH, and with aberrations set to NONE. The output is transformed to SI units and casted to an (N, 3) numpy array.
+
+    :param target: Target name known by spice
+    :param epochs: Array of epochs in ET (TDB seconds past J2000)
+    :return position_vector: Aberrated, cartesian position vector of the target with respect to EARTH in J2000 (GCRF) frame.
+    """
+
+    # Get position and LT from spice
+    cpos_tuple, _ = spice.spkpos(target, epochs, "J2000", "NONE", "EARTH")
+
+    # Convert to meters and numpy array
+    return np.array(cpos_tuple, dtype=np.float64) * 1e3
+
+
+def get_gcrf_state_vector(target: str, epochs: np.ndarray) -> np.ndarray:
+    """Get cartesian state vector of a target in GCRF (SI units)
+
+    Uses spice.spkezr to get the position of a target in J2000, with respect to EARTH, and with aberrations set to NONE. The output is transformed to SI units and casted to an (N, 6) numpy array.
+
+    :param target: Target name known by spice
+    :param epochs: Array of epochs in ET (TDB seconds past J2000)
+    :return position_vector: Aberrated, cartesian state vector of the target with respect to EARTH in J2000 (GCRF) frame.
+    """
+
+    # Get position and LT from spice
+    cstate_tuple, _ = spice.spkezr(target, epochs, "J2000", "NONE", "EARTH")
+
+    # Convert to meters and numpy array
+    return np.array(cstate_tuple, dtype=np.float64) * 1e3
+
+
+def get_body_gravitational_parameter(target: str) -> float:
+    """Get gravitational parameter of target from SPICE kernels
+
+    Uses spice.bodvrd to retrieve the gravitational parameter of the target in the loaded version of the SPICE kernels. The result is transformed to SI units.
+
+    :param target: Target name known to SPICE
+    :return mu_target: Gravitational parameter of the target in SI units
+    """
+
+    return spice.bodvrd(target, "GM", 1)[1][0] * 1e9

src/pride/astro/lorentz_transformations.py

@@ -0,0 +1,40 @@
+import numpy as np
+from ..constants import L_C, CLIGHT, GAMMA_PPN
+from ..logger import log
+
+
+def transform_position_from_gcrf_to_bcrf(
+    x_target_gcrf: np.ndarray,
+    x_earth_bcrf: np.ndarray,
+    v_earth_bcrf: np.ndarray,
+    potential_at_geocenter: np.ndarray,
+) -> np.ndarray:
+    """Lorentz transformation of position vector from GCRF to BCRF
+
+    Source: 10.1051/0004-6361/201218885 (Duev 2012)
+
+    :param x_target_gcrf: GCRF position of the target as (N, 3) array
+    :param x_earth_bcrf: BCRF positon of the Earth as (N, 3) array
+    :param v_earth_bcrf: BCRF velocity of the Earth as (N, 3) array
+    :param potential_at_geocenter: Newtonian potential of all the solar system bodies, excluding the Earth, evaluated at the geocenter. Passed as (N,) array.
+    :return: BCRF position of the target as (N, 3) array
+    """
+
+    log.warning("Missing tests: transform_position_from_gcrf_to_bcrf")
+
+    # Precomputed quantities
+    clight2 = CLIGHT * CLIGHT
+
+    # Components of the expression
+    term1 = (1.0 - L_C - (GAMMA_PPN * potential_at_geocenter / clight2))[
+        :, None
+    ] * x_target_gcrf
+    term2 = (
+        0.5
+        * np.sum(v_earth_bcrf.T * x_target_gcrf.T, axis=0)[:, None]
+        * v_earth_bcrf
+        / clight2
+    )
+    term3 = x_earth_bcrf
+
+    return term1 - term2 + term3

src/pride/cli.py

@@ -1,5 +1,9 @@
 from .experiment import Experiment
 import argparse
+from . import coordinates as coord
+from .logger import log
+from astropy import time
+import numpy as np
 
 
 class PrideArgumentParser(argparse.ArgumentParser):
@@ -36,20 +40,79 @@ def process_experiment() -> None:
         for baseline in experiment.baselines:
 
             # Update baseline with data from observations
-            baseline.update_with_observations()
+            # Augment time stamps with +/- 1 second around each epoch
+            __augmented_tstamps: time.Time = (
+                baseline.tstamps[:, None]
+                + time.TimeDelta([-1, 0, 1], format="sec")[None, :]
+            ).ravel()  # type: ignore
+            a_tstamps = time.Time(
+                __augmented_tstamps,
+                scale="utc",
+                location=baseline.station.tectonic_corrected_location(
+                    __augmented_tstamps
+                ),
+            )
+            assert a_tstamps.location is not None
+
+            # Calculate geodetic coordinates of station
+            a_geodetic = a_tstamps.location.to_geodetic("GRS80")
+            a_lat = np.array(a_geodetic.lat.rad, dtype=float)
+            a_lon = np.array(a_geodetic.lon.rad, dtype=float)
+
+            # Calculate rotation matrices
+            augmented_eops = experiment.eops.at_epoch(a_tstamps, unit="arcsec")
+            # self.a_eops = eops.at_epoch(self.a_tstamps, unit="arcsec")
+            a_icrf2itrf = coord.icrf2itrf(augmented_eops, a_tstamps)
+            a_seu2itrf = coord.seu2itrf(a_lat, a_lon)
+
+            # Calculate derivative of ICRF to ITRF rotation matrix
+            dt_tdb: np.ndarray = (
+                (a_tstamps.tdb[2::3] - a_tstamps.tdb[0::3])
+                .to("s")  # type: ignore
+                .value
+            )
+            diff_icrf2itrf = a_icrf2itrf[2::3] - a_icrf2itrf[0::3]
+            dot_icrf2itrf = diff_icrf2itrf / dt_tdb[:, None, None]
+
+            # Get attributes at observation epochs
+            icrf2itrf = a_icrf2itrf[1::3]
+            seu2itrf = a_seu2itrf[1::3]
+
+            # Update observations with calculated data
+            for observation in baseline.observations:
+
+                # Get the index of the time stamps associated with the observation
+                flags = np.sum(
+                    observation.tstamps[:, None] == baseline.tstamps[None, :],
+                    axis=0,
+                    dtype=bool,
+                )
+
+                observation.icrf2itrf = icrf2itrf[flags]
+                observation.seu2itrf = seu2itrf[flags]
+                observation.dot_icrf2itrf = dot_icrf2itrf[flags]
 
             # Update station coordinates with geophysical displacements
-            baseline.update_station_with_geophysical_displacements()
+            baseline.update_station_with_geophysical_displacements(
+                experiment.displacement_models,
+                a_tstamps,
+                augmented_eops,
+                a_icrf2itrf,
+                a_seu2itrf,
+                a_lat,
+                a_lon,
+                icrf2itrf,
+                dot_icrf2itrf,
+            )
 
             for observation in baseline.observations:
 
                 # Calculate spherical coordinates of the source
                 observation.update_with_source_coordinates()
 
                 # Calculate delays for the observation
-                observation.calculate_delays()
+                observation.calculate_delays(experiment.delay_models)
 
-    # Save the output
     experiment.save_output()
 
     return None

src/pride/constants.py

@@ -1,6 +1,13 @@
 from astropy import time
+import spiceypy as spice
 
 J2000: time.Time = time.Time("2000-01-01T12:00:00", scale="tt").utc  # type: ignore
 """J2000 epoch in UTC"""
 L_C: float = 1.48082686741e-8
 """1 - d(TT)/d(TCG) [From IERS Conventions 2010]"""
+
+CLIGHT: float = spice.clight() * 1e3
+"""Speed of light in m/s [From SPICE]"""
+
+GAMMA_PPN: float = 1.0
+"""PPN parameter in general relativity"""

src/pride/coordinates/__init__.py

@@ -1,5 +1,23 @@
-"""Transformations between coordinate systems"""
+"""Transformations between coordinate systems
+
+About EOP bulletins
+--------------------
+
+According to the *Explanatory supplement to IERS bulletins A and B* [iers_bulletins]_, the IERS distributes two collections of EOPs: bulletins A and B. Bulletin A is published daily and is only meant to be used while bulletin B is not available. The IERS also distributes long-term collections of EOPs (e.g. C04), which are consistent with bulletin B.
+
+PRIDE accesses EOP information via Astropy, which includes both types of bulletins, and defaults to working with bulletin B unless the user specifies otherwise. As of May of 2025, Astropy's bulletin B corresponds to the C04 20 series, which is consistent with ITRF2020.
+
+.. warning:: As of version 1.2, the station coordinates in the internal catalogs are given in ITRF2000, while the default EOPs are consistent with ITRF2020. PRIDE is under development, and this issue will be fixed in the future.
+
+.. [iers_bulletins] https://hpiers.obspm.fr/iers/bul/bulb/explanatory.pdf
+
+"""
 
 from .core import EOP, seu2itrf, icrf2itrf, itrf2icrf
 
-__all__ = ["EOP", "seu2itrf", "icrf2itrf", "itrf2icrf"]
+__all__ = [
+    "EOP",
+    "seu2itrf",
+    "icrf2itrf",
+    "itrf2icrf",
+]

src/pride/coordinates/core.py

@@ -20,7 +20,17 @@ def __init__(
         bulletin: Literal["A", "B"],
         initial_epoch: time.Time,
         final_epoch: time.Time,
+        margin_in_days: int = 8,
     ) -> None:
+        """Constructor for EOP interface class
+
+        .. Note:: You can find more information about EOP bulletins at the API documentation of this module
+
+        :param bulletin: EOP bulletin to use [A or B]
+        :param initial_epoch: Epoch from which to load EOPs
+        :param final_epoch: Epoch until which to load EOPs
+        :param margin_in_days: Number of days by which to augment the time range defined by the initial and final epochs.
+        """
 
         # Select EOP bulletin
         match bulletin:
@@ -34,32 +44,46 @@ def __init__(
                 )
                 exit(1)
 
-        # Calculate range of MJDs
-        initial_mjd = (initial_epoch.mjd // 1) - 8  # type: ignore
-        final_mjd = (final_epoch.mjd // 1) + 8  # type: ignore
-        idx = (
+        # Augment initial and final epochs with margin
+        initial_mjd = (initial_epoch.mjd // 1) - margin_in_days  # type: ignore
+        final_mjd = (final_epoch.mjd // 1) + margin_in_days  # type: ignore
+        self.__mjd_range = (initial_mjd, final_mjd)
+
+        # Find entries of EOP table within the time range
+        eop_table_mask: np.ndarray = (
             np.argwhere(
-                (eop_table["MJD"].data >= initial_mjd)
-                & (eop_table["MJD"].data <= final_mjd)
+                (eop_table["MJD"].data >= initial_mjd)  # type: ignore
+                & (eop_table["MJD"].data <= final_mjd)  # type: ignore
             )
             .ravel()
             .astype(int)
         )
-        mjd_list = eop_table["MJD"][idx].value
-        self.__mjd_range = (initial_mjd, final_mjd)
+        assert isinstance(eop_table_mask, np.ndarray)
+
+        # Get MJD for table entries in the time range
+        relevant_mjds = eop_table["MJD"][eop_table_mask].data  # type: ignore
 
         # EOPs at integer epochs
         self.__eops_dict = {
-            "xp": interpolate.interp1d(mjd_list, eop_table["PM_x"][idx].value),
-            "yp": interpolate.interp1d(mjd_list, eop_table["PM_y"][idx].value),
+            "xp": interpolate.interp1d(
+                relevant_mjds,
+                eop_table["PM_x"][eop_table_mask].value,  # type: ignore
+            ),
+            "yp": interpolate.interp1d(
+                relevant_mjds,
+                eop_table["PM_y"][eop_table_mask].value,  # type: ignore
+            ),
             "ut1_utc": interpolate.interp1d(
-                mjd_list, eop_table["UT1_UTC"][idx].value
+                relevant_mjds,
+                eop_table["UT1_UTC"][eop_table_mask].value,  # type: ignore
             ),
             "dx": interpolate.interp1d(
-                mjd_list, eop_table["dX_2000A"][idx].value
+                relevant_mjds,
+                eop_table["dX_2000A"][eop_table_mask].value,  # type: ignore
             ),
             "dy": interpolate.interp1d(
-                mjd_list, eop_table["dY_2000A"][idx].value
+                relevant_mjds,
+                eop_table["dY_2000A"][eop_table_mask].value,  # type: ignore
             ),
         }
 
@@ -85,6 +109,10 @@ def __eops(self, mjds) -> np.ndarray:
 
     @staticmethod
     def from_experiment(experiment: "Experiment") -> "EOP":
+
+        raise DeprecationWarning(
+            "Initializing EOPs from experiment is deprecated"
+        )
         return EOP(
             experiment.setup.internal["eop_bulletin"],
             experiment.initial_epoch,

src/pride/delays/__init__.py

@@ -7,10 +7,10 @@
     from .core import Delay
 
 DELAY_MODELS: dict[str, type["Delay"]] = {
-    Geometric.name: Geometric,
-    Tropospheric.name: Tropospheric,
-    Ionospheric.name: Ionospheric,
-    AntennaDelays.name: AntennaDelays,
+    Geometric.__name__: Geometric,
+    Tropospheric.__name__: Tropospheric,
+    Ionospheric.__name__: Ionospheric,
+    AntennaDelays.__name__: AntennaDelays,
 }
 
 __all__ = [