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);
}
}
|
