feature/angular_velocity

README.md

@@ -6,13 +6,22 @@ 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, 
 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,14 +48,10 @@ 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 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

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,55 @@ 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 = 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)

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)

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)

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)

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)
 
@@ -72,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():
@@ -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-8, f"val: {val}"