From 723f357716b3fb9dd821f5c5ceeb831db32e9662 Mon Sep 17 00:00:00 2001
From: MicheleCeresoli <michele.ceresoli@polimi.it>
Date: Sat, 10 Feb 2024 15:38:00 +0100
Subject: [PATCH 1/5] Fixed `mod` and `rem` operations

---
 src/delaunay.jl  | 86 +++++++++++++++++++++++++++++++-----------------
 src/era.jl       |  6 ++--
 src/planetary.jl | 30 ++++++++---------
 src/rate.jl      | 11 ++-----
 src/rotations.jl |  9 +++++
 src/sidereal.jl  | 11 ++++---
 6 files changed, 91 insertions(+), 62 deletions(-)

diff --git a/src/delaunay.jl b/src/delaunay.jl
index 097d862..3a4ca60 100644
--- a/src/delaunay.jl
+++ b/src/delaunay.jl
@@ -49,20 +49,24 @@ centuries since J2000.
 """
 function delaunay_anomaly_moon(::IERSAModels, tdb_c::Number)
     # This expression is valid for both IERS2003 and IERS2010
-    mod2pi(
-        arcsec2rad(
-            @evalpoly(tdb_c, 485868.249036, 1717915923.2178, 31.8792, 0.051635, -0.00024470)
+    arcsec2rad(
+        jrem(
+            @evalpoly(tdb_c, 485868.249036, 1717915923.2178, 31.8792, 0.051635, -0.00024470), 
+            1296000
         )
     )
 end
 
 function delaunay_anomaly_moon(::IERSModel, tdb_c::Number)
-    arcsec2rad(mod(485868.249036 + 1717915923.2178tdb_c, 1296000))
+    arcsec2rad(jrem(485868.249036 + 1717915923.2178tdb_c, 1296000))
 end
 
 function delaunay_anomaly_moon(::IERS1996, tdb_c::Number)
-    r = mod(1325tdb_c, 1)
-    mod2pi(arcsec2rad(@evalpoly(tdb_c, 485866.733, 715922.633, 31.310, 0.064)) + 2π*r)
+    r = jrem(1325tdb_c, 1)
+    rem2pi(
+        arcsec2rad(@evalpoly(tdb_c, 485866.733, 715922.633, 31.310, 0.064)) + 2π*r, 
+        RoundNearest
+    )
 end
 
 
@@ -73,20 +77,24 @@ Compute the mean anomaly of the Sun, in radians, given time `tdb_c` expressed in
 centuries since J2000.
 """
 function delaunay_anomaly_sun(::IERSAModels, tdb_c::Number)
