diff --git a/test/j2.jl b/test/j2.jl
index 2e7705e..f49b1a5 100644
--- a/test/j2.jl
+++ b/test/j2.jl
@@ -154,6 +154,20 @@
         @test v[2] / 1000 ≈ -6.699518 atol = 1e-5
         @test v[3] / 1000 ≈ -1.266334 atol = 1e-5
         @test eltype(v) == T
+
+        r, v, j2d = j2((jd₁ - jd₀) * 86400, orb)
+
+        @test j2d isa J2Propagator{Float64, Float64}
+
+        @test r[1] / 1000 ≈ -6849.654348 atol = 5e-3
+        @test r[2] / 1000 ≈ -2253.059809 atol = 5e-3
+        @test r[3] / 1000 ≈ +3574.529667 atol = 5e-3
+        @test eltype(r) == T
+
+        @test v[1] / 1000 ≈ +1.656142 atol = 1e-5
+        @test v[2] / 1000 ≈ -6.699518 atol = 1e-5
+        @test v[3] / 1000 ≈ -1.266334 atol = 1e-5
+        @test eltype(v) == T
     end
 
     # == Float32 ===========================================================================
@@ -254,6 +268,20 @@
         @test v[2] / 1000 ≈ -6.699518 atol = 1e-3
         @test v[3] / 1000 ≈ -1.266334 atol = 1e-3
         @test eltype(v) == T
+
+        r, v, j2d = j2((jd₁ - jd₀) * 86400, orb; j2c = j2c_egm2008_f32)
+
+        @test j2d isa J2Propagator{Float64, Float32}
+
+        @test r[1] / 1000 ≈ -6849.654348 atol = 5e-1
+        @test r[2] / 1000 ≈ -2253.059809 atol = 5e-1
+        @test r[3] / 1000 ≈ +3574.529667 atol = 5e-1
+        @test eltype(r) == T
+
+        @test v[1] / 1000 ≈ +1.656142 atol = 1e-3
+        @test v[2] / 1000 ≈ -6.699518 atol = 1e-3
+        @test v[3] / 1000 ≈ -1.266334 atol = 1e-3
+        @test eltype(v) == T
     end
 end
 
diff --git a/test/j2osc.jl b/test/j2osc.jl
index d676314..52a4716 100644
--- a/test/j2osc.jl
+++ b/test/j2osc.jl
@@ -202,6 +202,19 @@
         @test results[end, 7] ≈ v[3] / 1000 atol = 1e-3
         @test eltype(r) == T
         @test eltype(v) == T
+
+        r, v, j2oscd = j2osc(results[end, 1], orb)
+
+        @test j2oscd isa J2OsculatingPropagator{Float64, Float64}
+
+        @test results[end, 2] ≈ r[1] / 1000 atol = 2e-1
+        @test results[end, 3] ≈ r[2] / 1000 atol = 2e-1
+        @test results[end, 4] ≈ r[3] / 1000 atol = 2e-1
+        @test results[end, 5] ≈ v[1] / 1000 atol = 1e-3
+        @test results[end, 6] ≈ v[2] / 1000 atol = 1e-3
+        @test results[end, 7] ≈ v[3] / 1000 atol = 1e-3
+        @test eltype(r) == T
+        @test eltype(v) == T
     end
 
     # == Float32 ===========================================================================
@@ -328,6 +341,19 @@
         @test results[end, 7] ≈ v[3] / 1000 atol = 1e-3
         @test eltype(r) == T
         @test eltype(v) == T
+
+        r, v, j2oscd = j2osc(results[end, 1], orb; j2c = j2c_egm2008_f32)
+
+        @test j2oscd isa J2OsculatingPropagator{Float64, Float32}
+
+        @test results[end, 2] ≈ r[1] / 1000 atol = 2e-1
+        @test results[end, 3] ≈ r[2] / 1000 atol = 2e-1
+        @test results[end, 4] ≈ r[3] / 1000 atol = 2e-1
+        @test results[end, 5] ≈ v[1] / 1000 atol = 1e-3
+        @test results[end, 6] ≈ v[2] / 1000 atol = 1e-3
+        @test results[end, 7] ≈ v[3] / 1000 atol = 1e-3
+        @test eltype(r) == T
+        @test eltype(v) == T
     end
 end
 
diff --git a/test/j4.jl b/test/j4.jl
index 253888d..641b78f 100644
--- a/test/j4.jl
+++ b/test/j4.jl
@@ -159,6 +159,22 @@
         @test orbk.f |> rad2deg ≈  245.617459 (atol = 4e-3)
 
         @test_broken orbk.Ω |> rad2deg ≈ 84.158846 (atol = 4e-3)
+
+        r, v, j4d = j4((jd₁ - jd₀) * 86400, orb)
+
+        @test eltype(r) == T
+        @test eltype(v) == T
+
+        orbk = j4d.orbk
+        M_k  = true_to_mean_anomaly(orbk.e, orbk.f)
+
+        @test orbk.a            ≈ 8000.0e3    (atol = 1e-3)
+        @test orbk.e            ≈    0.015    (atol = 1e-6)
+        @test orbk.i |> rad2deg ≈   28.5      (atol = 1e-6)
+        @test orbk.ω |> rad2deg ≈  225.864212 (atol = 4e-3)
+        @test orbk.f |> rad2deg ≈  245.617459 (atol = 4e-3)
+
+        @test_broken orbk.Ω |> rad2deg ≈ 84.158846 (atol = 4e-3)
     end
 
     # == Float32 ===========================================================================
