Various small fixes + add a simple test (computing solar coordinates).

5 files changed, 109 insertions, 3 deletions

diff --git a/org/madore/ephem/Ephem.java b/org/madore/ephem/Ephem.java
index 380835c..4afbc09 100644
--- a/org/madore/ephem/Ephem.java
+++ b/org/madore/ephem/Ephem.java
@@ -2,9 +2,11 @@ package org.madore.ephem;
 
 public final class Ephem {
 
-    public static double fromJD(double jd) {
+    public static double fromJd(double jd) {
 	// Convert Julian date to Julian centuries from J2000
 	return (jd-2451545)/36525;
     }
 
+    public static final double auSeconds = 499.0047838362;  // TDB seconds!
+
 }
diff --git a/org/madore/ephem/Frames.java b/org/madore/ephem/Frames.java
index a8a4667..3671b3c 100644
--- a/org/madore/ephem/Frames.java
+++ b/org/madore/ephem/Frames.java
@@ -16,7 +16,7 @@ public final class Frames {
 		throw new IllegalArgumentException("Vector constructor expects 3 coordinates");
 	    this.v = v;
 	}
-	public Vector(Double... v) {
+	public Vector(Double[] v) {
 	    this(unboxDoubleArray(v));
 	}
 	public double dotprod(Vector w) {
@@ -40,7 +40,7 @@ public final class Frames {
 		throw new IllegalArgumentException("Rotation constructor expects 4 coordinates");
 	    this.q = q;
 	}
-	public Rotation(Double... q) {
+	public Rotation(Double[] q) {
 	    this(unboxDoubleArray(q));
 	}
 	public Vector apply(double[] v) {
diff --git a/org/madore/ephem/Test.java b/org/madore/ephem/Test.java
new file mode 100644
index 0000000..70b6663
--- /dev/null
+++ b/org/madore/ephem/Test.java
@@ -0,0 +1,44 @@
+package org.madore.ephem;
+
+import java.util.Date;
+
+public final class Test {
+
+    public static void main(String[] args) {
+	Date now = new Date();
+	long nowTime = now.getTime();
+	int nowUtcMjd = (int)(nowTime/86400000) + 40587;
+	double nowUtcSeconds = (nowTime%86400000)/1000.;
+	double nowTtJd = Time.utcToTt(nowUtcMjd, nowUtcSeconds);
+	System.out.format("UTC = %dms after Unix Epoch\nJD(TT) = %.10f\n", nowTime, nowTtJd);
+	double earthLat = VSOP87.getVal(VSOP87.Planet.EARTH, VSOP87.Variable.LAT, nowTtJd);
+	double earthLong = VSOP87.getVal(VSOP87.Planet.EARTH, VSOP87.Variable.LONG, nowTtJd);
+	double earthDist = VSOP87.getVal(VSOP87.Planet.EARTH, VSOP87.Variable.DIST, nowTtJd);
+	System.out.format("VSOP87 gives: ♁ ecliptic(J2000) lat=%.6f° long=%.6f° dist=%.6fAU\n", earthLat/Comput.degree, earthLong/Comput.degree, earthDist);
+	double delayedJd = nowTtJd - earthDist*Ephem.auSeconds/86400;
+	double sunLat = -VSOP87.getVal(VSOP87.Planet.EARTH, VSOP87.Variable.LAT, delayedJd);
+	double sunLong = Math.PI + VSOP87.getVal(VSOP87.Planet.EARTH, VSOP87.Variable.LONG, delayedJd);
+	double sunDist = VSOP87.getVal(VSOP87.Planet.EARTH, VSOP87.Variable.DIST, delayedJd);
+	System.out.format("⇒ ☉ ecliptic(J2000) lat=%.6f° long=%.6f° dist=%.6fAU\n", sunLat/Comput.degree, sunLong/Comput.degree, sunDist);
+	Frames.Vector sun0 = new Frames.Vector(sunDist*Math.cos(sunLat)*Math.cos(sunLong),
+					       sunDist*Math.cos(sunLat)*Math.sin(sunLong),
+					       sunDist*Math.sin(sunLat));
+	System.out.format("Check: ☉ ecliptic(J2000) lat=%.6f° long=%.6f°\n",
+			  Math.atan(sun0.v[2]/Math.sqrt(sun0.v[0]*sun0.v[0]+sun0.v[1]*sun0.v[1]))/Comput.degree,
+			  Math.atan2(sun0.v[1], sun0.v[0])/Comput.degree);
+	Frames.Vector sun;
+	sun = Frames.Rotation.rotx(-Precession.epsilon0*Comput.arcsecond).apply(sun0);
+	System.out.format("☉ equatorial(J2000) lat=%.6f° long=%.6f°\n",
+			  Math.atan(sun.v[2]/Math.sqrt(sun.v[0]*sun.v[0]+sun.v[1]*sun.v[1]))/Comput.degree,
+			  Math.atan2(sun.v[1], sun.v[0])/Comput.degree);
+	sun = Precession.matrix(Ephem.fromJd(nowTtJd)).apply(sun0);
+	System.out.format("☉ equatorial(mean) lat=%.6f° long=%.6f°\n",
+			  Math.atan(sun.v[2]/Math.sqrt(sun.v[0]*sun.v[0]+sun.v[1]*sun.v[1]))/Comput.degree,
+			  Math.atan2(sun.v[1], sun.v[0])/Comput.degree);
+	sun = Nutation.matrix(Ephem.fromJd(nowTtJd)).apply(sun);
+	System.out.format("☉ equatorial(true) lat=%.6f° long=%.6f°\n",
+			  Math.atan(sun.v[2]/Math.sqrt(sun.v[0]*sun.v[0]+sun.v[1]*sun.v[1]))/Comput.degree,
+			  Math.atan2(sun.v[1], sun.v[0])/Comput.degree);
+    }
+
+}
diff --git a/org/madore/ephem/Time.java b/org/madore/ephem/Time.java
new file mode 100644
index 0000000..3d2cb3e
--- /dev/null
+++ b/org/madore/ephem/Time.java
@@ -0,0 +1,56 @@
+package org.madore.ephem;
+
+public final class Time {
+
+    public static final class LeapSecond {
+	public int mjd;
+	public int offset;
+	public LeapSecond(int mjd, int offset) {
+	    this.mjd = mjd;
+	    this.offset = offset;
+	}
+    }
+
+    public static final LeapSecond[] leapSecondsTable = new LeapSecond[] {
+	new LeapSecond(41317, 10),
+	new LeapSecond(41499, 11),
+	new LeapSecond(41683, 12),
+	new LeapSecond(42048, 13),
+	new LeapSecond(42413, 14),
+	new LeapSecond(42778, 15),
+	new LeapSecond(43144, 16),
+	new LeapSecond(43509, 17),
+	new LeapSecond(43874, 18),
+	new LeapSecond(44239, 19),
+	new LeapSecond(44786, 20),
+	new LeapSecond(45151, 21),
+	new LeapSecond(45516, 22),
+	new LeapSecond(46247, 23),
+	new LeapSecond(47161, 24),
+	new LeapSecond(47892, 25),
+	new LeapSecond(48257, 26),
+	new LeapSecond(48804, 27),
+	new LeapSecond(49169, 28),
+	new LeapSecond(49534, 29),
+	new LeapSecond(50083, 30),
+	new LeapSecond(50630, 31),
+	new LeapSecond(51179, 32),
+	new LeapSecond(53736, 33),
+	new LeapSecond(54832, 34),
+	new LeapSecond(56109, 35),
+    };
+
+    public static double utcToTt(int utcMjd, double utcSeconds) {
+	// Returns JD in TT
+	int i0 = 0;  int i1 = leapSecondsTable.length;
+	while ( i1-i0 > 1 ) {
+	    int i = (i0+i1)/2;
+	    if ( utcMjd >= leapSecondsTable[i].mjd )
+		i0 = i;
+	    else
+		i1 = i;
+	}
+	return ( utcSeconds + leapSecondsTable[i0].offset + 32.184 )/86400. + utcMjd + 2400000.5;
+    }
+
+}
diff --git a/org/madore/ephem/VSOP87.java b/org/madore/ephem/VSOP87.java
index 3f06cf8..623401c 100644
--- a/org/madore/ephem/VSOP87.java
+++ b/org/madore/ephem/VSOP87.java
@@ -72,4 +72,8 @@ public final class VSOP87 {
 	return data.get(pl).get(v);
     }
 
+    public static double getVal(Planet pl, Variable v, double tdbJd) {
+	return getFunc(pl, v).v(Ephem.fromJd(tdbJd)/10);
+    }
+
 }