Fix t_cg derivative (e.g., Q'(t)) calc & add test

Author: cmorrison31

Tests/test_cgrs_to_cirs.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):

Tests/test_cip_calculation.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

Tests/test_cirs_to_gcrs.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

src/TerraFrame/CelestialTerrestrial.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