jd2cal (#23)

* jd2cal

* tcbtdb (#25)

* Clean up and fix dtdb

* Fast and low accuracy TDB-TT difference computation. See #27

* Added test for new dtdb routine. Fix comment.

* Use new high-performance, low-precision TT<->TDB conversion

... and minor cleanups

* Add possibility to use custom offsets

* Add precision-preserving ops from AstroPy

* jd2cal

* fixed

src/constants.jl

@@ -4,7 +4,7 @@ export MJD, J2000, J1950,
     HOURS_PER_DAY, HOURS_PER_YEAR, HOURS_PER_CENTURY,
     DAYS_PER_YEAR, DAYS_PER_CENTURY,
     YEARS_PER_CENTURY,
-    OFFSET_TT_TAI, MOD_JD_77, ELG, fairhd, DAYS_PER_MILLENNIUM, TDB0, ELB
+    OFFSET_TT_TAI, MOD_JD_77, ELG, fairhd, DAYS_PER_MILLENNIUM, TDB0, ELB, JD_MAX, JD_MIN
 
 const MJD = 2400000.5
 const J2000 = Dates.datetime2julian(DateTime(2000, 1, 1, 12, 0, 0))
@@ -38,4 +38,8 @@ const DAYS_PER_MILLENNIUM = 365250.0
 
 const ELB = 1.550519768e-8
 const TDB0 = -6.55e-5
+
+const JD_MIN = -68569.5
+const JD_MAX = 1e9
+
 include("fairhd.jl")

src/conversions.jl

@@ -460,7 +460,6 @@ julia> AstroTime.Epochs.ut1tt(ut1.jd1, ut1.jd2, AstroTime.Epochs.deltat(ut1))
      date, date1
 end
 
-
 """
     tcbtdb(jd1, jd2)
 
@@ -490,6 +489,46 @@ julia> AstroTime.Epochs.tcbtdb(tcb.jd1, tcb.jd2)
     date, date1
 end
 
+@inline function jd2cal(jd1, jd2)
+    dj = jd1 + jd2
+    if dj < JD_MIN || dj > JD_MAX
+        throw(ArgumentError("Julian date is outside of the representable range ($JD_MIN, $JD_MAX)."))
+    end
+
+    if jd1 >= jd2
+        date = jd1
+        date1 = jd2
+    else
+        date = jd2
+        date1 = jd1
+    end
+
+    date1 -= 0.5
+
+    f1 = mod(date, 1.0)
+    f2 = mod(date1, 1.0)
+    f = mod(f1 + f2, 1.0)
+    if f < 0.0
+        f += 1.0
+    end
+    d = round(date-f1) + round(date1-f2) + round(f1+f2-f)
+    jd = round(d) + 1
+
+    l = jd + 68569
+    n = (4 * l) ÷ 146097
+    l -= (146097 * n + 3) ÷ 4
+    i = (4000 * (l + 1)) ÷ 1461001
+    l -= (1461 * i) ÷ 4 - 31
+    k = (80 * l) ÷ 2447.
+    id = Int(floor((l - (2447 * k) ÷ 80)))
+    l = k / 11
+    im = Int(floor((k + 2 - 12 * l)))
+    iy = Int(floor((100 * (n - 49) + i + l)))
+
+    iy, im, id, f
+
+end
+
 
 # TAI <-> UTC
 @transform UTC TAI ep begin

test/runtests.jl

@@ -182,8 +182,7 @@ AstroTime.update()
         @test_broken Epochs.utctai(julian1(leap), julian2(leap)) == ERFA.utctai(julian1(leap), julian2(leap))
         @test_broken Epochs.taiutc(tai1, tai2) == ERFA.taiutc(tai1, tai2)
 
-        @test Epochs.diff_tdb_tt(julian1(tdb), julian2(tdb), 1.0, 2.0, 3.0, 4.0) ==
-            ERFA.dtdb(julian1(tdb), julian2(tdb), 1.0, 2.0, 3.0, 4.0)
+        @test Epochs.diff_tdb_tt(julian1(tdb), julian2(tdb), 1.0, 2.0, 3.0, 4.0) == ERFA.dtdb(julian1(tdb), julian2(tdb), 1.0, 2.0, 3.0, 4.0)
         @test Epochs.tdbtt(julian1(tdb), julian2(tdb), Δtr(tdb)) == ERFA.tdbtt(julian1(tdb), julian2(tdb), Δtr(tdb))
         @test Epochs.tdbtt(julian2(tdb), julian1(tdb), Δtr(tdb)) == ERFA.tdbtt(julian2(tdb), julian1(tdb), Δtr(tdb))
         @test Epochs.tttdb(julian1(tt), julian2(tt), Δtr(tdb)) == ERFA.tttdb(julian1(tt), julian2(tt), Δtr(tdb))
@@ -195,13 +194,14 @@ AstroTime.update()
         @test Epochs.ut1tt(julian1(ut1), julian2(ut1), dt(ut1)) == ERFA.tttdb(julian1(ut1), julian2(ut1), dt(ut1))
         @test Epochs.ut1tt(julian2(ut1), julian1(ut1), dt(ut1)) == ERFA.tttdb(julian2(ut1), julian1(ut1), dt(ut1))
 
-
         @test Epochs.tcbtdb(julian1(tcb), julian2(tcb)) == ERFA.tcbtdb(julian1(tcb), julian2(tcb))
         @test Epochs.tcbtdb(julian2(tcb), julian1(tcb)) == ERFA.tcbtdb(julian2(tcb), julian1(tcb))
-
         for jd in 2414105.0:10.0:2488985.0
             @test Epochs.diff_tdb_tt(jd, 0.5) ≈ Epochs.diff_tdb_tt(jd,0.5,0.0,0.0,0.0,0.0) atol=40e-6
         end
+
+        @test Epochs.jd2cal(julian1(tt), julian2(tt)) == ERFA.jd2cal(julian1(tt), julian2(tt))
+        @test Epochs.jd2cal(julian2(tt), julian1(tt)) == ERFA.jd2cal(julian2(tt), julian1(tt))
     end
     @testset "Leap Seconds" begin
         @test leapseconds(TTEpoch(1959,1,1)) == 0