source: trunk/MagicSoft/slalib/moon.c

Last change on this file was 731, checked in by tbretz, 24 years ago
*** empty log message ***
  • Property svn:executable set to *
File size: 12.4 KB
Line 
1#include "slalib.h"
2#include "slamac.h"
3void slaMoon ( int iy, int id, float fd, float pv[6] )
4/*
5** - - - - - - - -
6** s l a M o o n
7** - - - - - - - -
8**
9** Approximate geocentric position and velocity of the Moon
10** (single precision).
11**
12** Given:
13** iy int year
14** id int day in year (1 = Jan 1st)
15** fd float fraction of day
16**
17** Returned:
18** pv float[6] Moon position & velocity vector
19**
20** Notes:
21**
22** 1 The date and time is TDB (loosely ET) in a Julian calendar
23** which has been aligned to the ordinary Gregorian
24** calendar for the interval 1900 March 1 to 2100 February 28.
25** The year and day can be obtained by calling slaCalyd or
26** slaClyd.
27**
28** 2 The Moon 6-vector is Moon centre relative to Earth centre,
29** mean equator and equinox of date. Position part, pv[0-2],
30** is in AU; velocity part, pv[3-5], is in AU/sec.
31**
32** 3 The position is accurate to better than 0.5 arcminute
33** in direction and 1000 km in distance. The velocity
34** is accurate to better than 0.5"/hour in direction and
35** 4 m/s in distance. (RMS figures with respect to JPL DE200
36** for the interval 1960-2025 are 14 arcsec and 0.2 arcsec/hour in
37** longitude, 9 arcsec and 0.2 arcsec/hour in latitude, 350 km and
38** 2 m/s in distance.) Note that the distance accuracy is
39** comparatively poor because this routine is principally intended
40** for computing topocentric direction.
41**
42** 4 This routine is only a partial implementation of the original
43** Meeus algorithm (reference below), which offers 4 times the
44** accuracy in direction and 30 times the accuracy in distance
45** when fully implemented (as it is in slaDmoon).
46**
47** Reference:
48** Meeus, l'Astronomie, June 1984, p348.
49**
50** Defined in slamac.h: dmod
51**
52** Last revision: 9 April 1998
53**
54** Copyright P.T.Wallace. All rights reserved.
55*/
56
57#define D2R 0.01745329252f /* Degrees to radians */
58
59#define RATCON 9.652743551e-12f /* Rate conversion factor: */
60 /* D2R * D2R / (86400 * 365.25) */
61
62#define ERADAU 4.2635212653763e-5f /* Earth equatorial radius in AU */
63 /* ( = 6378.137 / 149597870 ) */
64
65{
66 int iy4, n;
67 float yi, yf, t, elp, em, emp, d, f, v, dv, emn, empn, dn, fn, coeff,
68 theta, el, del, b, db, p, dp, sp, r, dr, x, y, z, xd, yd, zd,
69 sel, cel, sb, cb, rcb, rbd, w, eps, sineps, coseps;
70
71/*
72** Coefficients for fundamental arguments
73**
74** Fixed term (deg), term in t (deg & whole revs + fraction per year)
75**
76** Moon's mean longitude
77*/
78 static float elp0 = 270.434164f;
79 static float elp1 = 4812.678831f;
80 static float elp1i = 4680.0f;
81 static float elp1f = 132.678831f;
82
83/* Sun's mean anomaly */
84 static float em0 = 358.475833f;
85 static float em1 = 359.990498f;
86 static float em1f = 359.990498f;
87
88/* Moon's mean anomaly */
89 static float emp0 = 296.104608f;
90 static float emp1 = 4771.988491f;
91 static float emp1i = 4680.0f;
92 static float emp1f = 91.988491f;
93
94/* Moon's mean elongation */
95 static float d0 = 350.737486f;
96 static float d1 = 4452.671142f;
97 static float d1i = 4320.0f;
98 static float d1f = 132.671142f;
99
100/* Mean distance of the Moon from its ascending node */
101 static float f0 = 11.250889f;
102 static float f1 = 4832.020251f;
103 static float f1i = 4680.0f;
104 static float f1f = 152.020251f;
105
106/*
107** Coefficients for Moon longitude, latitude, parallax series
108*/
109 struct term {
110 float coef; /* coefficient of L, B or P term (deg) */
111 int nem; /* multiple of M in argument */
112 int nemp; /* " " M' " " */
113 int nd; /* " " D " " */
114 int nf; /* " " F " " */
115 };
116
117/*
118** Longitude coeff M M' D F
119*/
120 static struct term tl[] = { { 6.288750f, 0, 1, 0, 0 },
121 { 1.274018f, 0, -1, 2, 0 },
122 { 0.658309f, 0, 0, 2, 0 },
123 { 0.213616f, 0, 2, 0, 0 },
124 { -0.185596f, 1, 0, 0, 0 },
125 { -0.114336f, 0, 0, 0, 2 },
126 { 0.058793f, 0, -2, 2, 0 },
127 { 0.057212f, -1, -1, 2, 0 },
128 { 0.053320f, 0, 1, 2, 0 },
129 { 0.045874f, -1, 0, 2, 0 },
130 { 0.041024f, -1, 1, 0, 0 },
131 { -0.034718f, 0, 0, 1, 0 },
132 { -0.030465f, 1, 1, 0, 0 },
133 { 0.015326f, 0, 0, 2, -2 },
134 { -0.012528f, 0, 1, 0, 2 },
135 { -0.010980f, 0, -1, 0, 2 },
136 { 0.010674f, 0, -1, 4, 0 },
137 { 0.010034f, 0, 3, 0, 0 },
138 { 0.008548f, 0, -2, 4, 0 },
139 { -0.007910f, 1, -1, 2, 0 },
140 { -0.006783f, 1, 0, 2, 0 },
141 { 0.005162f, 0, 1, -1, 0 },
142 { 0.005000f, 1, 0, 1, 0 },
143 { 0.004049f, -1, 1, 2, 0 },
144 { 0.003996f, 0, 2, 2, 0 },
145 { 0.003862f, 0, 0, 4, 0 },
146 { 0.003665f, 0, -3, 2, 0 },
147 { 0.002695f, -1, 2, 0, 0 },
148 { 0.002602f, 0, 1, -2, -2 },
149 { 0.002396f, -1, -2, 2, 0 },
150 { -0.002349f, 0, 1, 1, 0 },
151 { 0.002249f, -2, 0, 2, 0 },
152 { -0.002125f, 1, 2, 0, 0 },
153 { -0.002079f, 2, 0, 0, 0 },
154 { 0.002059f, -2, -1, 2, 0 },
155 { -0.001773f, 0, 1, 2, -2 },
156 { -0.001595f, 0, 0, 2, 2 },
157 { 0.001220f, -1, -1, 4, 0 },
158 { -0.001110f, 0, 2, 0, 2 } };
159 static int NL = ( sizeof tl / sizeof ( struct term ) );
160
161/*
162** Latitude coeff M M' D F
163*/
164 static struct term tb[] = { { 5.128189f, 0, 0, 0, 1 },
165 { 0.280606f, 0, 1, 0, 1 },
166 { 0.277693f, 0, 1, 0, -1 },
167 { 0.173238f, 0, 0, 2, -1 },
168 { 0.055413f, 0, -1, 2, 1 },
169 { 0.046272f, 0, -1, 2, -1 },
170 { 0.032573f, 0, 0, 2, 1 },
171 { 0.017198f, 0, 2, 0, 1 },
172 { 0.009267f, 0, 1, 2, -1 },
173 { 0.008823f, 0, 2, 0, -1 },
174 { 0.008247f, -1, 0, 2, -1 },
175 { 0.004323f, 0, -2, 2, -1 },
176 { 0.004200f, 0, 1, 2, 1 },
177 { 0.003372f, -1, 0, -2, 1 },
178 { 0.002472f, -1, -1, 2, 1 },
179 { 0.002222f, -1, 0, 2, 1 },
180 { 0.002072f, -1, -1, 2, -1 },
181 { 0.001877f, -1, 1, 0, 1 },
182 { 0.001828f, 0, -1, 4, -1 },
183 { -0.001803f, 1, 0, 0, 1 },
184 { -0.001750f, 0, 0, 0, 3 },
185 { 0.001570f, -1, 1, 0, -1 },
186 { -0.001487f, 0, 0, 1, 1 },
187 { -0.001481f, 1, 1, 0, 1 },
188 { 0.001417f, -1, -1, 0, 1 },
189 { 0.001350f, -1, 0, 0, 1 },
190 { 0.001330f, 0, 0, -1, 1 },
191 { 0.001106f, 0, 3, 0, 1 },
192 { 0.001020f, 0, 0, 4, -1 } };
193 static int NB = ( sizeof tb / sizeof ( struct term ) );
194
195/*
196** Parallax coeff M M' D F
197*/
198 static struct term tp[] = { { 0.950724f, 0, 0, 0, 0 },
199 { 0.051818f, 0, 1, 0, 0 },
200 { 0.009531f, 0, -1, 2, 0 },
201 { 0.007843f, 0, 0, 2, 0 },
202 { 0.002824f, 0, 2, 0, 0 } };
203 static int NP = ( sizeof tp / sizeof ( struct term ) );
204
205
206
207/* Whole years & fraction of year, and years since J1900.0 */
208 yi = (float) ( iy - 1900 );
209 iy4 = iy >= 4 ? iy % 4 : 3 - ( -iy - 1 ) % 4 ;
210 yf = ( (float) ( 4 * ( id - 1 / ( iy4 + 1 ) )
211 - iy4 - 2 ) + ( 4.0f * fd ) ) / 1461.0f;
212 t = yi + yf;
213
214/* Moon's mean longitude */
215 elp = D2R * (float) dmod ( (double) ( elp0 + elp1i * yf + elp1f * t ),
216 360.0 );
217
218/* Sun's mean anomaly */
219 em = D2R * (float) dmod ( (double) ( em0 + em1f * t ), 360.0 );
220
221/* Moon's mean anomaly */
222 emp = D2R * (float) dmod ( (double) ( emp0 + emp1i * yf + emp1f * t ),
223 360.0 );
224
225/* Moon's mean elongation */
226 d = D2R * (float) dmod ( (double) ( d0 + d1i * yf + d1f * t ), 360.0 );
227
228/* Mean distance of the Moon from its ascending node */
229 f = D2R * (float) dmod ( (double) ( f0 + f1i * yf + f1f * t ), 360.0 );
230
231/* Longitude */
232 v = 0.0f;
233 dv = 0.0f;
234 for ( n = NL -1; n >= 0; n-- ) {
235 coeff = tl[n].coef;
236 emn = (float) tl[n].nem;
237 empn = (float) tl[n].nemp;
238 dn = (float) tl[n].nd;
239 fn = (float) tl[n].nf;
240 theta = emn * em + empn * emp + dn * d + fn * f;
241 v += coeff * ( (float) sin ( (double) theta ) );
242 dv += coeff * ( (float) cos ( (double) theta ) ) *
243 ( emn * em1 + empn * emp1 + dn * d1 + fn * f1 );
244 }
245 el = elp + D2R * v;
246 del = RATCON * ( elp1 / D2R + dv );
247
248/* Latitude */
249 v = 0.0f;
250 dv = 0.0f;
251 for ( n = NB - 1; n >= 0; n-- ) {
252 coeff = tb[n].coef;
253 emn = (float) tb[n].nem;
254 empn = (float) tb[n].nemp;
255 dn = (float) tb[n].nd;
256 fn = (float) tb[n].nf;
257 theta = emn * em + empn * emp + dn * d + fn * f;
258 v += coeff * ( (float) sin ( (double) theta ) );
259 dv += coeff * ( (float) cos ( (double) theta ) ) *
260 ( emn * em1 + empn * emp1 + dn * d1 + fn * f1 );
261 }
262 b = D2R * v;
263 db = RATCON * dv;
264
265/* Parallax */
266 v = 0.0f;
267 dv = 0.0f;
268 for ( n = NP - 1; n >= 0; n-- ) {
269 coeff = tp[n].coef;
270 emn = (float) tp[n].nem;
271 empn = (float) tp[n].nemp;
272 dn = (float) tp[n].nd;
273 fn = (float) tp[n].nf;
274 theta = emn * em + empn * emp + dn * d + fn * f;
275 v += coeff * ( (float) cos ( (double) theta ) );
276 dv += coeff * ( - (float) sin ( (double) theta ) ) *
277 ( emn * em1 + empn * emp1 + dn * d1 + fn * f1 );
278 }
279 p = D2R * v;
280 dp = RATCON * dv;
281
282/* Parallax to distance (AU, AU/sec) */
283 sp = (float) sin ( (double) p );
284 r = ERADAU / sp;
285 dr = - r * dp * (float) ( cos ( (double) p ) ) / sp;
286
287/* Longitude, latitude to x, y, z (AU) */
288 sel = (float) sin ( (double) el );
289 cel = (float) cos ( (double) el );
290 sb = (float) sin ( (double) b );
291 cb = (float) cos ( (double) b );
292 rcb = r * cb;
293 rbd = r * db;
294 w = rbd * sb - cb * dr;
295 x = rcb * cel;
296 y = rcb * sel;
297 z = r * sb;
298 xd = - y * del - w * cel;
299 yd = x * del - w * sel;
300 zd = rbd * cb + sb * dr;
301
302/* Mean obliquity */
303 eps = D2R * ( 23.45229f - 0.00013f * t );
304 sineps = (float) sin ( (double) eps );
305 coseps = (float) cos ( (double) eps );
306
307/* To the equatorial system, mean of date */
308 pv[0] = x;
309 pv[1] = y * coseps - z * sineps;
310 pv[2] = y * sineps + z * coseps;
311 pv[3] = xd;
312 pv[4] = yd * coseps - zd * sineps;
313 pv[5] = yd * sineps + zd * coseps;
314}
Note: See TracBrowser for help on using the repository browser.