-    mod2pi(
-        arcsec2rad(
-            @evalpoly(tdb_c, 1287104.793048, 129596581.0481, -0.5532, 0.000136, -0.00001149)
+    arcsec2rad(
+        jrem(
+            @evalpoly(tdb_c, 1287104.793048, 129596581.0481, -0.5532, 0.000136, -0.00001149), 
+            1296000
         )
     )
 end
 
 function delaunay_anomaly_sun(::IERSModel, tdb_c::Number)
-    arcsec2rad(mod(1287104.79305 + 129596581.0481tdb_c, 1296000))
+    arcsec2rad(jrem(1287104.79305 + 129596581.0481tdb_c, 1296000))
 end
 
 function delaunay_anomaly_sun(::IERS1996, tdb_c::Number)
-    r = mod(99tdb_c, 1)
-    mod2pi(arcsec2rad(@evalpoly(tdb_c, 1287099.804, 1292581.224, -0.577, -0.012)) + 2π*r)
+    r = jrem(99tdb_c, 1)
+    rem2pi(
+        arcsec2rad(@evalpoly(tdb_c, 1287099.804, 1292581.224, -0.577, -0.012)) + 2π*r, 
+        RoundNearest
+    )
 end
 
 
@@ -97,20 +105,24 @@ Compute the difference between the longitude of the Moon and the longitude of th
 node, in radians, given time `tdb_c` expressed in TDB Julian centuries since J2000.
 """
 function delaunay_longitude_diff(::IERSAModels, tdb_c::Number)
-    mod2pi(
-        arcsec2rad(
-            @evalpoly(tdb_c, 335779.526232, 1739527262.8478, -12.7512, -0.001037, +0.00000417)
+    arcsec2rad(
+        jrem(
+            @evalpoly(tdb_c, 335779.526232, 1739527262.8478, -12.7512, -0.001037, 0.00000417), 
+            1296000
         )
     )
 end
 
 function delaunay_longitude_diff(::IERSModel, tdb_c::Number)
-    arcsec2rad(mod(335779.526232 + 1739527262.8478tdb_c, 1296000))
+    arcsec2rad(jrem(335779.526232 + 1739527262.8478tdb_c, 1296000))
 end
 
 function delaunay_longitude_diff(::IERS1996, tdb_c::Number)
-    r = mod(1342tdb_c, 1)
-    mod2pi(arcsec2rad(@evalpoly(tdb_c, 335778.877, 295263.137, -13.257, 0.011)) + 2π*r)
+    r = jrem(1342tdb_c, 1)
+    rem2pi(
+        arcsec2rad(@evalpoly(tdb_c, 335778.877, 295263.137, -13.257, 0.011)) + 2π*r, 
+        RoundNearest
+    )
 end
 
 """
@@ -120,20 +132,24 @@ Compute the mean elongation of the Moon from the Sun, in radians, given time `td
 TDB Julian centuries since J2000.
 """
 function delaunay_elongation_moon(::IERSAModels, tdb_c::Number)
-    mod2pi(
-        arcsec2rad(
-            @evalpoly(tdb_c, 1072260.703692, 1602961601.2090, -6.3706, 0.006593, -0.00003169)
+    arcsec2rad(
+        jrem(
+            @evalpoly(tdb_c, 1072260.703692, 1602961601.2090, -6.3706, 0.006593, -0.00003169), 
+            1296000
         )
     )
 end
 
 function delaunay_elongation_moon(::IERSModel, tdb_c::Number)
-    arcsec2rad(mod(1072260.70369 + 1602961601.2090tdb_c, 1296000))
+    arcsec2rad(jrem(1072260.70369 + 1602961601.2090tdb_c, 1296000))
 end
 
 function delaunay_elongation_moon(::IERS1996, tdb_c::Number)
-    r = mod(1236tdb_c, 1)
-    mod2pi(arcsec2rad(@evalpoly(tdb_c, 1072261.307, 1105601.328, -6.891, 0.019)) + 2π*r)
+    r = jrem(1236tdb_c, 1)
+    rem2pi(
+        arcsec2rad(@evalpoly(tdb_c, 1072261.307, 1105601.328, -6.891, 0.019)) + 2π*r, 
+        RoundNearest
+    )
 end
 
 """
@@ -144,18 +160,26 @@ measured from the mean equinox of date, in radians, given time `tdb_c` expressed
 centuries since J2000.
 """
 function delaunay_longitude_node(::IERSAModels, tdb_c::Number)
-    mod2pi(
-        arcsec2rad(
-            @evalpoly(tdb_c, 450160.398036, -6962890.5431, 7.4722, 0.007702, -0.00005939)
+    arcsec2rad(
+        jrem(
+            @evalpoly(tdb_c, 450160.398036, -6962890.5431, 7.4722, 0.007702, -0.00005939), 
+            1296000
         )
     )
 end
 
 function delaunay_longitude_node(::IERSModel, tdb_c::Number)
-    arcsec2rad(mod(450160.398036 - 6962890.5431tdb_c, 1296000))
+    arcsec2rad(jrem(450160.398036 - 6962890.5431tdb_c, 1296000))
 end
 
 function delaunay_longitude_node(::IERS1996, tdb_c::Number)
-    r = mod(-5tdb_c, 1)
-    mod2pi(arcsec2rad(@evalpoly(tdb_c, 450160.280, -482890.539, 7.455, 0.008)) + 2π*r)
-end
\ No newline at end of file
+    r = jrem(-5tdb_c, 1)
+    rem2pi(
+        arcsec2rad(@evalpoly(tdb_c, 450160.280, -482890.539, 7.455, 0.008)) + 2π*r, 
+        RoundNearest
+    )
+end
+
+
+jrem(x, y) = x - y*(x ÷ y)
+jrem2pi(x) = jrem(x, 2π)
\ No newline at end of file
diff --git a/src/era.jl b/src/era.jl
index 4f3b4cc..cf0f8b5 100644
--- a/src/era.jl
+++ b/src/era.jl
@@ -18,10 +18,12 @@ since `J2000`, according to the IERS convention `m`.
 See also [`iers_era_rotm`](@ref).
 """
 function iers_era(::IERSModel, ut1_d::Number)
-
     # The function uses the fractional UT1 date to gain additional precision
     # in the computations
-    return mod2pi(2π * (mod(ut1_d, 1) + 0.7790572732640 + 0.00273781191135448ut1_d))
+
+    # Uses integer divsion because the `mod` and `rem` functions don't work with ForwardDiff
+    f = ut1_d - (ut1_d ÷ 1)
+    return rem2pi(2π * (f + 0.7790572732640 + 0.00273781191135448ut1_d), RoundDown)
 
 end
 
diff --git a/src/planetary.jl b/src/planetary.jl
index dd8923e..01b2ccc 100644
--- a/src/planetary.jl
+++ b/src/planetary.jl
@@ -68,7 +68,7 @@ Julian centuries since `J2000`.
 - [ERFA](https://github.com/liberfa/erfa/blob/master/src/fame03.c) software library
 """
 function pa_mercury(::IERSModel, tdb_c::Number)
-    return mod2pi(@evalpoly(tdb_c, 4.402608842, 2608.7903141574))
+    return rem2pi(@evalpoly(tdb_c, 4.402608842, 2608.7903141574), RoundNearest)
 end
 
 function pa_mercury(::IERS1996, ::Number)
@@ -93,11 +93,11 @@ Julian centuries since `J2000`.
 - [ERFA](https://github.com/liberfa/erfa/blob/master/src/fave03.c) software library
 """
 function pa_venus(::IERSModel, tdb_c::Number)
-    return mod2pi(@evalpoly(tdb_c, 3.176146697, 1021.3285546211))
+    return rem2pi(@evalpoly(tdb_c, 3.176146697, 1021.3285546211), RoundNearest)
 end
 
 function pa_venus(::IERS1996, tdb_c::Number)
-    return mod2pi(deg2rad(@evalpoly(tdb_c, 181.979800853, 58517.8156748)))
+    return rem2pi(deg2rad(@evalpoly(tdb_c, 181.979800853, 58517.8156748)), RoundNearest)
 end
 
 
@@ -113,11 +113,11 @@ Julian centuries since `J2000`.
 - [ERFA](https://github.com/liberfa/erfa/blob/master/src/fae03.c) software library
 """
 function pa_earth(::IERSModel, tdb_c::Number)
-    return mod2pi(@evalpoly(tdb_c, 1.753470314, 628.3075849991))
+    return rem2pi(@evalpoly(tdb_c, 1.753470314, 628.3075849991), RoundNearest)
 end
 
 function pa_earth(::IERS1996, tdb_c::Number)
-    return mod2pi(deg2rad(@evalpoly(tdb_c, 100.466448494, 35999.3728521)))
+    return rem2pi(deg2rad(@evalpoly(tdb_c, 100.466448494, 35999.3728521)), RoundNearest)
 end
 
 
@@ -133,11 +133,11 @@ Julian centuries since `J2000`.
 - [ERFA](https://github.com/liberfa/erfa/blob/master/src/fama03.c) software library
 """
 function pa_mars(::IERSModel, tdb_c::Number)
-    return mod2pi(@evalpoly(tdb_c, 6.203480913, 334.0612426700))
+    return rem2pi(@evalpoly(tdb_c, 6.203480913, 334.0612426700), RoundNearest)
 end
 
 function pa_mars(::IERS1996, tdb_c::Number)
-    return mod2pi(deg2rad(@evalpoly(tdb_c, 355.433274605, 19140.299314)))
+    return rem2pi(deg2rad(@evalpoly(tdb_c, 355.433274605, 19140.299314)), RoundNearest)
 end
 
 
@@ -153,11 +153,11 @@ Julian centuries since `J2000`.
 - [ERFA](https://github.com/liberfa/erfa/blob/master/src/faju03.c) software library
 """
 function pa_jupiter(::IERSModel, tdb_c::Number)
-    return mod2pi(@evalpoly(tdb_c, 0.599546497, 52.9690962641))
+    return rem2pi(@evalpoly(tdb_c, 0.599546497, 52.9690962641), RoundNearest)
 end
 
 function pa_jupiter(::IERS1996, tdb_c::Number)
-    return mod2pi(deg2rad(@evalpoly(tdb_c, 34.3514839, 3034.90567464)))
+    return rem2pi(deg2rad(@evalpoly(tdb_c, 34.3514839, 3034.90567464)), RoundNearest)
 end
 
 
@@ -173,11 +173,11 @@ Julian centuries since `J2000`.
 - [ERFA](https://github.com/liberfa/erfa/blob/master/src/fasa03.c) software library
 """
 function pa_saturn(::IERSModel, tdb_c::Number)
-    return mod2pi(@evalpoly(tdb_c, 0.874016757, 21.3299104960))
+    return rem2pi(@evalpoly(tdb_c, 0.874016757, 21.3299104960), RoundNearest)
 end
 
 function pa_saturn(::IERS1996, tdb_c::Number)
-    return mod2pi(deg2rad(@evalpoly(tdb_c, 50.0774713998, 1222.11379404)))
+    return rem2pi(deg2rad(@evalpoly(tdb_c, 50.0774713998, 1222.11379404)), RoundNearest)
 end
 
 
@@ -196,7 +196,7 @@ Julian centuries since `J2000`.
 - [ERFA](https://github.com/liberfa/erfa/blob/master/src/faur03.c) software library
 """
 function pa_uranus(::IERSModel, tdb_c::Number)
-    return mod2pi(@evalpoly(tdb_c, 5.481293872, 7.4781598567))
+    return rem2pi(@evalpoly(tdb_c, 5.481293872, 7.4781598567), RoundNearest)
 end
 
 function pa_uranus(::IERS1996, ::Number)
@@ -224,7 +224,7 @@ Julian centuries since `J2000`.
 - [ERFA](https://github.com/liberfa/erfa/blob/master/src/fane03.c) software library
 """
 function pa_neptune(::IERSModel, tdb_c::Number)
-    return mod2pi(@evalpoly(tdb_c, 5.311886287, 3.8133035638))
+    return rem2pi(@evalpoly(tdb_c, 5.311886287, 3.8133035638), RoundNearest)
 end
 
 function pa_neptune(::IERS1996, ::Number)
@@ -249,9 +249,9 @@ expressed in `TDB` Julian centuries since `J2000`.
 - [ERFA](https://github.com/liberfa/erfa/blob/master/src/fapa03.c) software library
 """
 function pa_precession(::IERSModel, tdb_c::Number)
-    return mod2pi(@evalpoly(tdb_c, 0, 0.024381750, 0.00000538691))
+    return jrem2pi(@evalpoly(tdb_c, 0, 0.024381750, 0.00000538691))
 end
 
 function pa_precession(::IERS1996, tdb_c::Number)
-    return mod2pi(deg2rad(@evalpoly(tdb_c, 0, 1.39697137214, 0.0003086)))
+    return jrem2pi(deg2rad(@evalpoly(tdb_c, 0, 1.39697137214, 0.0003086)))
 end
diff --git a/src/rate.jl b/src/rate.jl
index 4eb73d0..810db90 100644
--- a/src/rate.jl
+++ b/src/rate.jl
@@ -4,18 +4,11 @@ export iers_earth_rot_rate
 const ωₑ = 7.292_115_146_706_979e-5
 
 """ 
-    iers_earth_rot_rate(LOD::Number)
+    iers_earth_rot_rate(LOD::Number=0)
 
 Compute the true angular velocity of the Earth accounting for the Length of the Day, i.e., 
 the instantaneous rate of change of UT1 with respect to a uniform time scale. 
 """
-function iers_earth_rot_rate(LOD::Number)
+function iers_earth_rot_rate(LOD::Number=0)
     return ωₑ * (1 - LOD / 86400)
 end
-
-"""
-    iers_earth_rot_rate()
-    
-Compute the nominal Earth angular velocity. 
-"""
-iers_earth_rot_rate() = iers_earth_rot_rate(0.0)
\ No newline at end of file
diff --git a/src/rotations.jl b/src/rotations.jl
index 9f688e9..fb7c33f 100644
--- a/src/rotations.jl
+++ b/src/rotations.jl
@@ -376,6 +376,9 @@ in TT seconds since `J2000`, following the IERS Conventions `m`.
 !!! note 
     Polar motion is neglected in the [`CPNd`](@ref) model. 
 
+!!! note 
+    The time derivative of the nutation and precession effects is neglected.
+
 ### References 
 - IERS Technical Note No. [36](https://www.iers.org/IERS/EN/Publications/TechnicalNotes/tn36.html) 
 """
@@ -414,6 +417,9 @@ Reference Frame (GCRF) to the International Terrestrial Reference Frame (ITRF) a
 !!! note 
     Polar motion is neglected in the [`CPNd`](@ref) model. 
 
+!!! note 
+    The time derivative of the nutation and precession effects is neglected.
+
 ### References 
 - IERS Technical Note No. [36](https://www.iers.org/IERS/EN/Publications/TechnicalNotes/tn36.html) 
 """
@@ -453,6 +459,9 @@ at time `tt_s`, expressed in TT seconds since `J2000`, following the IERS Conven
 !!! note 
     Polar motion is neglected in the [`CPNd`](@ref) model. 
 
+!!! note 
+    The time derivative of the nutation and precession effects is neglected.
+
 ### References 
 - IERS Technical Note No. [36](https://www.iers.org/IERS/EN/Publications/TechnicalNotes/tn36.html) 
 """
diff --git a/src/sidereal.jl b/src/sidereal.jl
index 265e30e..e2d691e 100644
--- a/src/sidereal.jl
+++ b/src/sidereal.jl
@@ -118,8 +118,9 @@ end
 # t = TT Julian centuries since J2000 
 function iers_gmst(::IERS1996, tt_c::Number)
 
-    # Transform from TT centuries to UT1 centuries
+    # Transform from TT centuries to UT1 centuries and days
     ut1_c = tt_c + offset_tt2ut1(tt_c*Tempo.CENTURY2SEC)/Tempo.CENTURY2SEC
+    ut1_d = ut1_c*Tempo.CENTURY2DAY
 
     # Coefficients for the IAU 1982 GMST-UT1 model. The first component has been adjusted 
     # of 12 hours because the input is defined with respect to noon of 01-01-2000
@@ -129,10 +130,10 @@ function iers_gmst(::IERS1996, tt_c::Number)
     D = -6.2e-6
 
     # Fractional part of UT1, in seconds
-    f = Tempo.DAY2SEC*mod(ut1_c*Tempo.CENTURY2DAY, 1)
+    f = Tempo.DAY2SEC*(ut1_d - (ut1_d ÷ 1))
 
     # Compute GMST
-    return mod2pi(2π/86400*(@evalpoly(ut1_c, A, B, C, D) + f))
+    return rem2pi(2π/86400*(@evalpoly(ut1_c, A, B, C, D) + f), RoundDown)
 
 end
 
@@ -164,7 +165,7 @@ function iers_gmst(::IERS2003, tt_c::Number, θ::Number)
     )
 
     # Compute GMST 
-    return mod2pi(θ + arcsec2rad(p))
+    return rem2pi(θ + arcsec2rad(p), RoundDown)
 
 end
 
@@ -190,6 +191,6 @@ function iers_gmst(::IERS2010, tt_c::Number, θ::Number)
     )
 
     # Compute GMST 
-    return mod2pi(θ + arcsec2rad(p))
+    return rem2pi(θ + arcsec2rad(p), RoundDown)
     
 end
\ No newline at end of file

From 411671cafa212bf14f2d1bdb433c17137211ed47 Mon Sep 17 00:00:00 2001
From: MicheleCeresoli <michele.ceresoli@polimi.it>
Date: Sat, 10 Feb 2024 15:38:10 +0100
Subject: [PATCH 2/5] Updated tests with new remainer ops

---
 test/fa.jl       | 46 ++++++++++++++++++++++------------------------
 test/sidereal.jl |  3 ++-
 2 files changed, 24 insertions(+), 25 deletions(-)

diff --git a/test/fa.jl b/test/fa.jl
index 7078759..20444b3 100644
--- a/test/fa.jl
+++ b/test/fa.jl
@@ -3,8 +3,6 @@ r2a = 180 / π * 3600
 
 @testset "Delaunay Arguments" verbose=true begin
 
-    r2a = 180 / π * 3600
-
     da = [
         IERSConventions.DelaunayArgs(iers1996, 0.06567),
         IERSConventions.DelaunayArgs(iers1996, 1.0)
@@ -12,24 +10,24 @@ r2a = 180 / π * 3600
 
     @testset "1996" begin 
         # -- Testing Mean Anomaly of the Moon (< 1μas)
-        @test r2a*abs(da[1].M - 2.66359306441736)                   ≤ 1e-6
-        @test r2a*abs(da[2].M - mod(-4.56593937247721e-1, 2π))      ≤ 1e-6
+        @test r2a*abs(da[1].M - 2.66359306441736)         ≤ 1e-6
+        @test r2a*abs(da[2].M - -4.56593937247721e-1)     ≤ 1e-6
 
         # -- Testing Mean Anomaly of the Sun (< 1μas)
-        @test r2a*abs(da[1].S - mod(-2.76485707808279, 2π))         ≤ 1e-6
-        @test r2a*abs(da[2].S - mod(-5.97269171806332e-2, 2π))      ≤ 1e-6
+        @test r2a*abs(da[1].S - -2.76485707808279)        ≤ 1e-6
+        @test r2a*abs(da[2].S - -5.97269171806332e-2)     ≤ 1e-6
 
         # -- Testing Mean Argument of Latitude of the Moon (< 1μas)
-        @test r2a*abs(da[1].F - 2.53331724178132)                   ≤ 1e-6
-        @test r2a*abs(da[2].F - 3.05931379900095)                   ≤ 1e-6
+        @test r2a*abs(da[1].F - 2.53331724178132)         ≤ 1e-6
+        @test r2a*abs(da[2].F - 3.05931379900095)         ≤ 1e-6
 
         # -- Testing mean elongation of the moon from the sun (< 1μas)
-        @test r2a*abs(da[1].D - 3.23611369830128e-1)                ≤ 1e-6
-        @test r2a*abs(da[2].D - mod(-2.00782792050270, 2π))         ≤ 1e-6
+        @test r2a*abs(da[1].D - 3.23611369830128e-1)      ≤ 1e-6
+        @test r2a*abs(da[2].D - -2.00782792050270)        ≤ 1e-6
 
         # -- Testing Mean Longitude of the Moon (< 1μas)
-        @test r2a*abs(da[1].Ω - mod(-3.43864262297640e-02, 2π))     ≤ 1e-6
-        @test r2a*abs(da[2].Ω - mod(-1.58644591849564e-01, 2π))     ≤ 1e-6
+        @test r2a*abs(da[1].Ω - -3.43864262297640e-02)    ≤ 1e-6
+        @test r2a*abs(da[2].Ω - -1.58644591849564e-01)    ≤ 1e-6
     end
 
     t = [rand(49)..., 1.0]
@@ -42,7 +40,7 @@ r2a = 180 / π * 3600
             @test r2a*(abs(da.S - falp03(t[i])))                ≤ 1e-6
             @test r2a*(abs(da.F - faf03(t[i])))                 ≤ 1e-6
             @test r2a*(abs(da.D - fad03(t[i])))                 ≤ 1e-6
-            @test r2a*(abs(da.Ω - mod(faom03(t[i]), 2π)))       ≤ 1e-6
+            @test r2a*(abs(da.Ω - faom03(t[i])))                ≤ 1e-6
 
         end
     end
@@ -72,8 +70,8 @@ r2a = 180 / π * 3600
         @test r2a*abs(dab[2].D - 4.275387201216140)             ≤ 1e-6
 
         # -- Testing Mean Longitude of the Moon (< 1μas)
-        @test r2a*abs(dab[1].Ω - -3.438601115926876e-2 - 2π)    ≤ 1e-6
-        @test r2a*abs(dab[2].Ω - -1.586802211172697e-1 - 2π)    ≤ 1e-6
+        @test r2a*abs(dab[1].Ω - -3.438601115926876e-2)    ≤ 1e-6
+        @test r2a*abs(dab[2].Ω - -1.586802211172697e-1)    ≤ 1e-6
 
     end
 
@@ -87,15 +85,15 @@ end;
         pa = IERSConventions.PlanetaryArgs(rand((iers2003a, iers2010a)), t[i])
 
         # --- Planetary Arguments (accurate to 1μas)
-        @test r2a*(abs(pa.λ_Me - fame03(t[i]))) ≤ 1e-6
-        @test r2a*(abs(pa.λ_Ve - fave03(t[i]))) ≤ 1e-6
-        @test r2a*(abs(pa.λ_Ea - fae03(t[i]) )) ≤ 1e-6
-        @test r2a*(abs(pa.λ_Ma - fama03(t[i]))) ≤ 1e-6
-        @test r2a*(abs(pa.λ_Ju - faju03(t[i]))) ≤ 1e-6
-        @test r2a*(abs(pa.λ_Sa - fasa03(t[i]))) ≤ 1e-6
-        @test r2a*(abs(pa.λ_Ur - faur03(t[i]))) ≤ 1e-6
-        @test r2a*(abs(pa.λ_Ne - fane03(t[i]))) ≤ 1e-6
-        @test r2a*(abs(pa.pₐ   - fapa03(t[i]))) ≤ 1e-6
+        @test r2a*(abs(pa.λ_Me - rem2pi(fame03(t[i]), RoundNearest))) ≤ 1e-6
+        @test r2a*(abs(pa.λ_Ve - rem2pi(fave03(t[i]), RoundNearest))) ≤ 1e-6
+        @test r2a*(abs(pa.λ_Ea - rem2pi(fae03(t[i]), RoundNearest)))  ≤ 1e-6
+        @test r2a*(abs(pa.λ_Ma - rem2pi(fama03(t[i]), RoundNearest))) ≤ 1e-6
+        @test r2a*(abs(pa.λ_Ju - rem2pi(faju03(t[i]), RoundNearest))) ≤ 1e-6
+        @test r2a*(abs(pa.λ_Sa - rem2pi(fasa03(t[i]), RoundNearest))) ≤ 1e-6
+        @test r2a*(abs(pa.λ_Ur - rem2pi(faur03(t[i]), RoundNearest))) ≤ 1e-6
+        @test r2a*(abs(pa.λ_Ne - rem2pi(fane03(t[i]), RoundNearest))) ≤ 1e-6
+        @test r2a*(abs(pa.pₐ   - rem2pi(fapa03(t[i]), RoundNearest))) ≤ 1e-6
 
     end
 
diff --git a/test/sidereal.jl b/test/sidereal.jl
index 9fbc504..08efa3f 100644
--- a/test/sidereal.jl
+++ b/test/sidereal.jl
@@ -1,7 +1,8 @@
 
-r2a = 180 / π * 3600
 
 @testset "Sidereal Time" verbose=true begin 
+    
+    r2a = 180 / π * 3600
 
     for _ in 1:50 
 

From e8a49fbeee0e0d2c9c2fbbfdcfca583f9cd483bb Mon Sep 17 00:00:00 2001
From: MicheleCeresoli <michele.ceresoli@polimi.it>
Date: Sat, 10 Feb 2024 16:06:17 +0100
Subject: [PATCH 3/5] Added tests on GCRF-to-ITRF derivatives

---
 test/Project.toml |  1 +
 test/rotations.jl | 33 ++++++++++++++++++++++++++++++++-
 test/runtests.jl  |  1 +
 3 files changed, 34 insertions(+), 1 deletion(-)

diff --git a/test/Project.toml b/test/Project.toml
index c3d2347..f41d852 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,6 +1,7 @@
 [deps]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 ERFA = "17511681-8477-586a-8d98-4cfd5a1f2ec3"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LazyArtifacts = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ReferenceFrameRotations = "74f56ac7-18b3-5285-802d-d4bd4f104033"
diff --git a/test/rotations.jl b/test/rotations.jl
index a9a8573..41ad01b 100644
--- a/test/rotations.jl
+++ b/test/rotations.jl
@@ -179,13 +179,44 @@ v2as = (x, y) -> acosd(max(-1, min(1, dot(x / norm(x), y / norm(y))))) * 3600
                 Rc = iers_rot3_gcrf_to_itrf(j2000s(ep_tt), CPNc)'
                 Rd = iers_rot3_gcrf_to_itrf(j2000s(ep_tt), CPNd)' 
 
+                # Test higher order rotations 
+                R1, dR1 = iers_rot6_gcrf_to_itrf(j2000s(ep_tt), iers2010a)
+                R2, dR2, ddR2 = iers_rot9_gcrf_to_itrf(j2000s(ep_tt), iers2010a)
+                R3, dR3, ddR3, dddR3 = iers_rot12_gcrf_to_itrf(j2000s(ep_tt), iers2010a)
+
+                dR = derivative(t->iers_rot3_gcrf_to_itrf(t, iers2010a), j2000s(ep_tt))
+                ddR = derivative(
+                    t->derivative(τ->iers_rot3_gcrf_to_itrf(τ, iers2010a), t),
+                    j2000s(ep_tt)
+                )
+
+                dddR = derivative(
+                    t->derivative(
+                        τ->derivative(κ->iers_rot3_gcrf_to_itrf(κ, iers2010a), τ),
+                        t
+                    ),
+                    j2000s(ep_tt)
+                )
+
+                @test all(abs.(R1' - Ra) .≤ 1e-16)
+                @test all(abs.(R2' - Ra) .≤ 1e-16)
+                @test all(abs.(R3' - Ra) .≤ 1e-16)
+
+                @test all(abs.(dR1 - dR) .≤ 1e-11)
+                @test all(abs.(dR2 - dR) .≤ 1e-11)
+                @test all(abs.(dR3 - dR) .≤ 1e-11)
+
+                @test all(abs.(ddR2 - ddR) .≤ 1e-13)
+                @test all(abs.(ddR3 - ddR) .≤ 1e-13)
+
+                @test all(abs.(dddR3 - dddR) .≤ 1e-15)
+
                 for _ in 1:10
                     v = rand(BigFloat, 3)
                     @test v2as(R*v, Ra*v) ≤ 20e-6
                     @test v2as(Ra*v, Rb*v) ≤ 3e-3
                     @test v2as(Ra*v, Rc*v) ≤ 50e-3
                     @test v2as(Ra*v, Rd*v) ≤ 2
-
                 end
 
             end
diff --git a/test/runtests.jl b/test/runtests.jl
index 4bfcd04..c6c9ab7 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -4,6 +4,7 @@ using Test
 
 using DelimitedFiles
 using ERFA 
+using ForwardDiff: derivative
 using LazyArtifacts
 using LinearAlgebra
 using ReferenceFrameRotations

From 8558ac4448b33fc5228b099ab1f0be227f5d9acf Mon Sep 17 00:00:00 2001
From: MicheleCeresoli <michele.ceresoli@polimi.it>
Date: Sat, 10 Feb 2024 16:14:36 +0100
Subject: [PATCH 4/5] Updated derivative tests

---
 test/rotations.jl | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/test/rotations.jl b/test/rotations.jl
index 41ad01b..44d9536 100644
--- a/test/rotations.jl
+++ b/test/rotations.jl
@@ -198,10 +198,6 @@ v2as = (x, y) -> acosd(max(-1, min(1, dot(x / norm(x), y / norm(y))))) * 3600
                     j2000s(ep_tt)
                 )
 
-                @test all(abs.(R1' - Ra) .≤ 1e-16)
-                @test all(abs.(R2' - Ra) .≤ 1e-16)
-                @test all(abs.(R3' - Ra) .≤ 1e-16)
-
                 @test all(abs.(dR1 - dR) .≤ 1e-11)
                 @test all(abs.(dR2 - dR) .≤ 1e-11)
                 @test all(abs.(dR3 - dR) .≤ 1e-11)
@@ -212,11 +208,18 @@ v2as = (x, y) -> acosd(max(-1, min(1, dot(x / norm(x), y / norm(y))))) * 3600
                 @test all(abs.(dddR3 - dddR) .≤ 1e-15)
 
                 for _ in 1:10
+
                     v = rand(BigFloat, 3)
+
                     @test v2as(R*v, Ra*v) ≤ 20e-6
                     @test v2as(Ra*v, Rb*v) ≤ 3e-3
                     @test v2as(Ra*v, Rc*v) ≤ 50e-3
                     @test v2as(Ra*v, Rd*v) ≤ 2
+
+                    @test v2as(R1'*v, Ra*v) ≤ 20e-6  
+                    @test v2as(R2'*v, Ra*v) ≤ 20e-6  
+                    @test v2as(R3'*v, Ra*v) ≤ 20e-6  
+
                 end
 
             end

From 13cbada0608292d642b780543ac70e97e8d5f1ed Mon Sep 17 00:00:00 2001
From: Michele Ceresoli <85893254+MicheleCeresoli@users.noreply.github.com>
Date: Thu, 15 Feb 2024 10:20:25 +0100
Subject: [PATCH 5/5] Updated to v.1.1.2

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 17adba4..163e814 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "IERSConventions"
 uuid = "4e86e20e-879b-40dc-9e12-cee74f4cd199"
 authors = ["JSMD Team"]
-version = "1.2.0"
+version = "1.1.2"
 
 [deps]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"