From 71a02de1931bd61b42a5f943e9ee5d16b1a55b02 Mon Sep 17 00:00:00 2001 From: Caleb Morrison Date: Mon, 15 Sep 2025 17:23:44 -0700 Subject: [PATCH 01/17] Add R3 derivative function --- Tests/test_matrix_derivatives.py | 25 ++++++++++++++++ .../Utilities/TransformationMatrices.py | 29 +++++++++++++++++++ 2 files changed, 54 insertions(+) create mode 100644 Tests/test_matrix_derivatives.py diff --git a/Tests/test_matrix_derivatives.py b/Tests/test_matrix_derivatives.py new file mode 100644 index 0000000..2ac5da0 --- /dev/null +++ b/Tests/test_matrix_derivatives.py @@ -0,0 +1,25 @@ +import numpy as np + +from TerraFrame.Utilities import TransformationMatrices + + +def test_r3_derivative(): + n = 25 + t = np.linspace(0, 24 * 60 * 60, n) + psi = [x ** 1.2 for x in t] + dpsidt = 1.2 + dt = 1e-4 + dpsi = dpsidt * dt + + error = np.zeros((n,)) + + for i in reversed(range(n)): + r1 = TransformationMatrices.r3(psi[i]) + r2 = TransformationMatrices.r3(psi[i] + dpsi) + + drdt_e = (r2 - r1) / dt + drdt = TransformationMatrices.dr3dt(psi[i], dpsidt) + + error[i] = np.max(np.abs(drdt_e - drdt)) + + assert (max(error) < 1e-4) diff --git a/src/TerraFrame/Utilities/TransformationMatrices.py b/src/TerraFrame/Utilities/TransformationMatrices.py index 470189f..27c0aa3 100644 --- a/src/TerraFrame/Utilities/TransformationMatrices.py +++ b/src/TerraFrame/Utilities/TransformationMatrices.py @@ -90,6 +90,35 @@ def r3(psi): return r +def dr3dt(psi, dpsi_dt): + """ + This function computes the time derivative of the R3 rotation matrix. As + per Kaplan (2005), R3 is defined as: + + [A] rotation matrix to transform column 3-vectors from one cartesian + coordinate system to another. Final system is formed by rotating original + system about its own z-axis by angle φ (counterclockwise as viewed from + the +z direction): + + Source: + Kaplan, G. H., 2005, U.S. Naval Observatory Circular No. 179 (Washington: + USNO), page xi + + :param psi: Input rotation angle in radians + :param dpsi_dt: The time derivative of psi where dt is in seconds + :return: dR3dt matrix + :type psi: float + :type dpsi_dt: float + :rtype: np.ndarray + """ + + dr_dt = np.array([[-np.sin(psi) * dpsi_dt, np.cos(psi) * dpsi_dt, 0.0], + [-np.cos(psi) * dpsi_dt, -np.sin(psi) * dpsi_dt, 0.0], + [0.0, 0.0, 0.0]]) + + return dr_dt + + def euler_angles_from_transformation(t_m): """ This function takes a transformation matrix and calculates the corresponding From ec89a6a45b98417b772298329a0a50fceadc4515 Mon Sep 17 00:00:00 2001 From: Caleb Morrison Date: Sun, 21 Sep 2025 17:45:20 -0700 Subject: [PATCH 02/17] Add some derivative calculation functions --- src/TerraFrame/CelestialTerrestrial.py | 106 +++++++++ src/TerraFrame/Earth.py | 17 +- .../PrecessionNutation/SeriesExpansion.py | 9 +- src/TerraFrame/Utilities/Conversions.py | 14 ++ src/TerraFrame/Utilities/Interpolation.py | 9 +- .../Utilities/TransformationMatrices.py | 211 +++++++++++++++++- 6 files changed, 355 insertions(+), 11 deletions(-) diff --git a/src/TerraFrame/CelestialTerrestrial.py b/src/TerraFrame/CelestialTerrestrial.py index 803e582..0047c86 100644 --- a/src/TerraFrame/CelestialTerrestrial.py +++ b/src/TerraFrame/CelestialTerrestrial.py @@ -4,6 +4,8 @@ import datetime +from scipy.differentiate import derivative + import TerraFrame.Utilities.Conversions from TerraFrame.PrecessionNutation import SeriesExpansion from TerraFrame.Utilities import (Conversions, Time, BulletinData, @@ -116,3 +118,107 @@ def gcrs_to_itrs(self, time): t_gi = self.itrs_to_gcrs(time) return t_gi.T + + + def itrs_to_gcrs_angular_vel(self, time): + if isinstance(time, datetime.datetime): + time = Time.JulianDate.julian_date_from_pydatetime(time) + + assert isinstance(time, JulianDate) + + jd_tt = Conversions.any_to_tt(time) + + if time.time_scale == Time.TimeScales.UTC: + jd_utc = time + else: + jd_utc = Conversions.tt_to_utc(jd_tt) + + # We also need time in Modified Julian Date (MJD) for the Bulletin + # corrections lookup table. + mjd_utc = Time.JulianDate.julian_date_to_modified_julian_date(jd_utc) + + # We also need time in UT1 for the ERA + jd_ut1 = Conversions.tt_to_ut1(jd_tt) + + # Time needs to be in Julian centuries + jdc_tt = Time.JulianDate.julian_day_datetime_to_century_datetime(jd_tt) + + # For the given terrestrial time (TT), call the routines to obtain the + # IAU 2006/2000A X and Y from series. Then calculate "s" which is the + # CIO locator + (cip_x, dcip_x_dt) = self.se_cip_x.compute(jdc_tt, derivative=True) + (cip_y, dcip_y_dt) = self.se_cip_y.compute(jdc_tt, derivative=True) + (sxy2, dsxy2_dt) = self.se_cip_sxy2.compute(jdc_tt, derivative=True) + cip_s = sxy2 - cip_x * cip_y / 2.0 + + dcip_s_dt = dsxy2_dt - 0.5 * cip_y * dcip_x_dt - 0.5 * cip_x * dcip_y_dt + + # Any CIP corrections ∆X, ∆Y can now be applied, and the corrected + # X, Y, and s can be used to construct the Celestial Intermediate + # Reference System (CIRS) to Geocentric Celestial Reference System + # (GCRS) matrix: CIRS -> GCRS. + # Get corrections by interpolating in the IERS Bulletin A data + if self._user_nutation_corrections: + dx = self.bd.f_nc_dx(float(mjd_utc)) + dy = self.bd.f_nc_dy(float(mjd_utc)) + + cip_x += Conversions.mas_to_rad(dx) + cip_y += Conversions.mas_to_rad(dy) + + ddx_dt = self.bd.f_nc_dx(float(mjd_utc), derivative=True) + ddy_dt = self.bd.f_nc_dy(float(mjd_utc), derivative=True) + + dcip_x_dt += ddx_dt + dcip_y_dt += ddy_dt + + # Create the first transformation matrix + t_gc = TransformationMatrices.cirs_to_gcrs(cip_x, cip_y, cip_s) + dt_gc_dt = (TransformationMatrices + .cirs_to_gcrs_derivative(cip_x, cip_y, cip_s, + dcip_x_dt, dcip_y_dt, dcip_s_dt)) + + # The Earth rotation matrix is the transformation from the Terrestrial + # Intermediate Reference System (TIRS) to the Celestial Intermediate + # Reference System (CIRS): TIRS -> CIRS. + # This function uses normal JD time in UT1. + t_ct = TransformationMatrices.earth_rotation_matrix(jd_ut1) + + dt_ct_dt = (TransformationMatrices + .earth_rotation_matrix_derivative(jd_ut1)) + + # Given polar motion offsets pm_x and pm_y, along with the Terrestrial + # Intermediate Origin (TIO) locator (s prime or sp), the International + # Terrestrial Reference System (ITRS) to Terrestrial Intermediate + # Reference System (TIRS) transformation matrix can be constructed: + # ITRS -> TIRS. + if self._user_polar_motion: + pm_x = self.bd.f_pm_x(float(mjd_utc)) + pm_y = self.bd.f_pm_y(float(mjd_utc)) + dpm_x_dt = self.bd.f_pm_x(float(mjd_utc), derivative=True) + dpm_y_dt = self.bd.f_pm_y(float(mjd_utc), derivative=True) + else: + pm_x = 0.0 + pm_y = 0.0 + dpm_x_dt = 0.0 + dpm_y_dt = 0.0 + + sp = TransformationMatrices.calculate_s_prime(jdc_tt) + dsp_dt = TransformationMatrices.calculate_s_prime_derivative() + + pm_x = TerraFrame.Utilities.Conversions.arcsec_to_rad(pm_x) + pm_y = TerraFrame.Utilities.Conversions.arcsec_to_rad(pm_y) + + t_ti = TransformationMatrices.itrs_to_tirs(pm_x, pm_y, sp) + dt_ti_dt = TransformationMatrices.itrs_to_tirs_derivative( + pm_x, pm_y, sp, + dpm_x_dt, dpm_y_dt, dsp_dt) + + # Construct the final transformation matrix: ITRS -> GCRS + t_gi = t_gc @ t_ct @ t_ti + + self.t_gi = t_gi + self.t_gc = t_gc + self.t_ct = t_ct + self.t_ti = t_ti + + return t_gi diff --git a/src/TerraFrame/Earth.py b/src/TerraFrame/Earth.py index f10e7eb..7181d01 100644 --- a/src/TerraFrame/Earth.py +++ b/src/TerraFrame/Earth.py @@ -5,7 +5,7 @@ import math from abc import ABC, abstractmethod -from TerraFrame.Utilities import Time +from TerraFrame.Utilities import Time, Conversions from TerraFrame.Utilities.Time import JulianDate @@ -308,3 +308,18 @@ def earth_rotation_angle(time): era = math.fmod(era, 2.0 * math.pi) return era + + +def earth_rotation_angle_derivative(): + """ + This function computes the time derivative of the earth rotation angle at a + given datetime in UT1. + + :return: Earth rotation angle in radians/sec + :rtype: float + """ + + dera_dt = 2.0 * math.pi * 1.00273781191135448 + dera_dt = Conversions.seconds_to_days(dera_dt) + + return dera_dt diff --git a/src/TerraFrame/PrecessionNutation/SeriesExpansion.py b/src/TerraFrame/PrecessionNutation/SeriesExpansion.py index 7b0a62d..299adb6 100644 --- a/src/TerraFrame/PrecessionNutation/SeriesExpansion.py +++ b/src/TerraFrame/PrecessionNutation/SeriesExpansion.py @@ -52,15 +52,22 @@ def __init__(self, data_file_path, polynomial_coefficients): super().__init__(data_file_path) self._polynomial_coefficients = polynomial_coefficients - def compute(self, t): + def compute(self, t, derivative=False): t = float(t) # units are micro-arcseconds poly_part = 0.0 + poly_part_dt = 0.0 for j in range(len(self._polynomial_coefficients)): poly_part += self._polynomial_coefficients[j] * t ** j + if derivative: + for j in range(1, len(self._polynomial_coefficients)): + poly_part_dt += self._polynomial_coefficients[j] * t ** (j - 1) + + # TODO: Finish adding derivative stuff below + # Initialize all the argument parameters. "argument" is the term that # IERS uses to refer to the input to the trigonometric functions. The # order is tightly coupled with the file format. diff --git a/src/TerraFrame/Utilities/Conversions.py b/src/TerraFrame/Utilities/Conversions.py index a811f2a..65bd0f1 100644 --- a/src/TerraFrame/Utilities/Conversions.py +++ b/src/TerraFrame/Utilities/Conversions.py @@ -295,3 +295,17 @@ def seconds_to_days(value): """ return value / 86400 + + +def days_to_centuries(value): + """ + This function takes a numeric input representing days and scales it to + centuries. The convertion assumes there is exactly 365 days in a year. + + :param value: Input value in days + :type value: int | float + :return: Output value in centuries + :rtype: int | float + """ + + return value / (365 * 100) diff --git a/src/TerraFrame/Utilities/Interpolation.py b/src/TerraFrame/Utilities/Interpolation.py index 822f988..2d21a3e 100644 --- a/src/TerraFrame/Utilities/Interpolation.py +++ b/src/TerraFrame/Utilities/Interpolation.py @@ -40,7 +40,7 @@ def _get_index(self, xv): return index - def __call__(self, xv): + def __call__(self, xv, derivative=False): xv = ensure_iterable(xv) yv = len(xv) * [0.0] @@ -62,7 +62,12 @@ def __call__(self, xv): x2 = self._x[index] x1 = self._x[index - 1] - yv[i] = (y2 - y1) / (x2 - x1) * (xv[i] - x1) + y1 + m = (y2 - y1) / (x2 - x1) + + if not derivative: + yv[i] = m * (xv[i] - x1) + y1 + else: + yv[i] = m if len(yv) == 1: return yv[0] diff --git a/src/TerraFrame/Utilities/TransformationMatrices.py b/src/TerraFrame/Utilities/TransformationMatrices.py index 27c0aa3..5ad5def 100644 --- a/src/TerraFrame/Utilities/TransformationMatrices.py +++ b/src/TerraFrame/Utilities/TransformationMatrices.py @@ -37,6 +37,37 @@ def r1(phi): return r +def dr1dt(phi, dphi_dt): + """ + This function computes the time derivative of the R2 rotation matrix. As + per Kaplan (2005), R2 is defined as: + + [A] rotation matrix to transform column 3-vectors from one cartesian + coordinate system to another. Final system is formed by rotating original + system about its own y-axis by angle φ (counterclockwise as viewed from + the +y direction): + + Source: + Kaplan, G. H., 2005, U.S. Naval Observatory Circular No. 179 (Washington: + USNO), page xi + + :param phi: Input rotation angle in radians + :param dphi_dt: The time derivative of phi where dt is in seconds + :return: dR3dt matrix + :type phi: float + :type dphi_dt: float + :rtype: np.ndarray + """ + + dr_dt = np.array( + [[0.0, 0.0, 0.0], + [0.0, -np.sin(phi) * dphi_dt, np.cos(phi) * dphi_dt], + [0.0, -np.cos(phi) * dphi_dt, -np.sin(phi) * dphi_dt]] + ) + + return dr_dt + + def r2(theta): """ This function computes the R2 rotation matrix. As per Kaplan (2005), R2 is @@ -57,12 +88,46 @@ def r2(theta): :rtype: np.ndarray """ - r = np.array([[np.cos(theta), 0.0, -np.sin(theta)], [0.0, 1.0, 0.0], - [np.sin(theta), 0.0, np.cos(theta)]]) + r = np.array( + [[np.cos(theta), 0.0, -np.sin(theta)], + [0.0, 1.0, 0.0], + [np.sin(theta), 0.0, np.cos(theta)]] + ) return r +def dr2dt(theta, dtheta_dt): + """ + This function computes the time derivative of the R2 rotation matrix. As + per Kaplan (2005), R2 is defined as: + + [A] rotation matrix to transform column 3-vectors from one cartesian + coordinate system to another. Final system is formed by rotating original + system about its own y-axis by angle φ (counterclockwise as viewed from + the +y direction): + + Source: + Kaplan, G. H., 2005, U.S. Naval Observatory Circular No. 179 (Washington: + USNO), page xi + + :param theta: Input rotation angle in radians + :param dtheta_dt: The time derivative of theta where dt is in seconds + :return: dR3dt matrix + :type theta: float + :type dtheta_dt: float + :rtype: np.ndarray + """ + + dr_dt = np.array( + [[-np.sin(theta) * dtheta_dt, 0.0, -np.cos(theta) * dtheta_dt], + [0.0, 0.0, 0.0], + [np.cos(theta) * dtheta_dt, 0.0, -np.sin(theta) * dtheta_dt]] + ) + + return dr_dt + + def r3(psi): """ This function computes the R3 rotation matrix. As per Kaplan (2005), R3 is @@ -84,8 +149,10 @@ def r3(psi): """ r = np.array( - [[np.cos(psi), np.sin(psi), 0.0], [-np.sin(psi), np.cos(psi), 0.0], - [0.0, 0.0, 1.0]]) + [[np.cos(psi), np.sin(psi), 0.0], + [-np.sin(psi), np.cos(psi), 0.0], + [0.0, 0.0, 1.0]] + ) return r @@ -112,9 +179,11 @@ def dr3dt(psi, dpsi_dt): :rtype: np.ndarray """ - dr_dt = np.array([[-np.sin(psi) * dpsi_dt, np.cos(psi) * dpsi_dt, 0.0], - [-np.cos(psi) * dpsi_dt, -np.sin(psi) * dpsi_dt, 0.0], - [0.0, 0.0, 0.0]]) + dr_dt = np.array( + [[-np.sin(psi) * dpsi_dt, np.cos(psi) * dpsi_dt, 0.0], + [-np.cos(psi) * dpsi_dt, -np.sin(psi) * dpsi_dt, 0.0], + [0.0, 0.0, 0.0]] + ) return dr_dt @@ -226,6 +295,31 @@ def calculate_s_prime(time): return s_prime +def calculate_s_prime_derivative(): + """ + This function computes the time derivative of the Terrestrial Intermediate + Origin (TIO) locator called s' (or s prime) per IERS Conventions (2010). + + Note that technically, per IERS Conventions (2010), the input time should + be Barycentric Dynamical Time (TDB) but the difference between TDB and TT is + already small. Additionally, we do not consider effects outside the + Geocentric Celestial Reference System (GCRS) and the primary driver of the + TDB vs TT difference is earth's mean anomaly in its orbit. The error from + this simplication is less than a microarcsecond in nutation. + + :return: s prime time derivative + :rtype: float + """ + + # This is an approximation good for the next century. See section 5.5.2 of + # IERS Conventions (2010) for more context. + dt_s_prime_dt = Conversions.days_to_centuries(-47e-6) + dt_s_prime_dt = Conversions.seconds_to_days(dt_s_prime_dt) + dt_s_prime_dt = Conversions.arcsec_to_rad(dt_s_prime_dt) + + return dt_s_prime_dt + + def cirs_to_gcrs(x, y, s): """ This function computes the transformation matrix from the Celestial @@ -258,6 +352,52 @@ def cirs_to_gcrs(x, y, s): return t_gc +def cirs_to_gcrs_derivative(x, y, s, dx_dt, dy_dt, ds_dt): + """ + This function computes the time derivative of the transformation matrix + from the Celestial Intermediate Reference System (CIRS) to the Geocentric + Celestial Reference System (GCRS) per IERS Conventions (2010). + + x and y are coordinates of the Celestial Intermediate Pole (CIP) and s is + the Celestial Intermediate Origin (CIO) locator parameter which provides + the position of the CIO on the equator of the CIP. + + :param x: X coordinate of the CIP + :type x: float + :param y: Y coordinate of the CIP + :type y: float + :param s: CIO location parameter + :type s: float + :param dx_dt: X coordinate of the CIP time derivative + :type dx_dt: float + :param dy_dt: Y coordinate of the CIP time derivative + :type dy_dt: float + :param ds_dt: CIO location parameter time derivative + :type ds_dt: float + :return: CGRS to CIRS transformation matrix time derivative + :rtype: np.ndarray + """ + + # This should never be true in reality + assert (1.0 - x ** 2 - y ** 2 > 0.0) + + # e and d formulas from Capitaine (2003) + e = np.atan2(y, x) + + de_dt = (x * dy_dt - y * dx_dt)/(x**2 + y**2) + + d = np.atan2(np.sqrt(x ** 2 + y ** 2), np.sqrt(1 - x ** 2 - y ** 2)) + + dd_dt = (np.sqrt(-1 - 1/(x**2 + y**2 - 1)) * (x * dx_dt + y * dy_dt) / + (x**2 + y**2)) + + dt_gc_dt = (r3(-e) @ r2(-d) @ dr3dt(s, ds_dt) + + r3(-e) @ dr2dt(-d, -dd_dt) @ r3(s) + + dr3dt(-e, -de_dt) @ r2(-d) @ r3(s)) + + return dt_gc_dt + + def earth_rotation_matrix(time): """ This function computes the earth rotation matrix at a given datetime in UT1. @@ -277,6 +417,27 @@ def earth_rotation_matrix(time): return r_era +def earth_rotation_matrix_derivative(time): + """ + This function computes the time derivative of the earth rotation matrix + at a given datetime in UT1. + + :param time: JulianDate in UT1 + :return: TIRS to ITRS transformation matrix + :type time: JulianDate + :rtype: np.ndarray + """ + + assert (time.time_scale == Time.TimeScales.UT1) + + era = TerraFrame.Earth.earth_rotation_angle(time) + dera_dt = TerraFrame.Earth.earth_rotation_angle_derivative() + + dr_era_dt = dr3dt(-era, -dera_dt) + + return dr_era_dt + + def itrs_to_tirs(pm_x, pm_y, sp): """ This function computes the transformation matrix from the International @@ -300,3 +461,39 @@ def itrs_to_tirs(pm_x, pm_y, sp): t_ti = (r3(-sp) @ r2(pm_x) @ r1(pm_y)) return t_ti + + +def itrs_to_tirs_derivative(pm_x, pm_y, sp, dpm_x_dt, dpm_y_dt, dsp_dt): + """ + This function computes the transformation matrix from the International + Terrestrial Reference System (ITRS) to the Terrestrial Intermediate + Reference System (TIRS) per IERS Conventions (2010). + + pm_x and pm_y are coordinates of polar motion and sp (s') is the + Terrestrial Intermediate Origin (TIO) locator parameter which provides + the position of the TIO on the equator of the CIP. + + :type pm_x: float + :param pm_x: Polar motion x coordinate + :type pm_y: float + :param pm_y: Polar motion y coordinate + :type sp: float + :param sp: TIO location parameter + :type dpm_x_dt: float + :param pm_x: Polar motion x coordinate time derivative + :type dpm_y_dt: float + :param pm_y: Polar motion y coordinate time derivative + :type dsp_dt: float + :param sp: TIO location parameter time derivative + :return: TIRS to ITRS transformation matrix time derivative + :rtype: np.ndarray + """ + + dr1_dt = dr1dt(pm_y, dpm_y_dt) + dr2_dt = dr2dt(pm_x, dpm_x_dt) + dr3_dt = dr2dt(-sp, -dsp_dt) + + dt_ti_dt = (r3(-sp) @ r2(pm_x) @ dr1_dt + r3(-sp) @ dr2_dt @ r1(pm_y) + + dr3_dt @ r2(pm_x) @ r1(pm_y)) + + return dt_ti_dt From b1c001b3228fbc03ff483c121a07415576754c82 Mon Sep 17 00:00:00 2001 From: Caleb Morrison Date: Sun, 21 Sep 2025 17:48:16 -0700 Subject: [PATCH 03/17] Remove scipy import that got automatically added --- src/TerraFrame/CelestialTerrestrial.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/TerraFrame/CelestialTerrestrial.py b/src/TerraFrame/CelestialTerrestrial.py index 0047c86..e3e6693 100644 --- a/src/TerraFrame/CelestialTerrestrial.py +++ b/src/TerraFrame/CelestialTerrestrial.py @@ -4,8 +4,6 @@ import datetime -from scipy.differentiate import derivative - import TerraFrame.Utilities.Conversions from TerraFrame.PrecessionNutation import SeriesExpansion from TerraFrame.Utilities import (Conversions, Time, BulletinData, From b429591cf59c9c7040ca25de3a5335dbf65559cb Mon Sep 17 00:00:00 2001 From: Caleb Morrison Date: Mon, 22 Sep 2025 21:15:55 -0700 Subject: [PATCH 04/17] Add mean arguments time derivatives --- .../PrecessionNutation/Arguments.py | 405 +++++++++++++++++- .../PrecessionNutation/SeriesExpansion.py | 35 +- src/TerraFrame/Utilities/Conversions.py | 16 + 3 files changed, 440 insertions(+), 16 deletions(-) diff --git a/src/TerraFrame/PrecessionNutation/Arguments.py b/src/TerraFrame/PrecessionNutation/Arguments.py index 570c603..90e7767 100644 --- a/src/TerraFrame/PrecessionNutation/Arguments.py +++ b/src/TerraFrame/PrecessionNutation/Arguments.py @@ -21,7 +21,8 @@ def mean_anomaly_of_the_moon(time): :type time: float :param time: Terrestrial time measured in Julian centuries. - :return: + :return: Mean anomaly of the Moon + :rtype: float """ # The polynomial is in arcseconds @@ -37,6 +38,35 @@ def mean_anomaly_of_the_moon(time): return value +def mean_anomaly_of_the_moon_derivative(time): + """ + This function computes the mean anomaly of the Moon time derivative per + IERS Conventions (2010) + + Note that technically, per IERS Conventions (2010), the input time should + be Barycentric Dynamical Time (TDB) but the difference between TDB and TT is + already small. Additionally, we do not consider effects outside the + Geocentric Celestial Reference System (GCRS) and the primary driver of the + TDB vs TT difference is earth's mean anomaly in its orbit. The error from + this simplication is less than a microarcsecond in nutation. + + :type time: float + :param time: Terrestrial time measured in Julian centuries. + :return: Mean anomaly of the Moon time derivative + :rtype: float + """ + + # The polynomial is in arcseconds + value = (1717915923.2178 + 63.7584 * time + 0.15490500000000001 * + time ** 2 - 0.0009788 * time ** 3) + + value = Conversions.arcsec_to_rad(value) + + value = Conversions.seconds_to_centuries(value) + + return value + + def mean_anomaly_of_the_sun(time): """ This function computes the mean anomaly of the Sun per IERS Conventions @@ -51,7 +81,8 @@ def mean_anomaly_of_the_sun(time): :type time: float :param time: Terrestrial time measured in Julian centuries. - :return: + :return: Mean anomaly of the Sun + :rtype: float """ # The polynomial is in arcseconds @@ -66,6 +97,35 @@ def mean_anomaly_of_the_sun(time): return value +def mean_anomaly_of_the_sun_derivative(time): + """ + This function computes the mean anomaly of the Sun time derivative per + IERS Conventions (2010) + + Note that technically, per IERS Conventions (2010), the input time should + be Barycentric Dynamical Time (TDB) but the difference between TDB and TT is + already small. Additionally, we do not consider effects outside the + Geocentric Celestial Reference System (GCRS) and the primary driver of the + TDB vs TT difference is earth's mean anomaly in its orbit. The error from + this simplication is less than a microarcsecond in nutation. + + :type time: float + :param time: Terrestrial time measured in Julian centuries. + :return: Mean anomaly of the Sun time derivative + :rtype: float + """ + + # The polynomial is in arcseconds + value = (129596581.0481 - 1.1064 * time + 0.000408 * time ** 2 - + 0.00004596 * time ** 3) + + value = Conversions.arcsec_to_rad(value) + + value = Conversions.seconds_to_centuries(value) + + return value + + def mean_longitude_moon_minus_ascending_node(time): """ This function computes Mean longitude of the Moon minus the ascending node @@ -80,7 +140,8 @@ def mean_longitude_moon_minus_ascending_node(time): :type time: float :param time: Terrestrial time measured in Julian centuries. - :return: + :return: Mean longitude of the Moon minus the ascending node + :rtype: float """ # The polynomial is in arcseconds @@ -95,6 +156,35 @@ def mean_longitude_moon_minus_ascending_node(time): return value +def mean_longitude_moon_minus_ascending_node_derivative(time): + """ + This function computes Mean longitude of the Moon minus the ascending node + time derivative per IERS Conventions (2010). + + Note that technically, per IERS Conventions (2010), the input time should + be Barycentric Dynamical Time (TDB) but the difference between TDB and TT is + already small. Additionally, we do not consider effects outside the + Geocentric Celestial Reference System (GCRS) and the primary driver of the + TDB vs TT difference is earth's mean anomaly in its orbit. The error from + this simplication is less than a microarcsecond in nutation. + + :type time: float + :param time: Terrestrial time measured in Julian centuries. + :return: Mean longitude of the Moon minus the ascending node time derivative + :rtype: float + """ + + # The polynomial is in arcseconds + value = (1739527262.8478-25.5024 * time - 0.003111 * time ** 2 + + 0.00001668 * time ** 3) + + value = Conversions.arcsec_to_rad(value) + + value = Conversions.seconds_to_centuries(value) + + return value + + def mean_elongation_of_the_moon_from_the_sun(time): """ This function computes Mean Elongation of the Moon from the Sun per IERS @@ -109,7 +199,8 @@ def mean_elongation_of_the_moon_from_the_sun(time): :type time: float :param time: Terrestrial time measured in Julian centuries. - :return: + :return: Mean Elongation of the Moon from the Sun + :rtype: float """ # The polynomial is in arcseconds @@ -124,6 +215,35 @@ def mean_elongation_of_the_moon_from_the_sun(time): return value +def mean_elongation_of_the_moon_from_the_sun_derivative(time): + """ + This function computes Mean Elongation of the Moon from the Sun time + derivative per IERS Conventions (2010). + + Note that technically, per IERS Conventions (2010), the input time should + be Barycentric Dynamical Time (TDB) but the difference between TDB and TT is + already small. Additionally, we do not consider effects outside the + Geocentric Celestial Reference System (GCRS) and the primary driver of the + TDB vs TT difference is earth's mean anomaly in its orbit. The error from + this simplication is less than a microarcsecond in nutation. + + :type time: float + :param time: Terrestrial time measured in Julian centuries. + :return: Mean Elongation of the Moon from the Sun time derivative + :rtype: float + """ + + # The polynomial is in arcseconds + value = (1602961601.209 - 12.7412 * time + 0.019779 * time ** 2 - + 0.00012676 * time ** 3) + + value = Conversions.arcsec_to_rad(value) + + value = Conversions.seconds_to_centuries(value) + + return value + + def mean_longitude_of_the_ascending_node_of_the_moon(time): """ This function computes Mean Longitude of the Ascending Node of the Moon per @@ -138,7 +258,8 @@ def mean_longitude_of_the_ascending_node_of_the_moon(time): :type time: float :param time: Terrestrial time measured in Julian centuries. - :return: + :return: Mean longitude of the Ascending Node of the Moon + :rtype: float """ # The polynomial is in arcseconds @@ -153,6 +274,35 @@ def mean_longitude_of_the_ascending_node_of_the_moon(time): return value +def mean_longitude_of_the_ascending_node_of_the_moon_derivative(time): + """ + This function computes Mean Longitude of the Ascending Node of the Moon + time derivative per IERS Conventions (2010). + + Note that technically, per IERS Conventions (2010), the input time should + be Barycentric Dynamical Time (TDB) but the difference between TDB and TT is + already small. Additionally, we do not consider effects outside the + Geocentric Celestial Reference System (GCRS) and the primary driver of the + TDB vs TT difference is earth's mean anomaly in its orbit. The error from + this simplication is less than a microarcsecond in nutation. + + :type time: float + :param time: Terrestrial time measured in Julian centuries. + :return: Mean Longitude of the Ascending Node of the Moon time derivative + :rtype: float + """ + + # The polynomial is in arcseconds + value = (-6962890.5431 + 14.9444 * time + 0.023106 * time ** 2 - + 0.00023756 * time ** 3) + + value = Conversions.arcsec_to_rad(value) + + value = Conversions.seconds_to_centuries(value) + + return value + + def mean_longitude_of_mercury(time): """ This function computes the mean longitude of Mercury per IERS Conventions @@ -167,7 +317,8 @@ def mean_longitude_of_mercury(time): :type time: float :param time: Terrestrial time measured in Julian centuries. - :return: + :return: Mean longitude of Mercury + :rtype: float """ # The polynomial is in radians @@ -179,6 +330,30 @@ def mean_longitude_of_mercury(time): return value +def mean_longitude_of_mercury_derivative(): + """ + This function computes the mean longitude of Mercury time derivative per + IERS Conventions (2010). + + Note that technically, per IERS Conventions (2010), the input time should + be Barycentric Dynamical Time (TDB) but the difference between TDB and TT is + already small. Additionally, we do not consider effects outside the + Geocentric Celestial Reference System (GCRS) and the primary driver of the + TDB vs TT difference is earth's mean anomaly in its orbit. The error from + this simplication is less than a microarcsecond in nutation. + + :return: Mean longitude of Mercury time derivative + :rtype: float + """ + + # The polynomial is in radians + value = 2608.7903141574 + + value = Conversions.seconds_to_centuries(value) + + return value + + def mean_longitude_of_venus(time): """ This function computes the mean longitude of Venus per IERS Conventions @@ -193,7 +368,8 @@ def mean_longitude_of_venus(time): :type time: float :param time: Terrestrial time measured in Julian centuries. - :return: + :return: Mean longitude of Venus + :rtype: float """ # The polynomial is in radians @@ -205,6 +381,30 @@ def mean_longitude_of_venus(time): return value +def mean_longitude_of_venus_derivative(): + """ + This function computes the mean longitude of Venus time derivative per IERS + Conventions (2010). + + Note that technically, per IERS Conventions (2010), the input time should + be Barycentric Dynamical Time (TDB) but the difference between TDB and TT is + already small. Additionally, we do not consider effects outside the + Geocentric Celestial Reference System (GCRS) and the primary driver of the + TDB vs TT difference is earth's mean anomaly in its orbit. The error from + this simplication is less than a microarcsecond in nutation. + + :return: Mean longitude of Venus time derivative + :rtype: float + """ + + # The polynomial is in radians + value = 1021.3285546211 + + value = Conversions.seconds_to_centuries(value) + + return value + + def mean_longitude_of_earth(time): """ This function computes the mean longitude of Earth per IERS Conventions @@ -219,7 +419,8 @@ def mean_longitude_of_earth(time): :type time: float :param time: Terrestrial time measured in Julian centuries. - :return: + :return: Mean longitude of Earth + :rtype: float """ # The polynomial is in radians @@ -231,6 +432,30 @@ def mean_longitude_of_earth(time): return value +def mean_longitude_of_earth_derivative(): + """ + This function computes the mean longitude of Earth time derivative per IERS + Conventions (2010). + + Note that technically, per IERS Conventions (2010), the input time should + be Barycentric Dynamical Time (TDB) but the difference between TDB and TT is + already small. Additionally, we do not consider effects outside the + Geocentric Celestial Reference System (GCRS) and the primary driver of the + TDB vs TT difference is earth's mean anomaly in its orbit. The error from + this simplication is less than a microarcsecond in nutation. + + :return: Mean longitude of Earth time derivative + :rtype: float + """ + + # The polynomial is in radians + value = 628.3075849991 + + value = Conversions.seconds_to_centuries(value) + + return value + + def mean_longitude_of_mars(time): """ This function computes the mean longitude of Mars per IERS Conventions @@ -245,7 +470,8 @@ def mean_longitude_of_mars(time): :type time: float :param time: Terrestrial time measured in Julian centuries. - :return: + :return: Mean longitude of Mars + :rtype: float """ # The polynomial is in radians @@ -257,6 +483,30 @@ def mean_longitude_of_mars(time): return value +def mean_longitude_of_mars_derivative(): + """ + This function computes the mean longitude of Mars time derivative per IERS + Conventions (2010). + + Note that technically, per IERS Conventions (2010), the input time should + be Barycentric Dynamical Time (TDB) but the difference between TDB and TT is + already small. Additionally, we do not consider effects outside the + Geocentric Celestial Reference System (GCRS) and the primary driver of the + TDB vs TT difference is earth's mean anomaly in its orbit. The error from + this simplication is less than a microarcsecond in nutation. + + :return: Mean longitude of Mars time derivative + :rtype: float + """ + + # The polynomial is in radians + value = 334.0612426700 + + value = Conversions.seconds_to_centuries(value) + + return value + + def mean_longitude_of_jupiter(time): """ This function computes the mean longitude of Jupiter per IERS Conventions @@ -271,7 +521,8 @@ def mean_longitude_of_jupiter(time): :type time: float :param time: Terrestrial time measured in Julian centuries. - :return: + :return: Mean longitude of Jupiter + :rtype: float """ # The polynomial is in radians @@ -283,6 +534,30 @@ def mean_longitude_of_jupiter(time): return value +def mean_longitude_of_jupiter_derivative(): + """ + This function computes the mean longitude of Jupiter time derivative per + IERS Conventions (2010). + + Note that technically, per IERS Conventions (2010), the input time should + be Barycentric Dynamical Time (TDB) but the difference between TDB and TT is + already small. Additionally, we do not consider effects outside the + Geocentric Celestial Reference System (GCRS) and the primary driver of the + TDB vs TT difference is earth's mean anomaly in its orbit. The error from + this simplication is less than a microarcsecond in nutation. + + :return: Mean longitude of Jupiter time derivative + :rtype: float + """ + + # The polynomial is in radians + value = 52.9690962641 + + value = Conversions.seconds_to_centuries(value) + + return value + + def mean_longitude_of_saturn(time): """ This function computes the mean longitude of Saturn per IERS Conventions @@ -297,7 +572,8 @@ def mean_longitude_of_saturn(time): :type time: float :param time: Terrestrial time measured in Julian centuries. - :return: + :return: Mean longitude of Saturn + :rtype: float """ # The polynomial is in radians @@ -309,6 +585,30 @@ def mean_longitude_of_saturn(time): return value +def mean_longitude_of_saturn_derivative(): + """ + This function computes the mean longitude of Saturn time derivative per + IERS Conventions (2010). + + Note that technically, per IERS Conventions (2010), the input time should + be Barycentric Dynamical Time (TDB) but the difference between TDB and TT is + already small. Additionally, we do not consider effects outside the + Geocentric Celestial Reference System (GCRS) and the primary driver of the + TDB vs TT difference is earth's mean anomaly in its orbit. The error from + this simplication is less than a microarcsecond in nutation. + + :return: Mean longitude of Saturn time derivative + :rtype: float + """ + + # The polynomial is in radians + value = 21.3299104960 + + value = Conversions.seconds_to_centuries(value) + + return value + + def mean_longitude_of_uranus(time): """ This function computes the mean longitude of Uranus per IERS Conventions @@ -323,7 +623,8 @@ def mean_longitude_of_uranus(time): :type time: float :param time: Terrestrial time measured in Julian centuries. - :return: + :return: Mean longitude of Uranus + :rtype: float """ # The polynomial is in radians @@ -335,6 +636,30 @@ def mean_longitude_of_uranus(time): return value +def mean_longitude_of_uranus_derivative(): + """ + This function computes the mean longitude of Uranus time derivative per + IERS Conventions (2010). + + Note that technically, per IERS Conventions (2010), the input time should + be Barycentric Dynamical Time (TDB) but the difference between TDB and TT is + already small. Additionally, we do not consider effects outside the + Geocentric Celestial Reference System (GCRS) and the primary driver of the + TDB vs TT difference is earth's mean anomaly in its orbit. The error from + this simplication is less than a microarcsecond in nutation. + + :return: Mean longitude of Uranus time derivative + :rtype: float + """ + + # The polynomial is in radians + value = 7.4781598567 + + value = Conversions.seconds_to_centuries(value) + + return value + + def mean_longitude_of_neptune(time): """ This function computes the mean longitude of Neptune per IERS Conventions @@ -349,7 +674,8 @@ def mean_longitude_of_neptune(time): :type time: float :param time: Terrestrial time measured in Julian centuries. - :return: + :return: Mean longitude of Neptune + :rtype: float """ # The polynomial is in radians @@ -361,6 +687,30 @@ def mean_longitude_of_neptune(time): return value +def mean_longitude_of_neptune_derivative(): + """ + This function computes the mean longitude of Neptune time derivative per + IERS Conventions (2010). + + Note that technically, per IERS Conventions (2010), the input time should + be Barycentric Dynamical Time (TDB) but the difference between TDB and TT is + already small. Additionally, we do not consider effects outside the + Geocentric Celestial Reference System (GCRS) and the primary driver of the + TDB vs TT difference is earth's mean anomaly in its orbit. The error from + this simplication is less than a microarcsecond in nutation. + + :return: Mean longitude of Neptune time derivative + :rtype: float + """ + + # The polynomial is in radians + value = 3.8133035638 + + value = Conversions.seconds_to_centuries(value) + + return value + + def general_precession_in_longitude(time): """ This function computes the general precession in longitude per IERS @@ -375,10 +725,37 @@ def general_precession_in_longitude(time): :type time: float :param time: Terrestrial time measured in Julian centuries. - :return: + :return: General precession in longitude + :rtype: float """ # The polynomial is in radians value = 0.02438175 * time + 0.00000538691 * time ** 2 return value + + +def general_precession_in_longitude_derivative(time): + """ + This function computes the general precession time derivative in longitude + per IERS Conventions (2010). + + Note that technically, per IERS Conventions (2010), the input time should + be Barycentric Dynamical Time (TDB) but the difference between TDB and TT is + already small. Additionally, we do not consider effects outside the + Geocentric Celestial Reference System (GCRS) and the primary driver of the + TDB vs TT difference is earth's mean anomaly in its orbit. The error from + this simplication is less than a microarcsecond in nutation. + + :type time: float + :param time: Terrestrial time measured in Julian centuries. + :return: General precession in longitude time derivative + :rtype: float + """ + + # The polynomial is in radians + value = 0.02438175 + 0.00001077382 * time + + value = Conversions.seconds_to_centuries(value) + + return value diff --git a/src/TerraFrame/PrecessionNutation/SeriesExpansion.py b/src/TerraFrame/PrecessionNutation/SeriesExpansion.py index 299adb6..84e1dda 100644 --- a/src/TerraFrame/PrecessionNutation/SeriesExpansion.py +++ b/src/TerraFrame/PrecessionNutation/SeriesExpansion.py @@ -66,12 +66,11 @@ def compute(self, t, derivative=False): for j in range(1, len(self._polynomial_coefficients)): poly_part_dt += self._polynomial_coefficients[j] * t ** (j - 1) - # TODO: Finish adding derivative stuff below - # Initialize all the argument parameters. "argument" is the term that # IERS uses to refer to the input to the trigonometric functions. The # order is tightly coupled with the file format. arguments = np.zeros((14, )) + d_arguments_dt = np.zeros((14, )) arguments[0] = Arguments.mean_anomaly_of_the_moon(t) # l arguments[1] = Arguments.mean_anomaly_of_the_sun(t) # l' @@ -89,7 +88,39 @@ def compute(self, t, derivative=False): arguments[12] = Arguments.mean_longitude_of_neptune(t) # L_Ne arguments[13] = Arguments.general_precession_in_longitude(t) # p_A + if derivative: + d_arguments_dt[0] = Arguments.mean_anomaly_of_the_moon_derivative(t) + d_arguments_dt[1] = Arguments.mean_anomaly_of_the_sun_derivative(t) + d_arguments_dt[2] = ( + Arguments + .mean_longitude_moon_minus_ascending_node_derivative(t) + ) + d_arguments_dt[3] = ( + Arguments + .mean_elongation_of_the_moon_from_the_sun_derivative(t) + ) + d_arguments_dt[4] = ( + Arguments + .mean_longitude_of_the_ascending_node_of_the_moon_derivative(t) + ) + d_arguments_dt[5] = Arguments.mean_longitude_of_mercury_derivative() + d_arguments_dt[6] = Arguments.mean_longitude_of_venus_derivative() + d_arguments_dt[7] = Arguments.mean_longitude_of_earth_derivative() + d_arguments_dt[8] = Arguments.mean_longitude_of_mars_derivative() + d_arguments_dt[9] = Arguments.mean_longitude_of_jupiter_derivative() + d_arguments_dt[10] = Arguments.mean_longitude_of_saturn_derivative() + d_arguments_dt[11] = Arguments.mean_longitude_of_uranus_derivative() + d_arguments_dt[12] = ( + Arguments.mean_longitude_of_neptune_derivative()) + d_arguments_dt[13] = ( + Arguments + .general_precession_in_longitude_derivative(t) + ) + + # TODO: Finish adding derivative stuff below + non_poly_part = 0.0 + d_non_poly_part_dt = 0.0 for row in self.data: j = row[0] diff --git a/src/TerraFrame/Utilities/Conversions.py b/src/TerraFrame/Utilities/Conversions.py index 65bd0f1..928051c 100644 --- a/src/TerraFrame/Utilities/Conversions.py +++ b/src/TerraFrame/Utilities/Conversions.py @@ -309,3 +309,19 @@ def days_to_centuries(value): """ return value / (365 * 100) + + +def seconds_to_centuries(value): + """ + This function takes a numeric input representing seconds and scales it to + centuries. The convertion assumes there is exactly 365 days in a year. + + :param value: Input value in seconds + :type value: int | float + :return: Output value in centuries + :rtype: int | float + """ + + value = days_to_centuries(seconds_to_days(value)) + + return value \ No newline at end of file From 7f86d51bc1f3acdb6ed3f4c0677f86c6b74ae8cf Mon Sep 17 00:00:00 2001 From: Caleb Morrison Date: Sat, 27 Sep 2025 14:45:23 -0700 Subject: [PATCH 05/17] Fix incorrect units in argument derivatives --- .../PrecessionNutation/Arguments.py | 29 ------------------- 1 file changed, 29 deletions(-) diff --git a/src/TerraFrame/PrecessionNutation/Arguments.py b/src/TerraFrame/PrecessionNutation/Arguments.py index 90e7767..2d8d4a8 100644 --- a/src/TerraFrame/PrecessionNutation/Arguments.py +++ b/src/TerraFrame/PrecessionNutation/Arguments.py @@ -62,8 +62,6 @@ def mean_anomaly_of_the_moon_derivative(time): value = Conversions.arcsec_to_rad(value) - value = Conversions.seconds_to_centuries(value) - return value @@ -121,8 +119,6 @@ def mean_anomaly_of_the_sun_derivative(time): value = Conversions.arcsec_to_rad(value) - value = Conversions.seconds_to_centuries(value) - return value @@ -180,8 +176,6 @@ def mean_longitude_moon_minus_ascending_node_derivative(time): value = Conversions.arcsec_to_rad(value) - value = Conversions.seconds_to_centuries(value) - return value @@ -239,8 +233,6 @@ def mean_elongation_of_the_moon_from_the_sun_derivative(time): value = Conversions.arcsec_to_rad(value) - value = Conversions.seconds_to_centuries(value) - return value @@ -298,8 +290,6 @@ def mean_longitude_of_the_ascending_node_of_the_moon_derivative(time): value = Conversions.arcsec_to_rad(value) - value = Conversions.seconds_to_centuries(value) - return value @@ -324,7 +314,6 @@ def mean_longitude_of_mercury(time): # The polynomial is in radians value = 4.402608842 + 2608.7903141574 * time - # Take the modulus before convertion to maintain accuracy value = math.fmod(value, 2.0 * math.pi) return value @@ -349,8 +338,6 @@ def mean_longitude_of_mercury_derivative(): # The polynomial is in radians value = 2608.7903141574 - value = Conversions.seconds_to_centuries(value) - return value @@ -400,8 +387,6 @@ def mean_longitude_of_venus_derivative(): # The polynomial is in radians value = 1021.3285546211 - value = Conversions.seconds_to_centuries(value) - return value @@ -451,8 +436,6 @@ def mean_longitude_of_earth_derivative(): # The polynomial is in radians value = 628.3075849991 - value = Conversions.seconds_to_centuries(value) - return value @@ -502,8 +485,6 @@ def mean_longitude_of_mars_derivative(): # The polynomial is in radians value = 334.0612426700 - value = Conversions.seconds_to_centuries(value) - return value @@ -553,8 +534,6 @@ def mean_longitude_of_jupiter_derivative(): # The polynomial is in radians value = 52.9690962641 - value = Conversions.seconds_to_centuries(value) - return value @@ -604,8 +583,6 @@ def mean_longitude_of_saturn_derivative(): # The polynomial is in radians value = 21.3299104960 - value = Conversions.seconds_to_centuries(value) - return value @@ -655,8 +632,6 @@ def mean_longitude_of_uranus_derivative(): # The polynomial is in radians value = 7.4781598567 - value = Conversions.seconds_to_centuries(value) - return value @@ -706,8 +681,6 @@ def mean_longitude_of_neptune_derivative(): # The polynomial is in radians value = 3.8133035638 - value = Conversions.seconds_to_centuries(value) - return value @@ -756,6 +729,4 @@ def general_precession_in_longitude_derivative(time): # The polynomial is in radians value = 0.02438175 + 0.00001077382 * time - value = Conversions.seconds_to_centuries(value) - return value From 47a2e8b207f441507c314216f8dedafad8436dbf Mon Sep 17 00:00:00 2001 From: Caleb Morrison Date: Sat, 27 Sep 2025 14:46:09 -0700 Subject: [PATCH 06/17] Add CIP series derivative calculations --- Tests/test_cip_calculation.py | 54 +++++++++++++++++++ .../PrecessionNutation/SeriesExpansion.py | 44 +++++++++++++-- 2 files changed, 94 insertions(+), 4 deletions(-) diff --git a/Tests/test_cip_calculation.py b/Tests/test_cip_calculation.py index e589830..c517ac8 100644 --- a/Tests/test_cip_calculation.py +++ b/Tests/test_cip_calculation.py @@ -7,6 +7,7 @@ from TerraFrame.PrecessionNutation import SeriesExpansion from TerraFrame.Utilities.Time import JulianDate +from TerraFrame.Utilities import Conversions def test_cip_calculation(): @@ -46,3 +47,56 @@ def test_cip_calculation(): assert np.max(np.abs(cip_x - cip_x_a)) < 1e-10 assert np.max(np.abs(cip_y - cip_y_a)) < 1e-10 assert np.max(np.abs(cip_s - cip_s_a)) < 1e-10 + + +def test_cip_derivatives(): + se_cip_x = SeriesExpansion.cip_x() + se_cip_y = SeriesExpansion.cip_y() + se_cip_sxy2 = SeriesExpansion.cip_sxy2() + + n = 10 + frac = np.linspace(0, 365.25 * 100, n) + + d_cip_x_dt = np.zeros((n,)) + d_cip_x_dt_fd = np.zeros((n,)) + d_cip_y_dt = np.zeros((n,)) + d_cip_y_dt_fd = np.zeros((n,)) + d_cip_s_dt = np.zeros((n,)) + d_cip_s_dt_fd = np.zeros((n,)) + + #dt = Conversions.seconds_to_centuries(1e-4) + dt = 1e-6 + + for i, val in enumerate(frac): + jd_tt = (JulianDate.JulianDate.j2000( + time_scale=JulianDate.TimeScales.TT) + val) + jdc_tt = JulianDate.julian_day_datetime_to_century_datetime(jd_tt) + + _, d_cip_x_dt[i] = se_cip_x.compute(jdc_tt, derivative=True) + _, d_cip_y_dt[i] = se_cip_y.compute(jdc_tt, derivative=True) + _, d_cip_s_dt[i] = se_cip_sxy2.compute(jdc_tt, derivative=True) + + f2 = se_cip_x.compute(jdc_tt + 1 * dt) + f1 = se_cip_x.compute(jdc_tt - 1 * dt) + + d_cip_x_dt_fd[i] = (f2 - f1) / (2 * dt) + # noinspection PyTypeChecker + d_cip_x_dt_fd[i] = Conversions.seconds_to_centuries(d_cip_x_dt_fd[i]) + + f2 = se_cip_y.compute(jdc_tt + 1 * dt) + f1 = se_cip_y.compute(jdc_tt - 1 * dt) + + d_cip_y_dt_fd[i] = (f2 - f1) / (2 * dt) + # noinspection PyTypeChecker + d_cip_y_dt_fd[i] = Conversions.seconds_to_centuries(d_cip_y_dt_fd[i]) + + f2 = se_cip_sxy2.compute(jdc_tt + 1 * dt) + f1 = se_cip_sxy2.compute(jdc_tt - 1 * dt) + + d_cip_s_dt_fd[i] = (f2 - f1) / (2 * dt) + # noinspection PyTypeChecker + d_cip_s_dt_fd[i] = Conversions.seconds_to_centuries(d_cip_s_dt_fd[i]) + + assert np.max(np.abs(d_cip_x_dt - d_cip_x_dt_fd)) < 1e-13 + assert np.max(np.abs(d_cip_y_dt - d_cip_y_dt_fd)) < 1e-13 + assert np.max(np.abs(d_cip_s_dt - d_cip_s_dt_fd)) < 1e-13 diff --git a/src/TerraFrame/PrecessionNutation/SeriesExpansion.py b/src/TerraFrame/PrecessionNutation/SeriesExpansion.py index 84e1dda..72ff18a 100644 --- a/src/TerraFrame/PrecessionNutation/SeriesExpansion.py +++ b/src/TerraFrame/PrecessionNutation/SeriesExpansion.py @@ -53,6 +53,20 @@ def __init__(self, data_file_path, polynomial_coefficients): self._polynomial_coefficients = polynomial_coefficients def compute(self, t, derivative=False): + """ + This function evaluates a celestial intermediate pole (CIP) series at + the given Julian century date. + + :param t: Time as a Julian Century Date (jdc) + :type t: JulianDate + :param derivative: Compute the time derivative of the series + (default false) + :type derivative: bool + :return: The computed series value at time t and optionally the series + time derivative + :rtype: float | float, float + """ + t = float(t) # units are micro-arcseconds @@ -117,8 +131,6 @@ def compute(self, t, derivative=False): .general_precession_in_longitude_derivative(t) ) - # TODO: Finish adding derivative stuff below - non_poly_part = 0.0 d_non_poly_part_dt = 0.0 @@ -129,13 +141,37 @@ def compute(self, t, derivative=False): arg = np.linalg.vecdot(row[4:], arguments) - non_poly_part += (a_s * np.sin(arg) + a_c * np.cos(arg)) * t ** j + arg_p = np.linalg.vecdot(row[4:], d_arguments_dt) + + tmp = a_s * np.sin(arg) + a_c * np.cos(arg) + + non_poly_part += tmp * t ** j + + if derivative: + tmp2 = (a_s * np.cos(arg) - a_c * np.sin(arg)) * arg_p + + # Avoid a numerical issue with the generic formula when j = 0 + if j > 0: + d_non_poly_part_dt += j * t ** (j - 1) * tmp + t ** j * tmp2 + else: + d_non_poly_part_dt += tmp2 + total = poly_part + non_poly_part + total_dt = poly_part_dt + d_non_poly_part_dt + total = Conversions.muas_to_rad(total) + total_dt = Conversions.muas_to_rad(total_dt) - return total + # Up till this point, we've been computing the derivative with + # respect to centuries. Change it to be with respect to seconds. + total_dt = Conversions.seconds_to_centuries(total_dt) + + if derivative: + return total, total_dt + else: + return total def cip_x(file_name=r'tab5.2a.txt'): From 6e1f4c7e5c1f22c11a89da98650a03b6967bf418 Mon Sep 17 00:00:00 2001 From: Caleb Morrison Date: Sat, 27 Sep 2025 15:19:23 -0700 Subject: [PATCH 07/17] Adjust readme and add r1 and r2 derivative tests --- README.md | 4 +-- Tests/test_matrix_derivatives.py | 50 ++++++++++++++++++++++++++++++-- 2 files changed, 49 insertions(+), 5 deletions(-) diff --git a/README.md b/README.md index 8b0c752..5987e04 100644 --- a/README.md +++ b/README.md @@ -45,8 +45,8 @@ post-processing. # License -This project - except for the IERS data files - is covered under the Mozilla -Public License Version 2.0 (MPL2). See the LICENSE.txt file for more +This project - except for the IERS and WGS84 data files - is covered under the +Mozilla Public License Version 2.0 (MPL2). See the LICENSE.txt file for more information. # Acknowledgements and References diff --git a/Tests/test_matrix_derivatives.py b/Tests/test_matrix_derivatives.py index 2ac5da0..87562f3 100644 --- a/Tests/test_matrix_derivatives.py +++ b/Tests/test_matrix_derivatives.py @@ -3,6 +3,50 @@ from TerraFrame.Utilities import TransformationMatrices +def test_r1_derivative(): + n = 25 + t = np.linspace(0, 24 * 60 * 60, n) + psi = [x ** 1.2 for x in t] + dpsidt = 1.2 + dt = 1e-4 + dpsi = dpsidt * dt + + error = np.zeros((n,)) + + for i in reversed(range(n)): + r1 = TransformationMatrices.r1(psi[i] - dpsi) + r2 = TransformationMatrices.r1(psi[i] + dpsi) + + drdt_e = (r2 - r1) / (2 * dt) + drdt = TransformationMatrices.dr1dt(psi[i], dpsidt) + + error[i] = np.max(np.abs(drdt_e - drdt)) + + assert (max(error) < 1e-6) + + +def test_r2_derivative(): + n = 25 + t = np.linspace(0, 24 * 60 * 60, n) + psi = [x ** 1.2 for x in t] + dpsidt = 1.2 + dt = 1e-4 + dpsi = dpsidt * dt + + error = np.zeros((n,)) + + for i in reversed(range(n)): + r1 = TransformationMatrices.r2(psi[i] - dpsi) + r2 = TransformationMatrices.r2(psi[i] + dpsi) + + drdt_e = (r2 - r1) / (2 * dt) + drdt = TransformationMatrices.dr2dt(psi[i], dpsidt) + + error[i] = np.max(np.abs(drdt_e - drdt)) + + assert (max(error) < 1e-6) + + def test_r3_derivative(): n = 25 t = np.linspace(0, 24 * 60 * 60, n) @@ -14,12 +58,12 @@ def test_r3_derivative(): error = np.zeros((n,)) for i in reversed(range(n)): - r1 = TransformationMatrices.r3(psi[i]) + r1 = TransformationMatrices.r3(psi[i] - dpsi) r2 = TransformationMatrices.r3(psi[i] + dpsi) - drdt_e = (r2 - r1) / dt + drdt_e = (r2 - r1) / (2 * dt) drdt = TransformationMatrices.dr3dt(psi[i], dpsidt) error[i] = np.max(np.abs(drdt_e - drdt)) - assert (max(error) < 1e-4) + assert (max(error) < 1e-6) From fda77c8486a02480d5257082b641f935b8badf42 Mon Sep 17 00:00:00 2001 From: Caleb Morrison Date: Sat, 27 Sep 2025 15:55:27 -0700 Subject: [PATCH 08/17] Fix t_cg derivative (e.g., Q'(t)) calc & add test --- Tests/test_cgrs_to_cirs.py | 45 ---------- Tests/test_cip_calculation.py | 1 - Tests/test_cirs_to_gcrs.py | 89 +++++++++++++++++++ src/TerraFrame/CelestialTerrestrial.py | 14 ++- .../Utilities/TransformationMatrices.py | 7 +- 5 files changed, 106 insertions(+), 50 deletions(-) delete mode 100644 Tests/test_cgrs_to_cirs.py create mode 100644 Tests/test_cirs_to_gcrs.py diff --git a/Tests/test_cgrs_to_cirs.py b/Tests/test_cgrs_to_cirs.py deleted file mode 100644 index 5e53025..0000000 --- a/Tests/test_cgrs_to_cirs.py +++ /dev/null @@ -1,45 +0,0 @@ -# This Source Code Form is subject to the terms of the Mozilla Public -# License, v. 2.0. If a copy of the MPL was not distributed with this -# file, You can obtain one at https://mozilla.org/MPL/2.0/. - -import random - -import erfa -import numpy as np - -from TerraFrame.PrecessionNutation import SeriesExpansion -from TerraFrame.Utilities import TransformationMatrices -from TerraFrame.Utilities.Time import JulianDate - - -def test_cirs_to_gcrs_calculation(): - se_cip_x = SeriesExpansion.cip_x() - se_cip_y = SeriesExpansion.cip_y() - se_cip_sxy2 = SeriesExpansion.cip_sxy2() - - val = random.uniform(0, 100.0) - jd_tt = (JulianDate.JulianDate.j2000( - time_scale=JulianDate.TimeScales.TT) + val) - jdc_tt = JulianDate.julian_day_datetime_to_century_datetime(jd_tt) - - cip_x = se_cip_x.compute(jdc_tt) - cip_y = se_cip_y.compute(jdc_tt) - sxy2 = se_cip_sxy2.compute(jdc_tt) - cip_s = sxy2 - cip_x * cip_y / 2.0 - - t_gc = TransformationMatrices.cirs_to_gcrs(cip_x, cip_y, cip_s) - - jd1, jd2 = jd_tt.integer_part(), jd_tt.fraction_part() - - # Get X, Y, s using IAU 2006/2000A model - x, y, s = erfa.xys06a(jd1, jd2) - - # ERFA/SOFA computes the inverse transform, so we need to take the transpose - t_cg_erfa = erfa.c2ixys(x, y, s) - t_gc_erfa = t_cg_erfa.T - - assert abs(cip_x - x) < 1e-8 - assert abs(cip_y - y) < 1e-8 - assert abs(cip_s - s) < 1e-8 - - assert (np.max(np.abs(t_gc - t_gc_erfa)) < 1e-10) diff --git a/Tests/test_cip_calculation.py b/Tests/test_cip_calculation.py index c517ac8..e4e3424 100644 --- a/Tests/test_cip_calculation.py +++ b/Tests/test_cip_calculation.py @@ -64,7 +64,6 @@ def test_cip_derivatives(): d_cip_s_dt = np.zeros((n,)) d_cip_s_dt_fd = np.zeros((n,)) - #dt = Conversions.seconds_to_centuries(1e-4) dt = 1e-6 for i, val in enumerate(frac): diff --git a/Tests/test_cirs_to_gcrs.py b/Tests/test_cirs_to_gcrs.py new file mode 100644 index 0000000..65e8cf2 --- /dev/null +++ b/Tests/test_cirs_to_gcrs.py @@ -0,0 +1,89 @@ +# This Source Code Form is subject to the terms of the Mozilla Public +# License, v. 2.0. If a copy of the MPL was not distributed with this +# file, You can obtain one at https://mozilla.org/MPL/2.0/. + +import random + +import erfa +import numpy as np + +from TerraFrame.PrecessionNutation import SeriesExpansion +from TerraFrame.Utilities import TransformationMatrices, Conversions +from TerraFrame.Utilities.Time import JulianDate + + +def test_cirs_to_gcrs_calculation(): + se_cip_x = SeriesExpansion.cip_x() + se_cip_y = SeriesExpansion.cip_y() + se_cip_sxy2 = SeriesExpansion.cip_sxy2() + + val = random.uniform(0, 100.0) + jd_tt = (JulianDate.JulianDate.j2000( + time_scale=JulianDate.TimeScales.TT) + val) + jdc_tt = JulianDate.julian_day_datetime_to_century_datetime(jd_tt) + + cip_x = se_cip_x.compute(jdc_tt) + cip_y = se_cip_y.compute(jdc_tt) + sxy2 = se_cip_sxy2.compute(jdc_tt) + cip_s = sxy2 - cip_x * cip_y / 2.0 + + t_gc = TransformationMatrices.cirs_to_gcrs(cip_x, cip_y, cip_s) + + jd1, jd2 = jd_tt.integer_part(), jd_tt.fraction_part() + + # Get X, Y, s using IAU 2006/2000A model + x, y, s = erfa.xys06a(jd1, jd2) + + # ERFA/SOFA computes the inverse transform, so we need to take the transpose + t_cg_erfa = erfa.c2ixys(x, y, s) + t_gc_erfa = t_cg_erfa.T + + assert abs(cip_x - x) < 1e-8 + assert abs(cip_y - y) < 1e-8 + assert abs(cip_s - s) < 1e-8 + + assert (np.max(np.abs(t_gc - t_gc_erfa)) < 1e-10) + + +def test_cirs_to_gcrs_derivative_calculation(): + se_cip_x = SeriesExpansion.cip_x() + se_cip_y = SeriesExpansion.cip_y() + se_cip_sxy2 = SeriesExpansion.cip_sxy2() + + dt = 1e-6 + + val = random.uniform(0, 100.0) + jd_tt = (JulianDate.JulianDate.j2000( + time_scale=JulianDate.TimeScales.TT) + val) + jdc_tt = JulianDate.julian_day_datetime_to_century_datetime(jd_tt) + + cip_x, d_cip_x_dt = se_cip_x.compute(jdc_tt, derivative=True) + cip_y, d_cip_y_dt = se_cip_y.compute(jdc_tt, derivative=True) + sxy2, d_cip_sxy2_dt = se_cip_sxy2.compute(jdc_tt, derivative=True) + cip_s = sxy2 - cip_x * cip_y / 2.0 + d_cip_s_dt = (d_cip_sxy2_dt - 0.5 * cip_y * d_cip_x_dt - + 0.5 * cip_x * d_cip_y_dt) + + d_t_gc_dt = (TransformationMatrices. + cirs_to_gcrs_derivative(cip_x, cip_y, cip_s, + d_cip_x_dt, d_cip_y_dt, d_cip_s_dt)) + + cip_x2 = se_cip_x.compute(jdc_tt + dt) + cip_y2 = se_cip_y.compute(jdc_tt + dt) + sxy22 = se_cip_sxy2.compute(jdc_tt + dt) + cip_s2 = sxy22 - cip_x2 * cip_y2 / 2.0 + + cip_x1 = se_cip_x.compute(jdc_tt - dt) + cip_y1 = se_cip_y.compute(jdc_tt - dt) + sxy21 = se_cip_sxy2.compute(jdc_tt - dt) + cip_s1 = sxy21 - cip_x1 * cip_y1 / 2.0 + + t_gc2 = TransformationMatrices.cirs_to_gcrs(cip_x2, cip_y2, cip_s2) + t_gc1 = TransformationMatrices.cirs_to_gcrs(cip_x1, cip_y1, cip_s1) + + d_t_gc_dt_fd = (t_gc2 - t_gc1) / (2 * dt) + d_t_gc_dt_fd *= Conversions.seconds_to_centuries(1.0) + + error = np.abs(d_t_gc_dt - d_t_gc_dt_fd) + + assert np.max(error) < 1e-15 diff --git a/src/TerraFrame/CelestialTerrestrial.py b/src/TerraFrame/CelestialTerrestrial.py index e3e6693..588003d 100644 --- a/src/TerraFrame/CelestialTerrestrial.py +++ b/src/TerraFrame/CelestialTerrestrial.py @@ -30,6 +30,7 @@ def __init__(self, user_polar_motion=True, user_nutation_corrections=True): self.t_gc = None self.t_ct = None self.t_ti = None + self.angular_velocity = None def itrs_to_gcrs(self, time): if isinstance(time, datetime.datetime): @@ -170,7 +171,9 @@ def itrs_to_gcrs_angular_vel(self, time): dcip_y_dt += ddy_dt # Create the first transformation matrix + # Q(t) t_gc = TransformationMatrices.cirs_to_gcrs(cip_x, cip_y, cip_s) + # Q'(t) dt_gc_dt = (TransformationMatrices .cirs_to_gcrs_derivative(cip_x, cip_y, cip_s, dcip_x_dt, dcip_y_dt, dcip_s_dt)) @@ -179,8 +182,10 @@ def itrs_to_gcrs_angular_vel(self, time): # Intermediate Reference System (TIRS) to the Celestial Intermediate # Reference System (CIRS): TIRS -> CIRS. # This function uses normal JD time in UT1. + # R(t) t_ct = TransformationMatrices.earth_rotation_matrix(jd_ut1) + # R'(t) dt_ct_dt = (TransformationMatrices .earth_rotation_matrix_derivative(jd_ut1)) @@ -206,17 +211,24 @@ def itrs_to_gcrs_angular_vel(self, time): pm_x = TerraFrame.Utilities.Conversions.arcsec_to_rad(pm_x) pm_y = TerraFrame.Utilities.Conversions.arcsec_to_rad(pm_y) + # W(t) t_ti = TransformationMatrices.itrs_to_tirs(pm_x, pm_y, sp) + # W'(t) dt_ti_dt = TransformationMatrices.itrs_to_tirs_derivative( pm_x, pm_y, sp, dpm_x_dt, dpm_y_dt, dsp_dt) # Construct the final transformation matrix: ITRS -> GCRS + # GCRS = Q(t) * R(t) * W(t) t_gi = t_gc @ t_ct @ t_ti + angular_vel = (t_gc @ t_ct @ dt_ti_dt + t_gc @ dt_ct_dt @ t_ti + + dt_ct_dt @ t_ct @ t_ti) + self.t_gi = t_gi self.t_gc = t_gc self.t_ct = t_ct self.t_ti = t_ti + self.angular_velocity = angular_vel - return t_gi + return t_gi, angular_vel diff --git a/src/TerraFrame/Utilities/TransformationMatrices.py b/src/TerraFrame/Utilities/TransformationMatrices.py index 5ad5def..35dcd8f 100644 --- a/src/TerraFrame/Utilities/TransformationMatrices.py +++ b/src/TerraFrame/Utilities/TransformationMatrices.py @@ -391,9 +391,10 @@ def cirs_to_gcrs_derivative(x, y, s, dx_dt, dy_dt, ds_dt): dd_dt = (np.sqrt(-1 - 1/(x**2 + y**2 - 1)) * (x * dx_dt + y * dy_dt) / (x**2 + y**2)) - dt_gc_dt = (r3(-e) @ r2(-d) @ dr3dt(s, ds_dt) + - r3(-e) @ dr2dt(-d, -dd_dt) @ r3(s) + - dr3dt(-e, -de_dt) @ r2(-d) @ r3(s)) + dt_gc_dt = (r3(-e) @ r2(-d) @ r3(e) @ dr3dt(s, ds_dt) + + r3(-e) @ r2(-d) @ dr3dt(e, de_dt) @ r3(s) + + r3(-e) @ dr2dt(-d, -dd_dt) @ r3(e) @ r3(s) + + dr3dt(-e, -de_dt) @ r2(-d) @ r3(e) @ r3(s)) return dt_gc_dt From 129750e35fab6d4992b680a809ec3f5b4d8e5623 Mon Sep 17 00:00:00 2001 From: Caleb Morrison Date: Sat, 27 Sep 2025 16:15:47 -0700 Subject: [PATCH 09/17] Add earth rotation derivative matrix test --- Tests/test_cip_calculation.py | 6 +++--- Tests/test_cirs_to_gcrs.py | 2 +- Tests/test_era_calculations.py | 22 +++++++++++++++++++++- 3 files changed, 25 insertions(+), 5 deletions(-) diff --git a/Tests/test_cip_calculation.py b/Tests/test_cip_calculation.py index e4e3424..8c09567 100644 --- a/Tests/test_cip_calculation.py +++ b/Tests/test_cip_calculation.py @@ -96,6 +96,6 @@ def test_cip_derivatives(): # noinspection PyTypeChecker d_cip_s_dt_fd[i] = Conversions.seconds_to_centuries(d_cip_s_dt_fd[i]) - assert np.max(np.abs(d_cip_x_dt - d_cip_x_dt_fd)) < 1e-13 - assert np.max(np.abs(d_cip_y_dt - d_cip_y_dt_fd)) < 1e-13 - assert np.max(np.abs(d_cip_s_dt - d_cip_s_dt_fd)) < 1e-13 + assert (np.max(np.abs(d_cip_x_dt - d_cip_x_dt_fd)) < 1e-13) + assert (np.max(np.abs(d_cip_y_dt - d_cip_y_dt_fd)) < 1e-13) + assert (np.max(np.abs(d_cip_s_dt - d_cip_s_dt_fd)) < 1e-13) diff --git a/Tests/test_cirs_to_gcrs.py b/Tests/test_cirs_to_gcrs.py index 65e8cf2..1853119 100644 --- a/Tests/test_cirs_to_gcrs.py +++ b/Tests/test_cirs_to_gcrs.py @@ -86,4 +86,4 @@ def test_cirs_to_gcrs_derivative_calculation(): error = np.abs(d_t_gc_dt - d_t_gc_dt_fd) - assert np.max(error) < 1e-15 + assert (np.max(error) < 1e-15) diff --git a/Tests/test_era_calculations.py b/Tests/test_era_calculations.py index 1d9693e..b61eb63 100644 --- a/Tests/test_era_calculations.py +++ b/Tests/test_era_calculations.py @@ -8,7 +8,7 @@ import numpy as np import TerraFrame.Earth -from TerraFrame.Utilities import TransformationMatrices +from TerraFrame.Utilities import TransformationMatrices, Conversions from TerraFrame.Utilities.Time import JulianDate @@ -38,3 +38,23 @@ def test_era_transformation_calculation(): t_e_erfa = TransformationMatrices.r3(-era_a) assert (np.max(np.abs(t_e - t_e_erfa)) < 1e-10) + + +def test_era_transformation_derivative_calculation(): + val = random.uniform(0, 100.0) + jd_ut1 = (JulianDate.JulianDate.j2000( + time_scale=JulianDate.TimeScales.UT1) + val) + + dt = 1e-6 + + t_e2 = TransformationMatrices.earth_rotation_matrix(jd_ut1 + dt) + t_e1 = TransformationMatrices.earth_rotation_matrix(jd_ut1 - dt) + + d_te_dt_fd = (t_e2 - t_e1) / (2 * dt) + d_te_dt_fd *= Conversions.seconds_to_days(1) + + d_te_dt = TransformationMatrices.earth_rotation_matrix_derivative(jd_ut1) + + error = np.abs(d_te_dt - d_te_dt_fd) + + assert (np.max(error) < 1e-13) From 4f8f5e207a880e3745c9f0bc24f39fad4bb62484 Mon Sep 17 00:00:00 2001 From: Caleb Morrison Date: Sat, 27 Sep 2025 16:40:13 -0700 Subject: [PATCH 10/17] Add t_it derivative test and unit error The pm_x and pm_y derivatives were accidentally in terms of days and not seconds. This commit changes the derivatives to be in terms of seconds. --- Tests/test_polar_motion.py | 48 ++++++++++++++++++- src/TerraFrame/CelestialTerrestrial.py | 2 + .../Utilities/TransformationMatrices.py | 3 +- 3 files changed, 50 insertions(+), 3 deletions(-) diff --git a/Tests/test_polar_motion.py b/Tests/test_polar_motion.py index fef7f29..69fc14f 100644 --- a/Tests/test_polar_motion.py +++ b/Tests/test_polar_motion.py @@ -7,7 +7,8 @@ import erfa import numpy as np -from TerraFrame.Utilities import TransformationMatrices +from TerraFrame.Utilities import (TransformationMatrices, BulletinData, + Conversions) from TerraFrame.Utilities.Time import JulianDate @@ -46,3 +47,48 @@ def test_itrs_to_tirs(): t_ti_erfa = t_it_erfa.T assert (np.max(np.abs(t_ti - t_ti_erfa)) < 1e-10) + + +def test_itrs_to_tirs_derivative_calculation(): + val = random.uniform(0, 100.0) + dt = 1e-6 + dtc = Conversions.seconds_to_centuries(dt) + bd = BulletinData.BulletinData() + + jd_utc = (JulianDate.JulianDate.j2000( + time_scale=JulianDate.TimeScales.UTC) + val) + mjd_utc = JulianDate.julian_date_to_modified_julian_date(jd_utc) + + jd_tt = Conversions.any_to_tt(jd_utc) + jdc_tt = JulianDate.julian_day_datetime_to_century_datetime(jd_tt) + + pm_x = bd.f_pm_x(float(mjd_utc)) + pm_y = bd.f_pm_y(float(mjd_utc)) + dpm_x_dt = bd.f_pm_x(float(mjd_utc), derivative=True) + dpm_x_dt = Conversions.seconds_to_days(dpm_x_dt) + dpm_y_dt = bd.f_pm_y(float(mjd_utc), derivative=True) + dpm_y_dt = Conversions.seconds_to_days(dpm_y_dt) + + sp = TransformationMatrices.calculate_s_prime(jdc_tt) + dsp_dt = TransformationMatrices.calculate_s_prime_derivative() + + dt_ti_dt = TransformationMatrices.itrs_to_tirs_derivative(pm_x, pm_y, sp, + dpm_x_dt, dpm_y_dt, dsp_dt) + + pm_x2 = bd.f_pm_x(float(mjd_utc) + dt) + pm_y2 = bd.f_pm_y(float(mjd_utc) + dt) + sp2 = TransformationMatrices.calculate_s_prime(jdc_tt + dtc) + + pm_x1 = bd.f_pm_x(float(mjd_utc) - dt) + pm_y1 = bd.f_pm_y(float(mjd_utc) - dt) + sp1 = TransformationMatrices.calculate_s_prime(jdc_tt - dtc) + + f2 = TransformationMatrices.itrs_to_tirs(pm_x2, pm_y2, sp2) + f1 = TransformationMatrices.itrs_to_tirs(pm_x1, pm_y1, sp1) + + dt_ti_dt_fd = (f2 - f1) / (2 * dt) + dt_ti_dt_fd *= Conversions.seconds_to_days(1.0) + + error = np.abs(dt_ti_dt - dt_ti_dt_fd) + + assert (np.max(error) < 1e-14) diff --git a/src/TerraFrame/CelestialTerrestrial.py b/src/TerraFrame/CelestialTerrestrial.py index 588003d..5d5f1a5 100644 --- a/src/TerraFrame/CelestialTerrestrial.py +++ b/src/TerraFrame/CelestialTerrestrial.py @@ -198,7 +198,9 @@ def itrs_to_gcrs_angular_vel(self, time): pm_x = self.bd.f_pm_x(float(mjd_utc)) pm_y = self.bd.f_pm_y(float(mjd_utc)) dpm_x_dt = self.bd.f_pm_x(float(mjd_utc), derivative=True) + dpm_x_dt = Conversions.seconds_to_days(dpm_x_dt) dpm_y_dt = self.bd.f_pm_y(float(mjd_utc), derivative=True) + dpm_y_dt = Conversions.seconds_to_days(dpm_y_dt) else: pm_x = 0.0 pm_y = 0.0 diff --git a/src/TerraFrame/Utilities/TransformationMatrices.py b/src/TerraFrame/Utilities/TransformationMatrices.py index 35dcd8f..bd55e5b 100644 --- a/src/TerraFrame/Utilities/TransformationMatrices.py +++ b/src/TerraFrame/Utilities/TransformationMatrices.py @@ -313,8 +313,7 @@ def calculate_s_prime_derivative(): # This is an approximation good for the next century. See section 5.5.2 of # IERS Conventions (2010) for more context. - dt_s_prime_dt = Conversions.days_to_centuries(-47e-6) - dt_s_prime_dt = Conversions.seconds_to_days(dt_s_prime_dt) + dt_s_prime_dt = Conversions.seconds_to_centuries(-47e-6) dt_s_prime_dt = Conversions.arcsec_to_rad(dt_s_prime_dt) return dt_s_prime_dt From bccb307edb92ca75f47091dddf179bd7f06cbd1e Mon Sep 17 00:00:00 2001 From: Caleb Morrison Date: Sat, 27 Sep 2025 17:30:07 -0700 Subject: [PATCH 11/17] Finish angular velocity calc and add test Additionally, fix a couple of typos and a unit error in the t_ti derivative test. --- Tests/test_itrs_to_gcrs.py | 43 +++++++++++++++++----- Tests/test_polar_motion.py | 8 +++++ src/TerraFrame/CelestialTerrestrial.py | 50 ++++++++++++++------------ 3 files changed, 71 insertions(+), 30 deletions(-) diff --git a/Tests/test_itrs_to_gcrs.py b/Tests/test_itrs_to_gcrs.py index 8ef8667..98a6a5c 100644 --- a/Tests/test_itrs_to_gcrs.py +++ b/Tests/test_itrs_to_gcrs.py @@ -26,8 +26,9 @@ def test_itrs_to_gcrs_calculation(): jd_ut1 = Conversions.utc_to_ut1(jd_utc) # No corrections since we're comparing against a bare ERFA routine - ct = (TerraFrame.CelestialTerrestrialTransformation(user_polar_motion=False, - user_nutation_corrections=False)) + ct = (TerraFrame. + CelestialTerrestrialTransformation(use_polar_motion=False, + use_nutation_corrections=False)) t_ig = ct.itrs_to_gcrs(jd_tt) @@ -49,8 +50,9 @@ def test_against_astropy(): # Apparently, Astropy doesn't use the dx and dy nutation corrections. # See: https://github.com/astropy/astropy/issues/11110 - ct = (TerraFrame.CelestialTerrestrialTransformation(user_polar_motion=True, - user_nutation_corrections=False)) + ct = (TerraFrame. + CelestialTerrestrialTransformation(use_polar_motion=True, + use_nutation_corrections=False)) t_ig = ct.gcrs_to_itrs(jd_utc) @@ -83,8 +85,9 @@ def test_leap_second_against_astropy(): # Apparently, Astropy doesn't use the dx and dy nutation corrections. # See: https://github.com/astropy/astropy/issues/11110 - ct = (TerraFrame.CelestialTerrestrialTransformation(user_polar_motion=True, - user_nutation_corrections=False)) + ct = (TerraFrame. + CelestialTerrestrialTransformation(use_polar_motion=True, + use_nutation_corrections=False)) bd = BulletinData.BulletinData() file_path = bd.dat_file_path() @@ -104,9 +107,33 @@ def test_leap_second_against_astropy(): t = astropy.time.Time(jd_tt.integer_part(), jd_tt.fraction_part(), format='jd', scale='tt') - gcrs_basis_astro = GCRS(CartesianRepresentation(gcrs_basis, unit=u.one), - obstime=t) + gcrs_basis_astro = GCRS(CartesianRepresentation(gcrs_basis, + unit=u.one), obstime=t) itrs_basis_astro = gcrs_basis_astro.transform_to(ITRS(obstime=t)) itrs_basis_astro = itrs_basis_astro.cartesian.xyz.to_value() assert (np.max(np.abs(itrs_basis - itrs_basis_astro)) < 1e-10) + +def test_t_gi_angular_velocity(): + val = random.uniform(0, 100.0) + jd_tt = (JulianDate.JulianDate.j2000( + time_scale=JulianDate.TimeScales.TT) + val) + + dt = 1e-6 + + ct = TerraFrame.CelestialTerrestrialTransformation(use_polar_motion=True, + use_nutation_corrections=True) + + f2, _ = ct.itrs_to_gcrs_angular_vel(jd_tt + dt) + f1, _ = ct.itrs_to_gcrs_angular_vel(jd_tt - dt) + + d_t_gi_dt_fd = (f2 - f1) / (2 * dt) + d_t_gi_dt_fd *= Conversions.seconds_to_days(1.0) + + t_gi, d_t_gi_dt = ct.itrs_to_gcrs_angular_vel(jd_tt) + + d_t_gi_dt_fd = d_t_gi_dt_fd.T @ t_gi + + error = np.abs(d_t_gi_dt - d_t_gi_dt_fd) + + assert (np.max(error) < 1e-12) diff --git a/Tests/test_polar_motion.py b/Tests/test_polar_motion.py index 69fc14f..af8d552 100644 --- a/Tests/test_polar_motion.py +++ b/Tests/test_polar_motion.py @@ -64,9 +64,13 @@ def test_itrs_to_tirs_derivative_calculation(): pm_x = bd.f_pm_x(float(mjd_utc)) pm_y = bd.f_pm_y(float(mjd_utc)) + pm_x = Conversions.arcsec_to_rad(pm_x) + pm_y = Conversions.arcsec_to_rad(pm_y) dpm_x_dt = bd.f_pm_x(float(mjd_utc), derivative=True) + dpm_x_dt = Conversions.arcsec_to_rad(dpm_x_dt) dpm_x_dt = Conversions.seconds_to_days(dpm_x_dt) dpm_y_dt = bd.f_pm_y(float(mjd_utc), derivative=True) + dpm_y_dt = Conversions.arcsec_to_rad(dpm_y_dt) dpm_y_dt = Conversions.seconds_to_days(dpm_y_dt) sp = TransformationMatrices.calculate_s_prime(jdc_tt) @@ -77,10 +81,14 @@ def test_itrs_to_tirs_derivative_calculation(): pm_x2 = bd.f_pm_x(float(mjd_utc) + dt) pm_y2 = bd.f_pm_y(float(mjd_utc) + dt) + pm_x2 = Conversions.arcsec_to_rad(pm_x2) + pm_y2 = Conversions.arcsec_to_rad(pm_y2) sp2 = TransformationMatrices.calculate_s_prime(jdc_tt + dtc) pm_x1 = bd.f_pm_x(float(mjd_utc) - dt) pm_y1 = bd.f_pm_y(float(mjd_utc) - dt) + pm_x1 = Conversions.arcsec_to_rad(pm_x1) + pm_y1 = Conversions.arcsec_to_rad(pm_y1) sp1 = TransformationMatrices.calculate_s_prime(jdc_tt - dtc) f2 = TransformationMatrices.itrs_to_tirs(pm_x2, pm_y2, sp2) diff --git a/src/TerraFrame/CelestialTerrestrial.py b/src/TerraFrame/CelestialTerrestrial.py index 5d5f1a5..790f0f4 100644 --- a/src/TerraFrame/CelestialTerrestrial.py +++ b/src/TerraFrame/CelestialTerrestrial.py @@ -12,15 +12,15 @@ class CelestialTerrestrialTransformation: - def __init__(self, user_polar_motion=True, user_nutation_corrections=True): + def __init__(self, use_polar_motion=True, use_nutation_corrections=True): self.se_cip_x = SeriesExpansion.cip_x() self.se_cip_y = SeriesExpansion.cip_y() self.se_cip_sxy2 = SeriesExpansion.cip_sxy2() - self._user_polar_motion = user_polar_motion - self._user_nutation_corrections = user_nutation_corrections + self._use_polar_motion = use_polar_motion + self._use_nutation_corrections = use_nutation_corrections - if not self._user_polar_motion and not self._user_nutation_corrections: + if not self._use_polar_motion and not self._use_nutation_corrections: self.bd = None else: self.bd = BulletinData.BulletinData() @@ -68,7 +68,7 @@ def itrs_to_gcrs(self, time): # Reference System (CIRS) to Geocentric Celestial Reference System # (GCRS) matrix: CIRS -> GCRS. # Get corrections by interpolating in the IERS Bulletin A data - if self._user_nutation_corrections: + if self._use_nutation_corrections: dx = self.bd.f_nc_dx(float(mjd_utc)) dy = self.bd.f_nc_dy(float(mjd_utc)) @@ -89,7 +89,7 @@ def itrs_to_gcrs(self, time): # Terrestrial Reference System (ITRS) to Terrestrial Intermediate # Reference System (TIRS) transformation matrix can be constructed: # ITRS -> TIRS. - if self._user_polar_motion: + if self._use_polar_motion: pm_x = self.bd.f_pm_x(float(mjd_utc)) pm_y = self.bd.f_pm_y(float(mjd_utc)) else: @@ -118,7 +118,6 @@ def gcrs_to_itrs(self, time): return t_gi.T - def itrs_to_gcrs_angular_vel(self, time): if isinstance(time, datetime.datetime): time = Time.JulianDate.julian_date_from_pydatetime(time) @@ -157,7 +156,7 @@ def itrs_to_gcrs_angular_vel(self, time): # Reference System (CIRS) to Geocentric Celestial Reference System # (GCRS) matrix: CIRS -> GCRS. # Get corrections by interpolating in the IERS Bulletin A data - if self._user_nutation_corrections: + if self._use_nutation_corrections: dx = self.bd.f_nc_dx(float(mjd_utc)) dy = self.bd.f_nc_dy(float(mjd_utc)) @@ -165,7 +164,11 @@ def itrs_to_gcrs_angular_vel(self, time): cip_y += Conversions.mas_to_rad(dy) ddx_dt = self.bd.f_nc_dx(float(mjd_utc), derivative=True) + ddx_dt = Conversions.mas_to_rad(ddx_dt) + ddx_dt = Conversions.seconds_to_days(ddx_dt) ddy_dt = self.bd.f_nc_dy(float(mjd_utc), derivative=True) + ddy_dt = Conversions.mas_to_rad(ddy_dt) + ddy_dt = Conversions.seconds_to_days(ddy_dt) dcip_x_dt += ddx_dt dcip_y_dt += ddy_dt @@ -174,9 +177,10 @@ def itrs_to_gcrs_angular_vel(self, time): # Q(t) t_gc = TransformationMatrices.cirs_to_gcrs(cip_x, cip_y, cip_s) # Q'(t) - dt_gc_dt = (TransformationMatrices - .cirs_to_gcrs_derivative(cip_x, cip_y, cip_s, - dcip_x_dt, dcip_y_dt, dcip_s_dt)) + dt_gc_dt = ( + TransformationMatrices.cirs_to_gcrs_derivative(cip_x, cip_y, cip_s, + dcip_x_dt, dcip_y_dt, + dcip_s_dt)) # The Earth rotation matrix is the transformation from the Terrestrial # Intermediate Reference System (TIRS) to the Celestial Intermediate @@ -186,20 +190,24 @@ def itrs_to_gcrs_angular_vel(self, time): t_ct = TransformationMatrices.earth_rotation_matrix(jd_ut1) # R'(t) - dt_ct_dt = (TransformationMatrices - .earth_rotation_matrix_derivative(jd_ut1)) + dt_ct_dt = ( + TransformationMatrices.earth_rotation_matrix_derivative(jd_ut1)) # Given polar motion offsets pm_x and pm_y, along with the Terrestrial # Intermediate Origin (TIO) locator (s prime or sp), the International # Terrestrial Reference System (ITRS) to Terrestrial Intermediate # Reference System (TIRS) transformation matrix can be constructed: # ITRS -> TIRS. - if self._user_polar_motion: + if self._use_polar_motion: pm_x = self.bd.f_pm_x(float(mjd_utc)) pm_y = self.bd.f_pm_y(float(mjd_utc)) + pm_x = Conversions.arcsec_to_rad(pm_x) + pm_y = Conversions.arcsec_to_rad(pm_y) dpm_x_dt = self.bd.f_pm_x(float(mjd_utc), derivative=True) + dpm_x_dt = Conversions.arcsec_to_rad(dpm_x_dt) dpm_x_dt = Conversions.seconds_to_days(dpm_x_dt) dpm_y_dt = self.bd.f_pm_y(float(mjd_utc), derivative=True) + dpm_y_dt = Conversions.arcsec_to_rad(dpm_y_dt) dpm_y_dt = Conversions.seconds_to_days(dpm_y_dt) else: pm_x = 0.0 @@ -210,22 +218,20 @@ def itrs_to_gcrs_angular_vel(self, time): sp = TransformationMatrices.calculate_s_prime(jdc_tt) dsp_dt = TransformationMatrices.calculate_s_prime_derivative() - pm_x = TerraFrame.Utilities.Conversions.arcsec_to_rad(pm_x) - pm_y = TerraFrame.Utilities.Conversions.arcsec_to_rad(pm_y) - # W(t) t_ti = TransformationMatrices.itrs_to_tirs(pm_x, pm_y, sp) # W'(t) - dt_ti_dt = TransformationMatrices.itrs_to_tirs_derivative( - pm_x, pm_y, sp, - dpm_x_dt, dpm_y_dt, dsp_dt) + dt_ti_dt = TransformationMatrices.itrs_to_tirs_derivative(pm_x, pm_y, + sp, dpm_x_dt, dpm_y_dt, dsp_dt) # Construct the final transformation matrix: ITRS -> GCRS # GCRS = Q(t) * R(t) * W(t) t_gi = t_gc @ t_ct @ t_ti - angular_vel = (t_gc @ t_ct @ dt_ti_dt + t_gc @ dt_ct_dt @ t_ti + - dt_ct_dt @ t_ct @ t_ti) + d_t_gi_dt = (t_gc @ t_ct @ dt_ti_dt + t_gc @ dt_ct_dt @ t_ti + + dt_gc_dt @ t_ct @ t_ti) + + angular_vel = d_t_gi_dt.T @ t_gi self.t_gi = t_gi self.t_gc = t_gc From 499a0d9680c5dd8afaf062f894f124f72a783b08 Mon Sep 17 00:00:00 2001 From: Caleb Morrison Date: Sat, 27 Sep 2025 17:38:15 -0700 Subject: [PATCH 12/17] Adjust angular velocity test assert threshold --- Tests/test_itrs_to_gcrs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Tests/test_itrs_to_gcrs.py b/Tests/test_itrs_to_gcrs.py index 98a6a5c..08432fb 100644 --- a/Tests/test_itrs_to_gcrs.py +++ b/Tests/test_itrs_to_gcrs.py @@ -136,4 +136,4 @@ def test_t_gi_angular_velocity(): error = np.abs(d_t_gi_dt - d_t_gi_dt_fd) - assert (np.max(error) < 1e-12) + assert (np.max(error) < 1e-11) From a2e512c371074eea4e313749d90fbeda06fe76dd Mon Sep 17 00:00:00 2001 From: Caleb Morrison Date: Sun, 28 Sep 2025 08:53:26 -0700 Subject: [PATCH 13/17] Add missing license headers --- Tests/test_lat_lot_cartesian.py | 4 ++++ Tests/test_matrix_derivatives.py | 4 ++++ src/TerraFrame/Earth.py | 2 +- 3 files changed, 9 insertions(+), 1 deletion(-) diff --git a/Tests/test_lat_lot_cartesian.py b/Tests/test_lat_lot_cartesian.py index 3116978..391ee34 100644 --- a/Tests/test_lat_lot_cartesian.py +++ b/Tests/test_lat_lot_cartesian.py @@ -1,3 +1,7 @@ +# This Source Code Form is subject to the terms of the Mozilla Public +# License, v. 2.0. If a copy of the MPL was not distributed with this +# file, You can obtain one at https://mozilla.org/MPL/2.0/. + import math import random diff --git a/Tests/test_matrix_derivatives.py b/Tests/test_matrix_derivatives.py index 87562f3..c9d0266 100644 --- a/Tests/test_matrix_derivatives.py +++ b/Tests/test_matrix_derivatives.py @@ -1,3 +1,7 @@ +# This Source Code Form is subject to the terms of the Mozilla Public +# License, v. 2.0. If a copy of the MPL was not distributed with this +# file, You can obtain one at https://mozilla.org/MPL/2.0/. + import numpy as np from TerraFrame.Utilities import TransformationMatrices diff --git a/src/TerraFrame/Earth.py b/src/TerraFrame/Earth.py index 7181d01..5370b63 100644 --- a/src/TerraFrame/Earth.py +++ b/src/TerraFrame/Earth.py @@ -174,7 +174,7 @@ def lat_lon_alt_from_cartesian(self, x, y, z): c = cn1 # It's symbol soup. No way around it. If you want this to make - # sense, you'll have to look at the source paper: Fukushima 2006. + # sense, you'll have to look at the reference paper: Fukushima 2006. a = math.sqrt(s ** 2 + c ** 2) b = 1.5 * e * s * c ** 2 * ((p * s - zb * c) * a - e * s * c) f = p * a ** 3 - e * c ** 3 From 5a98e4c3dfb08eb625abf116f6771f0cfa1a2742 Mon Sep 17 00:00:00 2001 From: Caleb Morrison Date: Sun, 28 Sep 2025 10:39:05 -0700 Subject: [PATCH 14/17] Add Pchip interp routine for bulletin data interp Pchip interpolation ensures C1 continuity while avoiding global constraints such as in piecewise cubic spline interpolation. The Pchip routine added by this commit is now used by the BulletinData class. The added curvature did increase the error in some tests when compared to the basic finite difference calculation, so the test thresholds had to be increased. --- Tests/test_interpolation.py | 26 +++- Tests/test_itrs_to_gcrs.py | 4 +- Tests/test_polar_motion.py | 2 +- src/TerraFrame/Utilities/BulletinData.py | 10 +- src/TerraFrame/Utilities/Interpolation.py | 142 +++++++++++++++++++++- 5 files changed, 171 insertions(+), 13 deletions(-) diff --git a/Tests/test_interpolation.py b/Tests/test_interpolation.py index a7a5dda..06e19bd 100644 --- a/Tests/test_interpolation.py +++ b/Tests/test_interpolation.py @@ -2,11 +2,12 @@ # License, v. 2.0. If a copy of the MPL was not distributed with this # file, You can obtain one at https://mozilla.org/MPL/2.0/. -from TerraFrame.Utilities import Interpolation import numpy as np +from TerraFrame.Utilities import Interpolation + -def test_interpolation(): +def test_linear_interpolation(): x = np.array([1, 2, 3]) y = np.array([1, 4, 9]) @@ -19,3 +20,24 @@ def test_interpolation(): for i, v in enumerate(yv): assert abs(v - y_answer[i]) < 1e-10 + + +def test_pchip_interpolation(): + x = np.array([-1, 1, 2, 3, 4]) + y = np.array([1, 1, 4, 9, 16]) + + f = Interpolation.InterpolationPchip(x, y) + + xv = np.linspace(-1.5, 4.5, 10) + yv = f(xv) + yvp = f(xv, derivative=True) + + # Answers calculated using Scipy's Pchip routine. + y_answer = [1.0, 1.0, 1.0, 1.0, 1.1354166666666663, 3.3437499999999982, + 6.239583333333333, 10.193672839506169, 15.616512345679013, 16.0] + yp_answer = [0.0, 0.0, 0.0, 0.0, 1.5624999999999987, 4.0625, + 5.104166666666666, 8.263888888888884, 4.3750000000000036, 0.0] + + for i, (v, vp) in enumerate(zip(yv, yvp)): + assert (abs(v - y_answer[i]) < 1e-10) + assert (abs(vp - yp_answer[i]) < 1e-10) diff --git a/Tests/test_itrs_to_gcrs.py b/Tests/test_itrs_to_gcrs.py index 08432fb..7caabf8 100644 --- a/Tests/test_itrs_to_gcrs.py +++ b/Tests/test_itrs_to_gcrs.py @@ -74,7 +74,7 @@ def test_against_astropy(): itrs_basis_astro = gcrs_basis_astro.transform_to(ITRS(obstime=t)) itrs_basis_astro = itrs_basis_astro.cartesian.xyz.to_value() - assert (np.max(np.abs(itrs_basis - itrs_basis_astro)) < 1e-10) + assert (np.max(np.abs(itrs_basis - itrs_basis_astro)) < 1e-9) def test_leap_second_against_astropy(): @@ -136,4 +136,4 @@ def test_t_gi_angular_velocity(): error = np.abs(d_t_gi_dt - d_t_gi_dt_fd) - assert (np.max(error) < 1e-11) + assert (np.max(error) < 1e-9) diff --git a/Tests/test_polar_motion.py b/Tests/test_polar_motion.py index af8d552..73af6d8 100644 --- a/Tests/test_polar_motion.py +++ b/Tests/test_polar_motion.py @@ -99,4 +99,4 @@ def test_itrs_to_tirs_derivative_calculation(): error = np.abs(dt_ti_dt - dt_ti_dt_fd) - assert (np.max(error) < 1e-14) + assert (np.max(error) < 1e-9) diff --git a/src/TerraFrame/Utilities/BulletinData.py b/src/TerraFrame/Utilities/BulletinData.py index 7f46c9f..07c9a2c 100644 --- a/src/TerraFrame/Utilities/BulletinData.py +++ b/src/TerraFrame/Utilities/BulletinData.py @@ -9,7 +9,7 @@ import numpy.typing as npt from .DataFiles import DataFileManager -from .Interpolation import Interpolation1D +from .Interpolation import InterpolationPchip class BulletinData: @@ -84,16 +84,16 @@ def ut1_utc_delta(index): def _init_interpolants(self): if BulletinData.f_pm_x is None: - BulletinData.f_pm_x = Interpolation1D(self.data[:, 0], + BulletinData.f_pm_x = InterpolationPchip(self.data[:, 0], self.data[:, 2]) - BulletinData.f_pm_y = Interpolation1D(self.data[:, 0], + BulletinData.f_pm_y = InterpolationPchip(self.data[:, 0], self.data[:, 3]) - BulletinData.f_nc_dx = Interpolation1D(self.data[:, 0], + BulletinData.f_nc_dx = InterpolationPchip(self.data[:, 0], self.data[:, 4]) - BulletinData.f_nc_dy = Interpolation1D(self.data[:, 0], + BulletinData.f_nc_dy = InterpolationPchip(self.data[:, 0], self.data[:, 5]) def _parse_file(self): diff --git a/src/TerraFrame/Utilities/Interpolation.py b/src/TerraFrame/Utilities/Interpolation.py index 2d21a3e..79e171f 100644 --- a/src/TerraFrame/Utilities/Interpolation.py +++ b/src/TerraFrame/Utilities/Interpolation.py @@ -3,6 +3,7 @@ # file, You can obtain one at https://mozilla.org/MPL/2.0/. import numpy as np + from .Helpers import clamp, ensure_iterable @@ -33,7 +34,6 @@ def _get_index(self, xv): if x1 < xv <= x2: index = self._index_cache - if index is None: index = np.searchsorted(self._x, xv) self._index_cache = index @@ -50,10 +50,16 @@ def __call__(self, xv, derivative=False): # If we're out of bounds, return the first or last value if index == 0: - yv[i] = self._y[index] + if derivative: + yv[i] = 0.0 + else: + yv[i] = self._y[index] continue elif index >= len(self._x): - yv[i] = self._y[-1] + if derivative: + yv[i] = 0.0 + else: + yv[i] = self._y[-1] continue y2 = self._y[index] @@ -73,3 +79,133 @@ def __call__(self, xv, derivative=False): return yv[0] else: return yv + + +class InterpolationPchip: + """ + This is not a general PCHIP interpolation class. It is specialized to the + specific usage requirements and patterns of TerraFrame. + + This Pchip routine is based on the method laid out by Butland et al. in [1]. + + Since most queries will be near each other, there is an index cache. If + the index is out of bounds, the first or last values are used. + + References: + [1] F., N., And, F., J., & Butland (1984). A Method for Constructing + Local Monotone Piecewise Cubic Interpolants. Siam Journal on + Scientific and Statistical Computing, 5, 300-304. + DOI: https://doi.org/10.1137/0905021 + """ + + def __init__(self, x, y): + self._x = x + self._y = y + + self._index_cache = None + self._coefficient_cache = None + + def _get_index(self, xv): + index = None + + # Under nominal usage patterns, most queries will use the same index + # with only the occasional change + if self._index_cache is not None: + x1 = self._x[clamp(self._index_cache - 1, 0, len(self._x) - 1)] + x2 = self._x[clamp(self._index_cache, 0, len(self._x) - 1)] + + if x1 < xv <= x2: + index = self._index_cache + + if index is None: + # If the index cache isn't correct, we need to clear the + # coefficient cache. + self._coefficient_cache = None + index = np.searchsorted(self._x, xv) + self._index_cache = index + + return index + + def _get_gi(self, index): + if index == 0 or index >= len(self._x) - 1: + return 0 + + h2 = self._x[index + 1] - self._x[index + 0] + h1 = self._x[index + 0] - self._x[index - 1] + + s2 = (self._y[index + 1] - self._y[index + 0]) / h2 + s1 = (self._y[index + 0] - self._y[index - 1]) / h1 + + if s1 * s2 > 0: + alpha = 1.0 / 3.0 * (1 + h2 / (h1 + h2)) + + gi = s1 * s2 / (alpha * s2 + (1 - alpha) * s1) + else: + gi = 0 + + return gi + + def __call__(self, xv, derivative=False): + xv = ensure_iterable(xv) + + yv = len(xv) * [0.0] + + for i, v in enumerate(xv): + index = self._get_index(v) + + # If we're out of bounds, return the first or last value + if index == 0: + if derivative: + yv[i] = 0.0 + else: + yv[i] = self._y[index] + continue + elif index >= len(self._x): + if derivative: + yv[i] = 0.0 + else: + yv[i] = self._y[-1] + continue + + if self._coefficient_cache is None: + g1 = self._get_gi(index - 1) + g2 = self._get_gi(index + 0) + + y2 = self._y[index] + y1 = self._y[index - 1] + + x2 = self._x[index] + x1 = self._x[index - 1] + + a = np.array(( + (1, x1, x1 ** 2, x1 ** 3), + (1, x2, x2 ** 2, x2 ** 3), + (0, 1, 2 * x1, 3 * x1 ** 2), + (0, 1, 2 * x2, 3 * x2 ** 2) + )) + + b = np.array(( + (y1, ), + (y2, ), + (g1, ), + (g2, ) + )) + + coefficients = np.linalg.solve(a, b) + + self._coefficient_cache = coefficients.flatten() + + if not derivative: + yv[i] = (self._coefficient_cache[0] + + self._coefficient_cache[1] * xv[i] + + self._coefficient_cache[2] * xv[i] ** 2 + + self._coefficient_cache[3] * xv[i] ** 3) + else: + yv[i] = (self._coefficient_cache[1] + + self._coefficient_cache[2] * 2 * xv[i] + + self._coefficient_cache[3] * 3 * xv[i] ** 2) + + if len(yv) == 1: + return yv[0] + else: + return yv From ede227ddbe27b78cd72f4e53f726018b621f4fe4 Mon Sep 17 00:00:00 2001 From: Caleb Morrison Date: Sun, 28 Sep 2025 10:53:45 -0700 Subject: [PATCH 15/17] Adjust test assert threshold --- Tests/test_itrs_to_gcrs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Tests/test_itrs_to_gcrs.py b/Tests/test_itrs_to_gcrs.py index 7caabf8..15d7a66 100644 --- a/Tests/test_itrs_to_gcrs.py +++ b/Tests/test_itrs_to_gcrs.py @@ -136,4 +136,4 @@ def test_t_gi_angular_velocity(): error = np.abs(d_t_gi_dt - d_t_gi_dt_fd) - assert (np.max(error) < 1e-9) + assert np.max(error) < 1e-8, f"val: {val}" From 1fbabeed2844df6629c6d4caadbb52a9f2f23d5b Mon Sep 17 00:00:00 2001 From: Caleb Morrison Date: Sun, 28 Sep 2025 10:58:10 -0700 Subject: [PATCH 16/17] Add inverse angular vel func to transform class --- src/TerraFrame/CelestialTerrestrial.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/TerraFrame/CelestialTerrestrial.py b/src/TerraFrame/CelestialTerrestrial.py index 790f0f4..279fbb7 100644 --- a/src/TerraFrame/CelestialTerrestrial.py +++ b/src/TerraFrame/CelestialTerrestrial.py @@ -30,7 +30,7 @@ def __init__(self, use_polar_motion=True, use_nutation_corrections=True): self.t_gc = None self.t_ct = None self.t_ti = None - self.angular_velocity = None + self.w_gi = None def itrs_to_gcrs(self, time): if isinstance(time, datetime.datetime): @@ -237,6 +237,11 @@ def itrs_to_gcrs_angular_vel(self, time): self.t_gc = t_gc self.t_ct = t_ct self.t_ti = t_ti - self.angular_velocity = angular_vel + self.w_gi = angular_vel return t_gi, angular_vel + + def gcrs_to_itrs_angular_vel(self, time): + t_gi, w_gi = self.itrs_to_gcrs_angular_vel(time) + + return t_gi.T, w_gi.T From e7276ae18d759eb23360e47ccf756a7e90930aca Mon Sep 17 00:00:00 2001 From: Caleb Morrison Date: Sun, 28 Sep 2025 11:02:38 -0700 Subject: [PATCH 17/17] Update readme with angular velocity info --- README.md | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index 5987e04..f74f78b 100644 --- a/README.md +++ b/README.md @@ -6,6 +6,10 @@ TerraFrame is a library designed to provided key Earth related calculations and associated functionality to modeling & simulation software. Primarily, TerraFrame provides Earth orientation and gravity routines. +![Animation of CGRS to ITRS Transformation](https://raw.githubusercontent.com/cmorrison31/TerraFrame/main/Animations%20and%20Plots/GCRS_to_ITRS.gif) + +![Example of Precession, Nutation, & Polar Motion](https://raw.githubusercontent.com/cmorrison31/TerraFrame/main/Animations%20and%20Plots/Earth%20Motion%20Example.gif) + ## Earth Orientation TerraFrame provides an implementation of the IAU 2006/2000A precession-nutation model which accounts for precession, nutation, and polar motion. Specifically, @@ -13,6 +17,11 @@ this implementation provides a transformation tensor between the Geocentric Celestial Reference System (GCRS) and the International Terrestrial Reference System (ITRS). +TerraFrame also provides routines for calculating the Earth's angular +velocity tensor. The angular velocity tensor is based on the partial +derivative of the numerous IAU 2006/2000A equations with respect to time (in +seconds). + IERS precession and nutation model data files are shipped with TerraFrame. Utility code is also provided which automates the downloading of IERS data for polar motion, UTC, and UT1 offsets. @@ -39,10 +48,6 @@ provided. The user is encouraged to not work in UTC directly to avoid leap second ambiguity. Conversion to UTC from TT or TAI can be safely done in post-processing. -![Animation of CGRS to ITRS Transformation](https://raw.githubusercontent.com/cmorrison31/TerraFrame/main/Animations%20and%20Plots/GCRS_to_ITRS.gif) - -![Example of Precession, Nutation, & Polar Motion](https://raw.githubusercontent.com/cmorrison31/TerraFrame/main/Animations%20and%20Plots/Earth%20Motion%20Example.gif) - # License This project - except for the IERS and WGS84 data files - is covered under the