@@ -245,6 +261,22 @@
         @test orbk.f |> rad2deg ≈  245.617459 (atol = 4e-3)
 
         @test_broken orbk.Ω |> rad2deg ≈ 84.158846 (atol = 4e-3)
+
+        r, v, j4d = j4((jd₁ - jd₀) * 86400, orb; j4c = j4c_egm2008_f32)
+
+        @test eltype(r) == T
+        @test eltype(v) == T
+
+        orbk = j4d.orbk
+        M_k  = true_to_mean_anomaly(orbk.e, orbk.f)
+
+        @test orbk.a            ≈ 8000.0e3    (atol = 1e-1)
+        @test orbk.e            ≈    0.015    (atol = 1e-6)
+        @test orbk.i |> rad2deg ≈   28.5      (atol = 2e-6)
+        @test orbk.ω |> rad2deg ≈  225.864212 (atol = 4e-3)
+        @test orbk.f |> rad2deg ≈  245.617459 (atol = 4e-3)
+
+        @test_broken orbk.Ω |> rad2deg ≈ 84.158846 (atol = 4e-3)
     end
 end
 
diff --git a/test/j4osc.jl b/test/j4osc.jl
index c24ccb0..333e8ac 100644
--- a/test/j4osc.jl
+++ b/test/j4osc.jl
@@ -202,6 +202,19 @@
         @test results[end, 7] ≈ v[3] / 1000 atol = 1e-3
         @test eltype(r) == T
         @test eltype(v) == T
+
+        r, v, j4oscd = j4osc(results[end, 1], orb)
+
+        @test j4oscd isa J4OsculatingPropagator{Float64, Float64}
+
+        @test results[end, 2] ≈ r[1] / 1000 atol = 2e-1
+        @test results[end, 3] ≈ r[2] / 1000 atol = 2e-1
+        @test results[end, 4] ≈ r[3] / 1000 atol = 2e-1
+        @test results[end, 5] ≈ v[1] / 1000 atol = 1e-3
+        @test results[end, 6] ≈ v[2] / 1000 atol = 1e-3
+        @test results[end, 7] ≈ v[3] / 1000 atol = 1e-3
+        @test eltype(r) == T
+        @test eltype(v) == T
     end
 
     # == Float32 ===========================================================================
@@ -328,6 +341,19 @@
         @test results[end, 7] ≈ v[3] / 1000 atol = 1e-3
         @test eltype(r) == T
         @test eltype(v) == T
+
+        r, v, j4oscd = j4osc(results[end, 1], orb; j4c = j4c_egm2008_f32)
+
+        @test j4oscd isa J4OsculatingPropagator{Float64, Float32}
+
+        @test results[end, 2] ≈ r[1] / 1000 atol = 2e-1
+        @test results[end, 3] ≈ r[2] / 1000 atol = 2e-1
+        @test results[end, 4] ≈ r[3] / 1000 atol = 2e-1
+        @test results[end, 5] ≈ v[1] / 1000 atol = 1e-3
+        @test results[end, 6] ≈ v[2] / 1000 atol = 1e-3
+        @test results[end, 7] ≈ v[3] / 1000 atol = 1e-3
+        @test eltype(r) == T
+        @test eltype(v) == T
     end
 end
 
diff --git a/test/twobody.jl b/test/twobody.jl
index 799ef06..ac85a84 100644
--- a/test/twobody.jl
+++ b/test/twobody.jl
@@ -152,6 +152,20 @@
         @test v[2] / 1000 ≈ -1.916735 atol = 1e-6
         @test v[3] / 1000 ≈ -6.112511 atol = 1e-6
         @test eltype(v) == T
+
+        r, v, tbd = twobody(40 * 60, orb)
+
+        @test tbd isa TwoBodyPropagator{Float64, Float64}
+
+        @test r[1] / 1000 ≈ -4219.7527 atol = 1e-3
+        @test r[2] / 1000 ≈ +4363.0292 atol = 1e-3
+        @test r[3] / 1000 ≈ -3958.7666 atol = 1e-3
+        @test eltype(r) == T
+
+        @test v[1] / 1000 ≈ +3.689866 atol = 1e-6
+        @test v[2] / 1000 ≈ -1.916735 atol = 1e-6
+        @test v[3] / 1000 ≈ -6.112511 atol = 1e-6
+        @test eltype(v) == T
     end
 
     # == Float32 ===========================================================================
@@ -252,6 +266,20 @@
         @test v[2] / 1000 ≈ -1.916735 atol = 5e-3
         @test v[3] / 1000 ≈ -6.112511 atol = 5e-3
         @test eltype(v) == T
+
+        r, v, tbd = twobody(40 * 60, orb; m0 = tbc_m0_f32)
+
+        @test tbd isa TwoBodyPropagator{Float64, Float32}
+
+        @test r[1] / 1000 ≈ -4219.7527 atol = 5e-1
+        @test r[2] / 1000 ≈ +4363.0292 atol = 5e-1
+        @test r[3] / 1000 ≈ -3958.7666 atol = 5e-1
+        @test eltype(r) == T
+
+        @test v[1] / 1000 ≈ +3.689866 atol = 5e-3
+        @test v[2] / 1000 ≈ -1.916735 atol = 5e-3
+        @test v[3] / 1000 ≈ -6.112511 atol = 5e-3
+        @test eltype(v) == T
     end
 end