summaryrefslogtreecommitdiffstats
path: root/org/madore/ephem/Test.java
blob: cb9dc826b8bf3d671a15dc1f434cb47d88af0405 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
package org.madore.ephem;

import java.util.Date;

public final class Test {

    public static void main(String[] args) {
	Date now;
	if ( args.length >= 1 )
	    now = new Date(Long.parseLong(args[0]));
	else
	    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);
	sun0 = Frames.Rotation.rotx(Precession.epsilon.v(Ephem.fromJd(nowTtJd))*Comput.arcsecond).apply(sun);
	System.out.format("☉ ecliptic(mean) 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);
	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);
	double frac = nowUtcSeconds/86400 - 0.5;
	double days = nowUtcMjd - 51544 + frac;
	double θ = (0.7790572732640 + 0.00273781191135448*days + frac)*2*Math.PI;
	Frames.Vector sunE = Frames.Rotation.rotz(θ + Nutation.eqorg(Ephem.fromJd(nowTtJd))*Comput.arcsecond).apply(sun);
	System.out.format("☉ astronomical lat=%.6f° long=%.6f°\n",
			  Math.atan(sunE.v[2]/Math.sqrt(sunE.v[0]*sunE.v[0]+sunE.v[1]*sunE.v[1]))/Comput.degree,
			  Math.atan2(sunE.v[1], sunE.v[0])/Comput.degree);
    }

}