diff --git a/README.md b/README.md
index 8b0c752..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,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
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 e589830..8c09567 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,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)
diff --git a/Tests/test_cirs_to_gcrs.py b/Tests/test_cirs_to_gcrs.py
new file mode 100644
index 0000000..1853119
--- /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/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)
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 8ef8667..15d7a66 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)
 
@@ -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}"
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
new file mode 100644
index 0000000..c9d0266
--- /dev/null
+++ b/Tests/test_matrix_derivatives.py
@@ -0,0 +1,73 @@
+# 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
+
+
+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)
+    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] - dpsi)
+        r2 = TransformationMatrices.r3(psi[i] + dpsi)
+
+        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-6)
diff --git a/Tests/test_polar_motion.py b/Tests/test_polar_motion.py
index fef7f29..73af6d8 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,56 @@ 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))
+    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)
+    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)
+    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)
+    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-9)
diff --git a/src/TerraFrame/CelestialTerrestrial.py b/src/TerraFrame/CelestialTerrestrial.py
index 803e582..279fbb7 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()
@@ -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.w_gi = None
 
     def itrs_to_gcrs(self, time):
         if isinstance(time, datetime.datetime):
@@ -67,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))
 
@@ -88,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:
@@ -116,3 +117,131 @@ 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._use_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)
+            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
+
+        # 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))
+
+        # 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.
+        # R(t)
+        t_ct = TransformationMatrices.earth_rotation_matrix(jd_ut1)
+
+        # R'(t)
+        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._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
+            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()
+
+        # 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
+
+        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
+        self.t_ct = t_ct
+        self.t_ti = t_ti
+        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
diff --git a/src/TerraFrame/Earth.py b/src/TerraFrame/Earth.py
index f10e7eb..5370b63 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
 
 
@@ -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
@@ -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/Arguments.py b/src/TerraFrame/PrecessionNutation/Arguments.py
index 570c603..2d8d4a8 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,33 @@ 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)
+
+    return value
+
+
 def mean_anomaly_of_the_sun(time):
     """
     This function computes the mean anomaly of the Sun per IERS Conventions
@@ -51,7 +79,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 +95,33 @@ 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)
+
+    return value
+
+
 def mean_longitude_moon_minus_ascending_node(time):
     """
     This function computes Mean longitude of the Moon minus the ascending node
@@ -80,7 +136,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 +152,33 @@ 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)
+
+    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 +193,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 +209,33 @@ 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)
+
+    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 +250,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 +266,33 @@ 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)
+
+    return value
+
+
 def mean_longitude_of_mercury(time):
     """
     This function computes the mean longitude of Mercury per IERS Conventions
@@ -167,18 +307,40 @@ 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
     value = 4.402608842 + 2608.7903141574 * time
 
-    # Take the modulus before convertion to maintain accuracy
     value = math.fmod(value, 2.0 * math.pi)
 
     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
+
+    return value
+
+
 def mean_longitude_of_venus(time):
     """
     This function computes the mean longitude of Venus per IERS Conventions
@@ -193,7 +355,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 +368,28 @@ 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
+
+    return value
+
+
 def mean_longitude_of_earth(time):
     """
     This function computes the mean longitude of Earth per IERS Conventions
@@ -219,7 +404,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 +417,28 @@ 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
+
+    return value
+
+
 def mean_longitude_of_mars(time):
     """
     This function computes the mean longitude of Mars per IERS Conventions
@@ -245,7 +453,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 +466,28 @@ 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
+
+    return value
+
+
 def mean_longitude_of_jupiter(time):
     """
     This function computes the mean longitude of Jupiter per IERS Conventions
@@ -271,7 +502,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 +515,28 @@ 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
+
+    return value
+
+
 def mean_longitude_of_saturn(time):
     """
     This function computes the mean longitude of Saturn per IERS Conventions
@@ -297,7 +551,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 +564,28 @@ 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
+
+    return value
+
+
 def mean_longitude_of_uranus(time):
     """
     This function computes the mean longitude of Uranus per IERS Conventions
@@ -323,7 +600,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 +613,28 @@ 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
+
+    return value
+
+
 def mean_longitude_of_neptune(time):
     """
     This function computes the mean longitude of Neptune per IERS Conventions
@@ -349,7 +649,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 +662,28 @@ 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
+
+    return value
+
+
 def general_precession_in_longitude(time):
     """
     This function computes the general precession in longitude per IERS
@@ -375,10 +698,35 @@ 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
+
+    return value
diff --git a/src/TerraFrame/PrecessionNutation/SeriesExpansion.py b/src/TerraFrame/PrecessionNutation/SeriesExpansion.py
index 7b0a62d..72ff18a 100644
--- a/src/TerraFrame/PrecessionNutation/SeriesExpansion.py
+++ b/src/TerraFrame/PrecessionNutation/SeriesExpansion.py
@@ -52,19 +52,39 @@ 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):
+        """
+        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
         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)
+
         # 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'
@@ -82,7 +102,37 @@ def compute(self, t):
         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)
+            )
+
         non_poly_part = 0.0
+        d_non_poly_part_dt = 0.0
 
         for row in self.data:
             j = row[0]
@@ -91,13 +141,37 @@ def compute(self, t):
 
             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)
+
+        # 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)
 
-        return total
+        if derivative:
+            return total, total_dt
+        else:
+            return total
 
 
 def cip_x(file_name=r'tab5.2a.txt'):
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/Conversions.py b/src/TerraFrame/Utilities/Conversions.py
index a811f2a..928051c 100644
--- a/src/TerraFrame/Utilities/Conversions.py
+++ b/src/TerraFrame/Utilities/Conversions.py
@@ -295,3 +295,33 @@ 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)
+
+
+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
diff --git a/src/TerraFrame/Utilities/Interpolation.py b/src/TerraFrame/Utilities/Interpolation.py
index 822f988..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,14 +34,13 @@ 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
 
         return index
 
-    def __call__(self, xv):
+    def __call__(self, xv, derivative=False):
         xv = ensure_iterable(xv)
 
         yv = len(xv) * [0.0]
@@ -50,10 +50,16 @@ def __call__(self, xv):
 
             # 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]
@@ -62,7 +68,142 @@ 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]
+        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]
diff --git a/src/TerraFrame/Utilities/TransformationMatrices.py b/src/TerraFrame/Utilities/TransformationMatrices.py
index 470189f..bd55e5b 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,12 +149,45 @@ 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
 
 
+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
@@ -197,6 +295,30 @@ 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.seconds_to_centuries(-47e-6)
+    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
@@ -229,6 +351,53 @@ 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) @ 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
+
+
 def earth_rotation_matrix(time):
     """
     This function computes the earth rotation matrix at a given datetime in UT1.
@@ -248,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
@@ -271,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