1
/****************************************************************************
2
Xplanet 0.94 - render an image of a planet into an X window
3
Copyright (C) 2002 Hari Nair <hari@alumni.caltech.edu>
5
This program is free software; you can redistribute it and/or modify
6
it under the terms of the GNU General Public License as published by
7
the Free Software Foundation; either version 2 of the License, or
8
(at your option) any later version.
10
This program is distributed in the hope that it will be useful,
11
but WITHOUT ANY WARRANTY; without even the implied warranty of
12
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13
GNU General Public License for more details.
15
You should have received a copy of the GNU General Public License
16
along with this program; if not, write to the Free Software
17
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
18
****************************************************************************/
20
/* Based on Chapter 30 of Astronomical Formulae for Calculators by Meeus */
27
static double S(const double T)
29
return sin(T*deg_to_rad);
32
static double C(const double T)
34
return cos(T*deg_to_rad);
38
moonpos(double &moonlon, double &moonlat, double &moondist, const double T)
40
double M = poly(358.475833, 35999.04975, -1.500E-4, -3.30E-6, T);
41
double Lp = poly(270.434164, 481267.8831, -1.133E-3, 1.90E-6, T);
42
double Mp = poly(296.104608, 477198.8491, 9.192E-3, 1.44E-5, T);
43
double D = poly(350.737486, 445267.1142, -1.436E-3, 1.90E-6, T);
44
double F = poly( 11.250889, 483202.0251, -3.211E-3, -3.00E-7, T);
45
double omega = poly(259.183275,-1934.142, 2.078E-3, 2.2E-6, T);
46
Lp = Lp + ( (2.330E-4) * S(51.2 + 20.2*T)
47
+ (1.964E-3) * S(omega)
48
+ (3.964E-3) * S(poly (346.56, 132.87, -9.1731E-3, 0, T)));
49
M = M - (1.778E-3) * S(51.2 + 20.2*T);
50
Mp = Mp + ( (8.170E-4) * S(51.2 + 20.2*T)
51
+ (2.541E-3) * S(omega)
52
+ (3.964E-3) * S(poly (346.56, 132.87, -9.1731E-3, 0, T)));
53
D = D + ( (2.011E-3) * S(51.2 + 20.2*T)
54
+ (1.964E-3) * S(omega)
55
+ (3.964E-3) * S(poly (346.56, 132.87, -9.1731E-3, 0, T)));
56
F = F - ( (2.4691E-2) * S(omega)
57
- (4.3280E-3) * S(omega + 275.05 - 2.3*T)
58
+ (3.9640E-3) * S(poly (346.56, 132.87, -9.1731E-3, 0, T)));
59
double E = 1 - T * (2.495E-3 + 7.52E-6 * T);
63
+ 1.274018 * S(2*D - Mp) + .658309 * S(2*D)
64
+ .213616 * S(2*Mp) - E * .185596 * S(M)
65
- .114336*S(2*F) + .58793E-1*S(2*(D - Mp))
66
+ E*.57212E-1*S(2*D - M - Mp)
67
+ .5332E-1*S(2*D + Mp)
68
+ E*.45874E-1*S(2*D - M)
69
+ E*.41024E-1*S(Mp - M) - .34718E-1*S(D)
70
- E*.30465E-1*S(M + Mp) + .15326E-1*S(2*(D - F))
71
- .12528E-1*S(2*F - Mp) - .1098E-1*S(2*F - Mp)
72
+ .10674E-1*S(4*D - Mp) + .10034E-1*S(3*Mp)
73
+ .8548E-2*S(4*D - 2*Mp)
74
- E*.791E-2*S(M - Mp + 2*D)
75
- E*.6783E-2*S(2*D + M) + .5162E-2*S(Mp - D)
77
+ E*.4049E-2*S(Mp - M + 2*D)
78
+ .3996E-2*S(2*(Mp + D)) + .3862E-2*S(4*D)
79
+ .3665E-2*S(2*D - 3*Mp)
80
+ E*.2695E-2*S(2*Mp - M)
81
+ .2602E-2*S(Mp - 2*F - 2*D)
82
+ E*.2396E-2*S(2*D - M - 2*Mp) - .2349E-2*S(Mp+D)
83
+ E2 * .2249E-2*S(2*(D - M))
84
- E*.2125E-2*S(2*Mp + M) - E2 * .2079E-2*S(2*M)
85
+ E2 * .2059E-2*S(2*D - Mp - 2*M)
86
- .1773E-2*S(Mp + 2*D - 2*F)
87
- .1595E-2*S(2*(F + D))
88
+ E*.122E-2*S(4*D - M - Mp)
89
- .111E-2*S(2*(Mp + F)) + .892E-3*S(Mp - 3*D)
90
- E*.811E-3*S(M + Mp + 2*D)
91
+ E*.761E-3*S(4*D - M - 2*Mp)
92
+ E2 * .717E-3*S(Mp - 2*M)
93
+ E2 * .704E-3*S(Mp - 2*M - 2*D)
94
+ E*.693E-3*S(M - 2*Mp + 2*D)
95
+ E*.598E-3*S(2*D - M - 2*F) + .55E-3*S(Mp + 4*D)
96
+ .538E-3*S(4*Mp) + E*.521E-3*S(4*D - M)
97
+ .486E-3*S(2*Mp - D));
98
double B = (5.128189*S(F) + .280606*S(Mp + F)
99
+ .277693*S(Mp - F) + .173238*S(2*D - F)
100
+ .55413E-1*S(2*D + F - Mp)
101
+ .46272E-1*S(2*D - F - Mp)
102
+ .32573E-1*S(2*D + F) + .17198E-1*S(2*Mp + F)
103
+ .9267E-2*S(2*D + Mp - F) + .8823E-2*S(2*Mp - F)
104
+ E*.8247E-2*S(2*D - M - F)
105
+ .4323E-2*S(2*D - F - 2*Mp)
106
+ .42E-2*S(2*D + F + Mp)
107
+ E*.3372E-2*S(F - M - 2*D)
108
+ E*.2472E-2*S(2*D + F - M - Mp)
109
+ E*.2222E-2*S(2*D + F - M)
110
+ E*.2072E-2*S(2*D - F - M - Mp)
111
+ E*.1877E-2*S(F - M + Mp)
112
+ .1828E-2*S(4*D - F - Mp)
113
- E*.1803E-2*S(F + M) - .175E-2*S(3*F)
114
+ E*.157E-2*S(Mp - M - F) - .1487E-2*S(F + D)
115
- E*.1481E-2*S(F + M + Mp)
116
+ E*.1417E-2*S(F - M - Mp)
117
+ E*.135E-2*S(F - M) + .133E-2*S(F - D)
118
+ .1106E-2*S(F + 3*Mp) + .102E-2*S(4*D - F)
119
+ .833E-3*S(F + 4*D - Mp) + .781E-3*S(Mp - 3*F)
120
+ .67E-3*S(F + 4*D - 2*Mp) + .606E-3*S(2*D - 3*F)
121
+ .597E-3*S(2*D + 2*Mp - F)
122
+ E *.492E-3*S(2*D + Mp - M - F)
123
+ .45E-3*S(2*Mp - F - 2*D) + .439E-3*S(3*Mp - F)
124
+ .423E-3*S(F + 2*D + 2*Mp)
125
+ .422E-3*S(2*D - F - 3*Mp)
126
- E*.367E-3*S(M + F + 2*D - Mp)
127
- E*.353E-3*S(M + F + 2*D) + .331E-3*S(F + 4*D)
128
+ E*.317E-3*S(2*D + F - M + Mp)
129
+ E2 * .306E-3*S(2*D - 2*M - F)
130
- .283E-3*S(Mp + 3*F));
131
double omega1 = .4664E-3*C(omega);
132
double omega2 = .754E-4*C(omega + 275.05 - 2.3*T);
133
double beta = B*(1 - omega1 - omega2);
134
double parallax = (.9507234 + .51818E-1*C(Mp)
135
+ .9531E-2*C(2*D - Mp) + .7843E-2*C(2*D)
136
+ .2824E-2*C(2*Mp) + .857E-3*C(2*D + Mp)
137
+ E*.533E-3*C(2*D - M)
138
+ E*.401E-3*C(2*D - M - Mp)
139
+ E*.32E-3*C(Mp - M) - .271E-3*C(D)
140
- .264E-3*C(M + Mp) - .198E-3*C(2*F - Mp)
141
+ .173E-3*C(3*Mp) + .167E-3*C(4*D - Mp)
142
- E*.111E-3*C(M) + .103E-3*C(4*D - 2*Mp)
143
- .84E-4*C(2*(Mp - D)) - E*.83E-4*C(2*D + M)
144
+ .79E-4*C(2*(D + Mp)) + .72E-4*C(4*D)
145
+ E*.64E-4*C(2*D - M + Mp)
146
- E*.63E-4*C(2*D + M - Mp)
147
+ E*.41E-4*C(M + D) + E*.35E-4*C(2*Mp - M)
148
- .33E-4*C(3*Mp - 2*D) - .3E-4*C(Mp + D)
149
- .29E-4*C(2*(F - D)) - E*.29E-4*C(2*Mp + M)
150
+ E2 * .26E-4*C(2*(D - M))
151
- .23E-4* C(2*F - 2*D + Mp)
152
+ E*.19E-4*C(4*D - M - Mp));
154
moonlon = lambda * deg_to_rad;
155
moonlat = beta * deg_to_rad;
156
moondist = 1/S(parallax);