From 12ab2221611336f6d209d643574fdc2d2b1fac16 Mon Sep 17 00:00:00 2001 From: "David A. Madore" Date: Tue, 17 Apr 2012 21:09:36 +0200 Subject: Various small fixes + add a simple test (computing solar coordinates). --- org/madore/ephem/Ephem.java | 4 +++- org/madore/ephem/Frames.java | 4 ++-- org/madore/ephem/Test.java | 44 ++++++++++++++++++++++++++++++++++ org/madore/ephem/Time.java | 56 ++++++++++++++++++++++++++++++++++++++++++++ org/madore/ephem/VSOP87.java | 4 ++++ 5 files changed, 109 insertions(+), 3 deletions(-) create mode 100644 org/madore/ephem/Test.java create mode 100644 org/madore/ephem/Time.java (limited to 'org/madore') 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); + } + } -- cgit v1.2.3