JPrec.java
3.88 KB
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
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
/*
* ESO Archive
*
* $Id: JPrec.java,v 1.2 2009/02/20 23:10:11 abrighto Exp $
*
* who when what
* -------------- ---------- ----------------------------------------
* Allan Brighton 1999/05/10 port of jprec.c, Francois Ochsenbein [ESO-IPG]
*/
package jsky.coords;
/**
* Precession of Coordinates in new IAU system.
* <p/>
* This class is based on C routintes by Francois Ochsenbein [ESO-IPG].
* It uses the IAU 76 precession constants and assumes the FK5 system.
* <p/>
* Precession constants are taken from Lederle and Schwan (Astron. Astrophys.
* <b>134</b>, 1, 1984), Liske J.H. (Astron. Astrophys. <b>73</b>, 282, 1979).
* Dates must be expressed in <em>Julian Years</em>.
* <p/>
* The precession may be applied on unit vectors (mnemonic <tt>u</tt>), or
* on equatorial coordinates (mnemonic <tt>q</tt>).
*/
public class JPrec {
/**
* Compute the precession matrix, using the new IAU constants.
* (IAU 1976). The resulting matrix is such that
* u(t_1) = R * u(t_0) (old to new frame).
*
* @param R OUT: rotation matrix
* @param eq0 IN: Initial equinox (Julian Years)
* @param eq1 IN: Final equinox (Julian Years)
*/
public static void prej_R(double[][] R, double eq0, double eq1) {
double t0, dt, w;
double[] euler_angles = new double[3];
int ze = 0, theta = 1, zeta = 2;
t0 = (eq0 - 2000.e0) / 100.e0; // Origin is J2000
dt = (eq1 - eq0) / 100.e0;
w = 2306.2181e0 + (1.39656e0 - 0.000139e0 * t0) * t0; // Arc seconds
euler_angles[zeta] = (w + ((0.30188e0 - 0.000344e0 * t0) + 0.017998e0 * dt) * dt)
* dt / 3600.e0; // Degrees
euler_angles[ze] = (w + ((1.09468e0 + 0.000066e0 * t0) + 0.018203e0 * dt) * dt)
* dt / 3600.e0; // Degrees
euler_angles[theta] = ((2004.3109e0 + (-0.85330e0 - 0.000217e0 * t0) * t0)
+ ((-0.42665e0 - 0.000217e0 * t0) - 0.041833e0 * dt) * dt) * dt / 3600.e0;
// Computation of rotation matrix
Cotr.tr_Euler(euler_angles, R);
}
/**
* Performs a complete precession between 2 equinoxes.
* Use the new IAU constants.
* Compute the precession rotation matrix if necessary,
* then apply the rotation.
*
* @param q0 IN: ra+dec at equinox eq0 in degrees
* @param q1 OUT: precessed to equinox eq1
* @param eq0 IN: Initial equinox (Julian Years)
* @param eq1 IN: Final equinox (Julian Years)
*/
public static void prej_q(double[] q0, double[] q1, double eq0, double eq1) {
double[] us = new double[3];
if (eq0 == eq1) { // No precession at all, same equinox!!!
q1[0] = q0[0];
q1[1] = q0[1];
return;
}
Cotr.tr_ou(q0, us); // Convert to unit vector...
prej_u(us, us, eq0, eq1); // precess on unit vectors...
Cotr.tr_uo(us, q1); // And finally -> coordinates
}
/**
* Performs a complete precession between 2 equinoxes.
* Use the new IAU constants.
* Compute the precession rotation matrix if necessary,
* then apply the rotation.
*
* @param u0 IN: Unit vector at equinox eq0
* @param u1 OUT: precessed to equinox eq1
* @param eq0 IN: Initial equinox (Julian Years)
* @param eq1 IN: Final equinox (Julian Years)
*/
public static void prej_u(double u0[], double u1[], double eq0, double eq1) {
double[][] _r = {
{1.e0, 0.e0, 0.e0},
{0.e0, 1.e0, 0.e0},
{0.e0, 0.e0, 1.e0}
};
if (eq0 == eq1) { // No precession at all, same equinox!!!
u1[0] = u0[0];
u1[1] = u0[1];
u1[2] = u0[2];
return;
}
prej_R(_r, eq0, eq1); // Compute precession matrix
Cotr.tr_uu(u0, u1, _r); // And finally rotate...
}
}