1 /* SWISSEPH
2 *
3 * Steve Moshier's analytical lunar ephemeris
4
5 **************************************************************/
6 /* Copyright (C) 1997 - 2021 Astrodienst AG, Switzerland. All rights reserved.
7
8 License conditions
9 ------------------
10
11 This file is part of Swiss Ephemeris.
12
13 Swiss Ephemeris is distributed with NO WARRANTY OF ANY KIND. No author
14 or distributor accepts any responsibility for the consequences of using it,
15 or for whether it serves any particular purpose or works at all, unless he
16 or she says so in writing.
17
18 Swiss Ephemeris is made available by its authors under a dual licensing
19 system. The software developer, who uses any part of Swiss Ephemeris
20 in his or her software, must choose between one of the two license models,
21 which are
22 a) GNU Affero General Public License (AGPL)
23 b) Swiss Ephemeris Professional License
24
25 The choice must be made before the software developer distributes software
26 containing parts of Swiss Ephemeris to others, and before any public
27 service using the developed software is activated.
28
29 If the developer choses the AGPL software license, he or she must fulfill
30 the conditions of that license, which includes the obligation to place his
31 or her whole software project under the AGPL or a compatible license.
32 See https://www.gnu.org/licenses/agpl-3.0.html
33
34 If the developer choses the Swiss Ephemeris Professional license,
35 he must follow the instructions as found in http://www.astro.com/swisseph/
36 and purchase the Swiss Ephemeris Professional Edition from Astrodienst
37 and sign the corresponding license contract.
38
39 The License grants you the right to use, copy, modify and redistribute
40 Swiss Ephemeris, but only under certain conditions described in the License.
41 Among other things, the License requires that the copyright notices and
42 this notice be preserved on all copies.
43
44 Authors of the Swiss Ephemeris: Dieter Koch and Alois Treindl
45
46 The authors of Swiss Ephemeris have no control or influence over any of
47 the derived works, i.e. over software or services created by other
48 programmers which use Swiss Ephemeris functions.
49
50 The names of the authors or of the copyright holder (Astrodienst) must not
51 be used for promoting any software, product or service which uses or contains
52 the Swiss Ephemeris. This copyright notice is the ONLY place where the
53 names of the authors can legally appear, except in cases where they have
54 given special permission in writing.
55
56 The trademarks 'Swiss Ephemeris' and 'Swiss Ephemeris inside' may be used
57 for promoting such software, products or services.
58 */
59
60
61 /*
62 * Expansions for the geocentric ecliptic longitude,
63 * latitude, and distance of the Moon referred to the mean equinox
64 * and ecliptic of date.
65 *
66 * This version of cmoon.c adjusts the ELP2000-85 analytical Lunar
67 * theory of Chapront-Touze and Chapront to fit the Jet Propulsion
68 * Laboratory's DE404 long ephemeris on the interval from 3000 B.C.
69 * to 3000 A.D.
70 *
71 * The fit is much better in the remote past and future if
72 * secular terms are included in the arguments of the oscillatory
73 * perturbations. Such adjustments cannot easily be incorporated
74 * into the 1991 lunar tables. In this program the traditional
75 * literal arguments are used instead, with mean elements adjusted
76 * for a best fit to the reference ephemeris.
77 *
78 * This program omits many oscillatory terms from the analytical
79 * theory which, if they were included, would yield a much higher
80 * accuracy for modern dates. Detailed statistics of the precision
81 * are given in the table below. Comparing at 64-day intervals
82 * over the period -3000 to +3000, the maximum discrepancies noted
83 * were 7" longitude, 5" latitude, and 5 x 10^-8 au radius.
84 * The expressions used for precession in this comparision were
85 * those of Simon et al (1994).
86 *
87 * The adjusted coefficients were found by an unweighted least squares
88 * fit to the numerical ephemeris in the mentioned test interval.
89 * The approximation error increases rapidly outside this interval.
90 * J. Chapront (1994) has described the basic fitting procedure.
91 *
92 * A major change from DE200 to DE404 is in the coefficient
93 * of tidal acceleration of the Moon, which causes the Moon's
94 * longitude to depart by about -0.9" per century squared
95 * from DE200. Uncertainty in this quantity continues to
96 * be the limiting factor in long term projections of the Moon's
97 * ephemeris.
98 *
99 * Since the Lunar theory is cast in the ecliptic of date, it makes
100 * some difference what formula you use for precession. The adjustment
101 * to DE404 was carried out relative to the mean equinox and ecliptic
102 * of date as defined in Williams (1994). An earlier version of this
103 * program used the precession given by Simon et al (1994). The difference
104 * between these two precession formulas amounts to about 12" in Lunar
105 * longitude at 3000 B.C.
106 *
107 * Maximum deviations between DE404 and this program
108 * in a set of 34274 samples spaced 64 days apart
109 *
110 * Interval Longitude Latitude Radius
111 * Julian Year arc sec arc sec 10^-8 au
112 * -3000 to -2500 5.66 4.66 4.93
113 * -2500 to -2000 5.49 3.98 4.56
114 * -2000 to -1500 6.98 4.17 4.81
115 * -1500 to -1000 5.74 3.53 4.87
116 * -1000 to -500 5.95 3.42 4.67
117 * -500 to 0 4.94 3.07 4.04
118 * 0 to 500 4.42 2.65 4.55
119 * 500 to 1000 5.68 3.30 3.99
120 * 1000 to 1500 4.32 3.21 3.83
121 * 1500 to 2000 2.70 2.69 3.71
122 * 2000 to 2500 3.35 2.32 3.85
123 * 2500 to 3000 4.62 2.39 4.11
124 *
125 *
126 *
127 * References:
128 *
129 * James G. Williams, "Contributions to the Earth's obliquity rate,
130 * precession, and nutation," Astron. J. 108, 711-724 (1994)
131 *
132 * DE403 and DE404 ephemerides by E. M. Standish, X. X. Newhall, and
133 * J. G. Williams are at the JPL computer site navigator.jpl.nasa.gov.
134 *
135 * J. L. Simon, P. Bretagnon, J. Chapront, M. Chapront-Touze', G. Francou,
136 * and J. Laskar, "Numerical Expressions for precession formulae and
137 * mean elements for the Moon and the planets," Astronomy and Astrophysics
138 * 282, 663-683 (1994)
139 *
140 * P. Bretagnon and Francou, G., "Planetary theories in rectangular
141 * and spherical variables. VSOP87 solutions," Astronomy and
142 * Astrophysics 202, 309-315 (1988)
143 *
144 * M. Chapront-Touze' and J. Chapront, "ELP2000-85: a semi-analytical
145 * lunar ephemeris adequate for historical times," Astronomy and
146 * Astrophysics 190, 342-352 (1988).
147 *
148 * M. Chapront-Touze' and J. Chapront, _Lunar Tables and
149 * Programs from 4000 B.C. to A.D. 8000_, Willmann-Bell (1991)
150 *
151 * J. Laskar, "Secular terms of classical planetary theories
152 * using the results of general theory," Astronomy and Astrophysics
153 * 157, 59070 (1986)
154 *
155 * S. L. Moshier, "Comparison of a 7000-year lunar ephemeris
156 * with analytical theory," Astronomy and Astrophysics 262,
157 * 613-616 (1992)
158 *
159 * J. Chapront, "Representation of planetary ephemerides by frequency
160 * analysis. Application to the five outer planets," Astronomy and
161 * Astrophysics Suppl. Ser. 109, 181-192 (1994)
162 *
163 *
164 * Entry swi_moshmoon2() returns the geometric position of the Moon
165 * relative to the Earth. Its calling procedure is as follows:
166 *
167 * double JD; input Julian Ephemeris Date
168 * double pol[3]; output ecliptic polar coordinatees in radians and au
169 * pol[0] longitude, pol[1] latitude, pol[2] radius
170 * swi_moshmoon2( JD, pol );
171 *
172 * - S. L. Moshier, August, 1991
173 * DE200 fit: July, 1992
174 * DE404 fit: October, 1995
175 *
176 * Dieter Koch: adaptation to SWISSEPH, April 1996
177 * 18-feb-2006 replaced LP by SWELP because of name collision
178 */
179
180 #include <string.h>
181 #include "swephexp.h"
182 #include "sweph.h"
183 #include "swephlib.h"
184
185 static void mean_elements(void);
186 static void mean_elements_pl(void);
187 static double mods3600(double x);
188 static void ecldat_equ2000(double tjd, double *xpm);
189 static void chewm(const short *pt, int nlines, int nangles,
190 int typflg, double *ans );
191 static void sscc(int k, double arg, int n );
192 static void moon1(void);
193 static void moon2(void);
194 static void moon3(void);
195 static void moon4(void);
196
197
198 #ifdef MOSH_MOON_200
199 /* The following coefficients were calculated by a simultaneous least
200 * squares fit between the analytical theory and the continued DE200
201 * numerically integrated ephemeris from 9000 BC to 13000 AD.
202 * See references to the array z[] later on in the program.
203 * The 71 coefficients were estimated from 42,529 Lunar positions.
204 */
205 static const double z[] = {
206 -1.225346551567e+001, /* F, t^2 */
207 -1.096676093208e-003, /* F, t^3 */
208 -2.165750777942e-006, /* F, t^4 */
209 -2.790392351314e-009, /* F, t^5 */
210 4.189032191814e-011, /* F, t^6 */
211 4.474984866301e-013, /* F, t^7 */
212 3.239398410335e+001, /* l, t^2 */
213 5.185305877294e-002, /* l, t^3 */
214 -2.536291235258e-004, /* l, t^4 */
215 -2.506365935364e-008, /* l, t^5 */
216 3.452144225877e-011, /* l, t^6 */
217 -1.755312760154e-012, /* l, t^7 */
218 -5.870522364514e+000, /* D, t^2 */
219 6.493037519768e-003, /* D, t^3 */
220 -3.702060118571e-005, /* D, t^4 */
221 2.560078201452e-009, /* D, t^5 */
222 2.555243317839e-011, /* D, t^6 */
223 -3.207663637426e-013, /* D, t^7 */
224 -4.776684245026e+000, /* L, t^2 */
225 6.580112707824e-003, /* L, t^3 */
226 -6.073960534117e-005, /* L, t^4 */
227 -1.024222633731e-008, /* L, t^5 */
228 2.235210987108e-010, /* L, t^6 */
229 7.200592540556e-014, /* L, t^7 */
230 -8.552017636339e+001, /* t^2 cos(18V - 16E - l) */
231 -2.055794304596e+002, /* t^2 sin(18V - 16E - l) */
232 -1.097555241866e+000, /* t^3 cos(18V - 16E - l) */
233 5.219423171002e-001, /* t^3 sin(18V - 16E - l) */
234 2.088802640755e-003, /* t^4 cos(18V - 16E - l) */
235 4.616541527921e-003, /* t^4 sin(18V - 16E - l) */
236 4.794930645807e+000, /* t^2 cos(10V - 3E - l) */
237 -4.595134364283e+001, /* t^2 sin(10V - 3E - l) */
238 -6.659812174691e-002, /* t^3 cos(10V - 3E - l) */
239 -2.570048828246e-001, /* t^3 sin(10V - 3E - l) */
240 6.229863046223e-004, /* t^4 cos(10V - 3E - l) */
241 5.504368344700e-003, /* t^4 sin(10V - 3E - l) */
242 -3.084830597278e+000, /* t^2 cos(8V - 13E) */
243 -1.000471012253e+001, /* t^2 sin(8V - 13E) */
244 6.590112074510e-002, /* t^3 cos(8V - 13E) */
245 -3.212573348278e-003, /* t^3 sin(8V - 13E) */
246 5.409038312567e-004, /* t^4 cos(8V - 13E) */
247 1.293377988163e-003, /* t^4 sin(8V - 13E) */
248 2.311794636111e+001, /* t^2 cos(4E - 8M + 3J) */
249 -3.157036220040e+000, /* t^2 sin(4E - 8M + 3J) */
250 -3.019293162417e+000, /* t^2 cos(18V - 16E) */
251 -9.211526858975e+000, /* t^2 sin(18V - 16E) */
252 -4.993704215784e-002, /* t^3 cos(18V - 16E) */
253 2.991187525454e-002, /* t^3 sin(18V - 16E) */
254 -3.827414182969e+000, /* t^2 cos(18V - 16E - 2l) */
255 -9.891527703219e+000, /* t^2 sin(18V - 16E - 2l) */
256 -5.322093802878e-002, /* t^3 cos(18V - 16E - 2l) */
257 3.164702647371e-002, /* t^3 sin(18V - 16E - 2l) */
258 7.713905234217e+000, /* t^2 cos(2J - 5S) */
259 -6.077986950734e+000, /* t^3 sin(2J - 5S) */
260 -1.278232501462e-001, /* t^2 cos(L - F) */
261 4.760967236383e-001, /* t^2 sin(L - F) */
262 -6.759005756460e-001, /* t^3 sin(l') */
263 1.655727996357e-003, /* t^4 sin(l') */
264 1.646526117252e-001, /* t^3 sin(2D - l') */
265 -4.167078100233e-004, /* t^4 sin(2D - l') */
266 2.067529538504e-001, /* t^3 sin(2D - l' - l) */
267 -5.219127398748e-004, /* t^4 sin(2D - l' - l) */
268 -1.526335222289e-001, /* t^3 sin(l' - l) */
269 -1.120545131358e-001, /* t^3 sin(l' + l) */
270 4.619472391553e-002, /* t^3 sin(2D - 2l') */
271 4.863621236157e-004, /* t^4 sin(2D - 2l') */
272 -4.280059182608e-002, /* t^3 sin(2l') */
273 -4.328378207833e-004, /* t^4 sin(2l') */
274 -8.371028286974e-003, /* t^3 sin(2D - l) */
275 4.089447328174e-002, /* t^3 sin(2D - 2l' - l) */
276 -1.238363006354e-002, /* t^3 sin(2D + 2l' - l) */
277 };
278 #else
279 /* The following coefficients were calculated by a simultaneous least
280 * squares fit between the analytical theory and DE404 on the finite
281 * interval from -3000 to +3000.
282 * The coefficients were estimated from 34,247 Lunar positions.
283 */
284 static const double z[] = {
285 /* The following are scaled in arc seconds, time in Julian centuries.
286 They replace the corresponding terms in the mean elements. */
287 -1.312045233711e+01, /* F, t^2 */
288 -1.138215912580e-03, /* F, t^3 */
289 -9.646018347184e-06, /* F, t^4 */
290 3.146734198839e+01, /* l, t^2 */
291 4.768357585780e-02, /* l, t^3 */
292 -3.421689790404e-04, /* l, t^4 */
293 -6.847070905410e+00, /* D, t^2 */
294 -5.834100476561e-03, /* D, t^3 */
295 -2.905334122698e-04, /* D, t^4 */
296 -5.663161722088e+00, /* L, t^2 */
297 5.722859298199e-03, /* L, t^3 */
298 -8.466472828815e-05, /* L, t^4 */
299 /* The following longitude terms are in arc seconds times 10^5. */
300 -8.429817796435e+01, /* t^2 cos(18V - 16E - l) */
301 -2.072552484689e+02, /* t^2 sin(18V - 16E - l) */
302 7.876842214863e+00, /* t^2 cos(10V - 3E - l) */
303 1.836463749022e+00, /* t^2 sin(10V - 3E - l) */
304 -1.557471855361e+01, /* t^2 cos(8V - 13E) */
305 -2.006969124724e+01, /* t^2 sin(8V - 13E) */
306 2.152670284757e+01, /* t^2 cos(4E - 8M + 3J) */
307 -6.179946916139e+00, /* t^2 sin(4E - 8M + 3J) */
308 -9.070028191196e-01, /* t^2 cos(18V - 16E) */
309 -1.270848233038e+01, /* t^2 sin(18V - 16E) */
310 -2.145589319058e+00, /* t^2 cos(2J - 5S) */
311 1.381936399935e+01, /* t^2 sin(2J - 5S) */
312 -1.999840061168e+00, /* t^3 sin(l') */
313 };
314 #endif /* ! MOSH_MOON_200 */
315
316 /* Perturbation tables
317 */
318 #define NLR 118
319 static const short LR[8*NLR] = {
320 /*
321 Longitude Radius
322 D l' l F 1" .0001" 1km .0001km */
323
324 0, 0, 1, 0, 22639, 5858,-20905,-3550,
325 2, 0,-1, 0, 4586, 4383, -3699,-1109,
326 2, 0, 0, 0, 2369, 9139, -2955,-9676,
327 0, 0, 2, 0, 769, 257, -569,-9251,
328 0, 1, 0, 0, -666,-4171, 48, 8883,
329 0, 0, 0, 2, -411,-5957, -3,-1483,
330 2, 0,-2, 0, 211, 6556, 246, 1585,
331 2,-1,-1, 0, 205, 4358, -152,-1377,
332 2, 0, 1, 0, 191, 9562, -170,-7331,
333 2,-1, 0, 0, 164, 7285, -204,-5860,
334 0, 1,-1, 0, -147,-3213, -129,-6201,
335 1, 0, 0, 0, -124,-9881, 108, 7427,
336 0, 1, 1, 0, -109,-3803, 104, 7552,
337 2, 0, 0,-2, 55, 1771, 10, 3211,
338 0, 0, 1, 2, -45, -996, 0, 0,
339 0, 0, 1,-2, 39, 5333, 79, 6606,
340 4, 0,-1, 0, 38, 4298, -34,-7825,
341 0, 0, 3, 0, 36, 1238, -23,-2104,
342 4, 0,-2, 0, 30, 7726, -21,-6363,
343 2, 1,-1, 0, -28,-3971, 24, 2085,
344 2, 1, 0, 0, -24,-3582, 30, 8238,
345 1, 0,-1, 0, -18,-5847, -8,-3791,
346 1, 1, 0, 0, 17, 9545, -16,-6747,
347 2,-1, 1, 0, 14, 5303, -12,-8314,
348 2, 0, 2, 0, 14, 3797, -10,-4448,
349 4, 0, 0, 0, 13, 8991, -11,-6500,
350 2, 0,-3, 0, 13, 1941, 14, 4027,
351 0, 1,-2, 0, -9,-6791, -7, -27,
352 2, 0,-1, 2, -9,-3659, 0, 7740,
353 2,-1,-2, 0, 8, 6055, 10, 562,
354 1, 0, 1, 0, -8,-4531, 6, 3220,
355 2,-2, 0, 0, 8, 502, -9,-8845,
356 0, 1, 2, 0, -7,-6302, 5, 7509,
357 0, 2, 0, 0, -7,-4475, 1, 657,
358 2,-2,-1, 0, 7, 3712, -4,-9501,
359 2, 0, 1,-2, -6,-3832, 4, 1311,
360 2, 0, 0, 2, -5,-7416, 0, 0,
361 4,-1,-1, 0, 4, 3740, -3,-9580,
362 0, 0, 2, 2, -3,-9976, 0, 0,
363 3, 0,-1, 0, -3,-2097, 3, 2582,
364 2, 1, 1, 0, -2,-9145, 2, 6164,
365 4,-1,-2, 0, 2, 7319, -1,-8970,
366 0, 2,-1, 0, -2,-5679, -2,-1171,
367 2, 2,-1, 0, -2,-5212, 2, 3536,
368 2, 1,-2, 0, 2, 4889, 0, 1437,
369 2,-1, 0,-2, 2, 1461, 0, 6571,
370 4, 0, 1, 0, 1, 9777, -1,-4226,
371 0, 0, 4, 0, 1, 9337, -1,-1169,
372 4,-1, 0, 0, 1, 8708, -1,-5714,
373 1, 0,-2, 0, -1,-7530, -1,-7385,
374 2, 1, 0,-2, -1,-4372, 0,-1357,
375 0, 0, 2,-2, -1,-3726, -4,-4212,
376 1, 1, 1, 0, 1, 2618, 0,-9333,
377 3, 0,-2, 0, -1,-2241, 0, 8624,
378 4, 0,-3, 0, 1, 1868, 0,-5142,
379 2,-1, 2, 0, 1, 1770, 0,-8488,
380 0, 2, 1, 0, -1,-1617, 1, 1655,
381 1, 1,-1, 0, 1, 777, 0, 8512,
382 2, 0, 3, 0, 1, 595, 0,-6697,
383 2, 0, 1, 2, 0,-9902, 0, 0,
384 2, 0,-4, 0, 0, 9483, 0, 7785,
385 2,-2, 1, 0, 0, 7517, 0,-6575,
386 0, 1,-3, 0, 0,-6694, 0,-4224,
387 4, 1,-1, 0, 0,-6352, 0, 5788,
388 1, 0, 2, 0, 0,-5840, 0, 3785,
389 1, 0, 0,-2, 0,-5833, 0,-7956,
390 6, 0,-2, 0, 0, 5716, 0,-4225,
391 2, 0,-2,-2, 0,-5606, 0, 4726,
392 1,-1, 0, 0, 0,-5569, 0, 4976,
393 0, 1, 3, 0, 0,-5459, 0, 3551,
394 2, 0,-2, 2, 0,-5357, 0, 7740,
395 2, 0,-1,-2, 0, 1790, 8, 7516,
396 3, 0, 0, 0, 0, 4042, -1,-4189,
397 2,-1,-3, 0, 0, 4784, 0, 4950,
398 2,-1, 3, 0, 0, 932, 0, -585,
399 2, 0, 2,-2, 0,-4538, 0, 2840,
400 2,-1,-1, 2, 0,-4262, 0, 373,
401 0, 0, 0, 4, 0, 4203, 0, 0,
402 0, 1, 0, 2, 0, 4134, 0,-1580,
403 6, 0,-1, 0, 0, 3945, 0,-2866,
404 2,-1, 0, 2, 0,-3821, 0, 0,
405 2,-1, 1,-2, 0,-3745, 0, 2094,
406 4, 1,-2, 0, 0,-3576, 0, 2370,
407 1, 1,-2, 0, 0, 3497, 0, 3323,
408 2,-3, 0, 0, 0, 3398, 0,-4107,
409 0, 0, 3, 2, 0,-3286, 0, 0,
410 4,-2,-1, 0, 0,-3087, 0,-2790,
411 0, 1,-1,-2, 0, 3015, 0, 0,
412 4, 0,-1,-2, 0, 3009, 0,-3218,
413 2,-2,-2, 0, 0, 2942, 0, 3430,
414 6, 0,-3, 0, 0, 2925, 0,-1832,
415 2, 1, 2, 0, 0,-2902, 0, 2125,
416 4, 1, 0, 0, 0,-2891, 0, 2445,
417 4,-1, 1, 0, 0, 2825, 0,-2029,
418 3, 1,-1, 0, 0, 2737, 0,-2126,
419 0, 1, 1, 2, 0, 2634, 0, 0,
420 1, 0, 0, 2, 0, 2543, 0, 0,
421 3, 0, 0,-2, 0,-2530, 0, 2010,
422 2, 2,-2, 0, 0,-2499, 0,-1089,
423 2,-3,-1, 0, 0, 2469, 0,-1481,
424 3,-1,-1, 0, 0,-2314, 0, 2556,
425 4, 0, 2, 0, 0, 2185, 0,-1392,
426 4, 0,-1, 2, 0,-2013, 0, 0,
427 0, 2,-2, 0, 0,-1931, 0, 0,
428 2, 2, 0, 0, 0,-1858, 0, 0,
429 2, 1,-3, 0, 0, 1762, 0, 0,
430 4, 0,-2, 2, 0,-1698, 0, 0,
431 4,-2,-2, 0, 0, 1578, 0,-1083,
432 4,-2, 0, 0, 0, 1522, 0,-1281,
433 3, 1, 0, 0, 0, 1499, 0,-1077,
434 1,-1,-1, 0, 0,-1364, 0, 1141,
435 1,-3, 0, 0, 0,-1281, 0, 0,
436 6, 0, 0, 0, 0, 1261, 0, -859,
437 2, 0, 2, 2, 0,-1239, 0, 0,
438 1,-1, 1, 0, 0,-1207, 0, 1100,
439 0, 0, 5, 0, 0, 1110, 0, -589,
440 0, 3, 0, 0, 0,-1013, 0, 213,
441 4,-1,-3, 0, 0, 998, 0, 0,
442 };
443
444
445 #ifdef MOSH_MOON_200
446 #define NMB 56
447 static const short MB[6*NMB] = {
448 /*
449 Latitude
450 D l' l F 1" .0001" */
451
452 0, 0, 0, 1,18461, 2387,
453 0, 0, 1, 1, 1010, 1671,
454 0, 0, 1,-1, 999, 6936,
455 2, 0, 0,-1, 623, 6524,
456 2, 0,-1, 1, 199, 4837,
457 2, 0,-1,-1, 166, 5741,
458 2, 0, 0, 1, 117, 2607,
459 0, 0, 2, 1, 61, 9120,
460 2, 0, 1,-1, 33, 3572,
461 0, 0, 2,-1, 31, 7597,
462 2,-1, 0,-1, 29, 5766,
463 2, 0,-2,-1, 15, 5663,
464 2, 0, 1, 1, 15, 1216,
465 2, 1, 0,-1, -12, -941,
466 2,-1,-1, 1, 8, 8681,
467 2,-1, 0, 1, 7, 9586,
468 2,-1,-1,-1, 7, 4346,
469 0, 1,-1,-1, -6,-7314,
470 4, 0,-1,-1, 6, 5796,
471 0, 1, 0, 1, -6,-4601,
472 0, 0, 0, 3, -6,-2965,
473 0, 1,-1, 1, -5,-6324,
474 1, 0, 0, 1, -5,-3684,
475 0, 1, 1, 1, -5,-3113,
476 0, 1, 1,-1, -5, -759,
477 0, 1, 0,-1, -4,-8396,
478 1, 0, 0,-1, -4,-8057,
479 0, 0, 3, 1, 3, 9841,
480 4, 0, 0,-1, 3, 6745,
481 4, 0,-1, 1, 2, 9985,
482 0, 0, 1,-3, 2, 7986,
483 4, 0,-2, 1, 2, 4139,
484 2, 0, 0,-3, 2, 1863,
485 2, 0, 2,-1, 2, 1462,
486 2,-1, 1,-1, 1, 7660,
487 2, 0,-2, 1, -1,-6244,
488 0, 0, 3,-1, 1, 5813,
489 2, 0, 2, 1, 1, 5198,
490 2, 0,-3,-1, 1, 5156,
491 2, 1,-1, 1, -1,-3178,
492 2, 1, 0, 1, -1,-2643,
493 4, 0, 0, 1, 1, 1919,
494 2,-1, 1, 1, 1, 1346,
495 2,-2, 0,-1, 1, 859,
496 0, 0, 1, 3, -1, -194,
497 2, 1, 1,-1, 0,-8227,
498 1, 1, 0,-1, 0, 8042,
499 1, 1, 0, 1, 0, 8026,
500 0, 1,-2,-1, 0,-7932,
501 2, 1,-1,-1, 0,-7910,
502 1, 0, 1, 1, 0,-6674,
503 2,-1,-2,-1, 0, 6502,
504 0, 1, 2, 1, 0,-6388,
505 4, 0,-2,-1, 0, 6337,
506 4,-1,-1,-1, 0, 5958,
507 1, 0, 1,-1, 0,-5889,
508 };
509 #else
510 #define NMB 77
511 static const short MB[6*NMB] = {
512 /*
513 Latitude
514 D l' l F 1" .0001" */
515
516 0, 0, 0, 1,18461, 2387,
517 0, 0, 1, 1, 1010, 1671,
518 0, 0, 1,-1, 999, 6936,
519 2, 0, 0,-1, 623, 6524,
520 2, 0,-1, 1, 199, 4837,
521 2, 0,-1,-1, 166, 5741,
522 2, 0, 0, 1, 117, 2607,
523 0, 0, 2, 1, 61, 9120,
524 2, 0, 1,-1, 33, 3572,
525 0, 0, 2,-1, 31, 7597,
526 2,-1, 0,-1, 29, 5766,
527 2, 0,-2,-1, 15, 5663,
528 2, 0, 1, 1, 15, 1216,
529 2, 1, 0,-1, -12, -941,
530 2,-1,-1, 1, 8, 8681,
531 2,-1, 0, 1, 7, 9586,
532 2,-1,-1,-1, 7, 4346,
533 0, 1,-1,-1, -6,-7314,
534 4, 0,-1,-1, 6, 5796,
535 0, 1, 0, 1, -6,-4601,
536 0, 0, 0, 3, -6,-2965,
537 0, 1,-1, 1, -5,-6324,
538 1, 0, 0, 1, -5,-3684,
539 0, 1, 1, 1, -5,-3113,
540 0, 1, 1,-1, -5, -759,
541 0, 1, 0,-1, -4,-8396,
542 1, 0, 0,-1, -4,-8057,
543 0, 0, 3, 1, 3, 9841,
544 4, 0, 0,-1, 3, 6745,
545 4, 0,-1, 1, 2, 9985,
546 0, 0, 1,-3, 2, 7986,
547 4, 0,-2, 1, 2, 4139,
548 2, 0, 0,-3, 2, 1863,
549 2, 0, 2,-1, 2, 1462,
550 2,-1, 1,-1, 1, 7660,
551 2, 0,-2, 1, -1,-6244,
552 0, 0, 3,-1, 1, 5813,
553 2, 0, 2, 1, 1, 5198,
554 2, 0,-3,-1, 1, 5156,
555 2, 1,-1, 1, -1,-3178,
556 2, 1, 0, 1, -1,-2643,
557 4, 0, 0, 1, 1, 1919,
558 2,-1, 1, 1, 1, 1346,
559 2,-2, 0,-1, 1, 859,
560 0, 0, 1, 3, -1, -194,
561 2, 1, 1,-1, 0,-8227,
562 1, 1, 0,-1, 0, 8042,
563 1, 1, 0, 1, 0, 8026,
564 0, 1,-2,-1, 0,-7932,
565 2, 1,-1,-1, 0,-7910,
566 1, 0, 1, 1, 0,-6674,
567 2,-1,-2,-1, 0, 6502,
568 0, 1, 2, 1, 0,-6388,
569 4, 0,-2,-1, 0, 6337,
570 4,-1,-1,-1, 0, 5958,
571 1, 0, 1,-1, 0,-5889,
572 4, 0, 1,-1, 0, 4734,
573 1, 0,-1,-1, 0,-4299,
574 4,-1, 0,-1, 0, 4149,
575 2,-2, 0, 1, 0, 3835,
576 3, 0, 0,-1, 0,-3518,
577 4,-1,-1, 1, 0, 3388,
578 2, 0,-1,-3, 0, 3291,
579 2,-2,-1, 1, 0, 3147,
580 0, 1, 2,-1, 0,-3129,
581 3, 0,-1,-1, 0,-3052,
582 0, 1,-2, 1, 0,-3013,
583 2, 0, 1,-3, 0,-2912,
584 2,-2,-1,-1, 0, 2686,
585 0, 0, 4, 1, 0, 2633,
586 2, 0,-3, 1, 0, 2541,
587 2, 0,-1, 3, 0,-2448,
588 2, 1, 1, 1, 0,-2370,
589 4,-1,-2, 1, 0, 2138,
590 4, 0, 1, 1, 0, 2126,
591 3, 0,-1, 1, 0,-2059,
592 4, 1,-1,-1, 0,-1719,
593 };
594 #endif /* ! MOSH_MOON_200 */
595
596 #define NLRT 38
597 static const short LRT[8*NLRT] = {
598 /*
599 Multiply by T
600 Longitude Radius
601 D l' l F .1" .00001" .1km .00001km */
602
603 0, 1, 0, 0, 16, 7680, -1,-2302,
604 2,-1,-1, 0, -5,-1642, 3, 8245,
605 2,-1, 0, 0, -4,-1383, 5, 1395,
606 0, 1,-1, 0, 3, 7115, 3, 2654,
607 0, 1, 1, 0, 2, 7560, -2,-6396,
608 2, 1,-1, 0, 0, 7118, 0,-6068,
609 2, 1, 0, 0, 0, 6128, 0,-7754,
610 1, 1, 0, 0, 0,-4516, 0, 4194,
611 2,-2, 0, 0, 0,-4048, 0, 4970,
612 0, 2, 0, 0, 0, 3747, 0, -540,
613 2,-2,-1, 0, 0,-3707, 0, 2490,
614 2,-1, 1, 0, 0,-3649, 0, 3222,
615 0, 1,-2, 0, 0, 2438, 0, 1760,
616 2,-1,-2, 0, 0,-2165, 0,-2530,
617 0, 1, 2, 0, 0, 1923, 0,-1450,
618 0, 2,-1, 0, 0, 1292, 0, 1070,
619 2, 2,-1, 0, 0, 1271, 0,-6070,
620 4,-1,-1, 0, 0,-1098, 0, 990,
621 2, 0, 0, 0, 0, 1073, 0,-1360,
622 2, 0,-1, 0, 0, 839, 0, -630,
623 2, 1, 1, 0, 0, 734, 0, -660,
624 4,-1,-2, 0, 0, -688, 0, 480,
625 2, 1,-2, 0, 0, -630, 0, 0,
626 0, 2, 1, 0, 0, 587, 0, -590,
627 2,-1, 0,-2, 0, -540, 0, -170,
628 4,-1, 0, 0, 0, -468, 0, 390,
629 2,-2, 1, 0, 0, -378, 0, 330,
630 2, 1, 0,-2, 0, 364, 0, 0,
631 1, 1, 1, 0, 0, -317, 0, 240,
632 2,-1, 2, 0, 0, -295, 0, 210,
633 1, 1,-1, 0, 0, -270, 0, -210,
634 2,-3, 0, 0, 0, -256, 0, 310,
635 2,-3,-1, 0, 0, -187, 0, 110,
636 0, 1,-3, 0, 0, 169, 0, 110,
637 4, 1,-1, 0, 0, 158, 0, -150,
638 4,-2,-1, 0, 0, -155, 0, 140,
639 0, 0, 1, 0, 0, 155, 0, -250,
640 2,-2,-2, 0, 0, -148, 0, -170,
641 };
642
643 #define NBT 16
644 static const short BT[5*NBT] = {
645 /*
646 Multiply by T
647 Latitude
648 D l' l F .00001" */
649
650 2,-1, 0,-1, -7430,
651 2, 1, 0,-1, 3043,
652 2,-1,-1, 1, -2229,
653 2,-1, 0, 1, -1999,
654 2,-1,-1,-1, -1869,
655 0, 1,-1,-1, 1696,
656 0, 1, 0, 1, 1623,
657 0, 1,-1, 1, 1418,
658 0, 1, 1, 1, 1339,
659 0, 1, 1,-1, 1278,
660 0, 1, 0,-1, 1217,
661 2,-2, 0,-1, -547,
662 2,-1, 1,-1, -443,
663 2, 1,-1, 1, 331,
664 2, 1, 0, 1, 317,
665 2, 0, 0,-1, 295,
666 };
667
668 #define NLRT2 25
669 static const short LRT2[6*NLRT2] = {
670 /*
671 Multiply by T^2
672 Longitude Radius
673 D l' l F .00001" .00001km */
674
675 0, 1, 0, 0, 487, -36,
676 2,-1,-1, 0, -150, 111,
677 2,-1, 0, 0, -120, 149,
678 0, 1,-1, 0, 108, 95,
679 0, 1, 1, 0, 80, -77,
680 2, 1,-1, 0, 21, -18,
681 2, 1, 0, 0, 20, -23,
682 1, 1, 0, 0, -13, 12,
683 2,-2, 0, 0, -12, 14,
684 2,-1, 1, 0, -11, 9,
685 2,-2,-1, 0, -11, 7,
686 0, 2, 0, 0, 11, 0,
687 2,-1,-2, 0, -6, -7,
688 0, 1,-2, 0, 7, 5,
689 0, 1, 2, 0, 6, -4,
690 2, 2,-1, 0, 5, -3,
691 0, 2,-1, 0, 5, 3,
692 4,-1,-1, 0, -3, 3,
693 2, 0, 0, 0, 3, -4,
694 4,-1,-2, 0, -2, 0,
695 2, 1,-2, 0, -2, 0,
696 2,-1, 0,-2, -2, 0,
697 2, 1, 1, 0, 2, -2,
698 2, 0,-1, 0, 2, 0,
699 0, 2, 1, 0, 2, 0,
700 };
701
702 #define NBT2 12
703 static const short BT2[5*NBT2] = {
704 /*
705 Multiply by T^2
706 Latitiude
707 D l' l F .00001" */
708
709 2,-1, 0,-1, -22,
710 2, 1, 0,-1, 9,
711 2,-1, 0, 1, -6,
712 2,-1,-1, 1, -6,
713 2,-1,-1,-1, -5,
714 0, 1, 0, 1, 5,
715 0, 1,-1,-1, 5,
716 0, 1, 1, 1, 4,
717 0, 1, 1,-1, 4,
718 0, 1, 0,-1, 4,
719 0, 1,-1, 1, 4,
720 2,-2, 0,-1, -2,
721 };
722
723 /* corrections for mean lunar node in degrees, from -13100 to 17200,
724 * in 100-year steps. corrections are set to 0 between the years 0 and 3000 */
725 static const double mean_node_corr[] = {
726 -2.56,
727 -2.473, -2.392347, -2.316425, -2.239639, -2.167764, -2.095100, -2.024810, -1.957622, -1.890097, -1.826389,
728 -1.763335, -1.701047, -1.643016, -1.584186, -1.527309, -1.473352, -1.418917, -1.367736, -1.317202, -1.267269,
729 -1.221121, -1.174218, -1.128862, -1.086214, -1.042998, -1.002491, -0.962635, -0.923176, -0.887191, -0.850403,
730 -0.814929, -0.782117, -0.748462, -0.717241, -0.686598, -0.656013, -0.628726, -0.600460, -0.573219, -0.548634,
731 -0.522931, -0.499285, -0.476273, -0.452978, -0.432663, -0.411386, -0.390788, -0.372825, -0.353681, -0.336230,
732 -0.319520, -0.302343, -0.287794, -0.272262, -0.257166, -0.244534, -0.230635, -0.218126, -0.206365, -0.194000,
733 -0.183876, -0.172782, -0.161877, -0.153254, -0.143371, -0.134501, -0.126552, -0.117932, -0.111199, -0.103716,
734 -0.096160, -0.090718, -0.084046, -0.078007, -0.072959, -0.067235, -0.062990, -0.058102, -0.053070, -0.049786,
735 -0.045381, -0.041317, -0.038165, -0.034501, -0.031871, -0.028844, -0.025701, -0.024018, -0.021427, -0.018881,
736 -0.017291, -0.015186, -0.013755, -0.012098, -0.010261, -0.009688, -0.008218, -0.006670, -0.005979, -0.004756,
737 -0.003991, -0.002996, -0.001974, -0.001975, -0.001213, -0.000377, -0.000356, 5.779e-05, 0.000378, 0.000710,
738 0.001092, 0.000767, 0.000985, 0.001443, 0.001069, 0.001141, 0.001321, 0.001462, 0.001695, 0.001319,
739 0.001567, 0.001873, 0.001376, 0.001336, 0.001347, 0.001330, 0.001256, 0.000813, 0.000946, 0.001079,
740 #if 0
741 0.000509, 0.000375, 0.000477, 0.000321, 0.000279, 5.998e-05, 0.000251, 0.000623, 0.000180, 0.000225,
742 0.000506, 0.000331, 0.000253, 4.156e-05, 0.000247, 0.000394, -9.294e-05, -2.738e-05, 0.000140, -6.193e-05,
743 -0.000232, -0.000361, -0.000152, -3.571e-05, -0.000395, -0.000218, 0.000127, -0.000125, -0.000254, -0.000339,
744 #endif
745 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
746 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
747 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
748 -0.000364, -0.000452, -0.001091, -0.001159, -0.001136, -0.001798, -0.002249, -0.002622, -0.002990, -0.003555,
749 -0.004425, -0.004758, -0.005134, -0.006065, -0.006839, -0.007474, -0.008283, -0.009411, -0.010786, -0.011810,
750 -0.012989, -0.014825, -0.016426, -0.017922, -0.019774, -0.021881, -0.024194, -0.026190, -0.028440, -0.031285,
751 -0.033817, -0.036318, -0.039212, -0.042456, -0.045799, -0.048994, -0.052710, -0.056948, -0.061017, -0.065181,
752 -0.069843, -0.074922, -0.079976, -0.085052, -0.090755, -0.096840, -0.102797, -0.108939, -0.115568, -0.122636,
753 -0.129593, -0.136683, -0.144641, -0.152825, -0.161044, -0.169758, -0.178916, -0.188712, -0.198401, -0.208312,
754 -0.219395, -0.230407, -0.241577, -0.253508, -0.265640, -0.278556, -0.291330, -0.304353, -0.318815, -0.332882,
755 -0.347316, -0.362895, -0.378421, -0.395061, -0.411748, -0.428666, -0.447477, -0.465636, -0.484277, -0.504600,
756 -0.524405, -0.545533, -0.567020, -0.588404, -0.612099, -0.634965, -0.658262, -0.683866, -0.708526, -0.734719,
757 -0.761800, -0.788562, -0.818092, -0.846885, -0.876177, -0.908385, -0.939371, -0.972027, -1.006149, -1.039634,
758 -1.076135, -1.112156, -1.148490, -1.188312, -1.226761, -1.266821, -1.309156, -1.350583, -1.395223, -1.440028,
759 -1.485047, -1.534104, -1.582023, -1.631506, -1.684031, -1.735687, -1.790421, -1.846039, -1.901951, -1.961872,
760 -2.021179, -2.081987, -2.146259, -2.210031, -2.276609, -2.344904, -2.413795, -2.486559, -2.559564, -2.634215,
761 -2.712692, -2.791289, -2.872533, -2.956217, -3.040965, -3.129234, -3.218545, -3.309805, -3.404827, -3.5008,
762 -3.601, -3.7, -3.8,
763 };
764
765 /* corrections for mean lunar apsides in degrees, from -13100 to 17200,
766 * in 100-year steps. corrections are set to 0 between the years 0 and 3000 */
767 static const double mean_apsis_corr[] = {
768 7.525,
769 7.290, 7.057295, 6.830813, 6.611723, 6.396775, 6.189569, 5.985968, 5.788342, 5.597304, 5.410167,
770 5.229946, 5.053389, 4.882187, 4.716494, 4.553532, 4.396734, 4.243718, 4.094282, 3.950865, 3.810366,
771 3.674978, 3.543284, 3.414270, 3.290526, 3.168775, 3.050904, 2.937541, 2.826189, 2.719822, 2.616193,
772 2.515431, 2.419193, 2.323782, 2.232545, 2.143635, 2.056803, 1.974913, 1.893874, 1.816201, 1.741957,
773 1.668083, 1.598335, 1.529645, 1.463016, 1.399693, 1.336905, 1.278097, 1.220965, 1.165092, 1.113071,
774 1.060858, 1.011007, 0.963701, 0.916523, 0.872887, 0.829596, 0.788486, 0.750017, 0.711177, 0.675589,
775 0.640303, 0.605303, 0.573490, 0.541113, 0.511482, 0.483159, 0.455210, 0.430305, 0.404643, 0.380782,
776 0.358524, 0.335405, 0.315244, 0.295131, 0.275766, 0.259223, 0.241586, 0.225890, 0.210404, 0.194775,
777 0.181573, 0.167246, 0.154514, 0.143435, 0.131131, 0.121648, 0.111835, 0.102474, 0.094284, 0.085204,
778 0.078240, 0.070697, 0.063696, 0.058894, 0.052390, 0.047632, 0.043129, 0.037823, 0.034143, 0.029188,
779 0.025648, 0.021972, 0.018348, 0.017127, 0.013989, 0.011967, 0.011003, 0.007865, 0.007033, 0.005574,
780 0.004060, 0.003699, 0.002465, 0.002889, 0.002144, 0.001018, 0.001757, -9.67e-05, -0.000734, -0.000392,
781 -0.001546, -0.000863, -0.001266, -0.000933, -0.000503, -0.001304, 0.000238, -0.000507, -0.000897, 0.000647,
782 #if 0
783 -0.000247, 0.000938, 0.001373, 0.001159, 0.001644, 0.000691, 0.001454, 0.000532, -0.000249, 0.000871,
784 -0.000210, 0.000171, 0.000702, 0.000389, 0.000609, -0.000250, 0.000426, 0.000123, -0.000339, 0.001200,
785 0.000413, 0.000612, 0.001169, 0.000163, 0.000553, -0.000330, -0.000498, -0.000224, -0.000948, 0.000863,
786 #endif
787 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
788 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
789 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
790 0.000514, 0.000683, 0.002228, 0.001974, 0.003485, 0.004280, 0.005409, 0.007468, 0.007938, 0.011012,
791 0.012525, 0.013757, 0.016757, 0.017932, 0.020780, 0.023416, 0.026386, 0.030428, 0.033512, 0.038789,
792 0.043126, 0.047778, 0.054175, 0.058891, 0.065878, 0.072345, 0.079668, 0.088238, 0.095307, 0.104873,
793 0.113533, 0.122336, 0.133205, 0.142922, 0.154871, 0.166488, 0.179234, 0.193928, 0.207262, 0.223089,
794 0.238736, 0.254907, 0.273232, 0.291085, 0.311046, 0.331025, 0.351955, 0.374422, 0.396341, 0.420772,
795 0.444867, 0.469984, 0.497448, 0.524717, 0.554752, 0.584581, 0.616272, 0.649744, 0.682947, 0.719405,
796 0.755834, 0.793780, 0.833875, 0.873893, 0.917340, 0.960429, 1.005471, 1.052384, 1.099317, 1.149508,
797 1.200130, 1.253038, 1.307672, 1.363480, 1.422592, 1.481900, 1.544111, 1.607982, 1.672954, 1.741025,
798 1.809727, 1.882038, 1.955243, 2.029956, 2.108428, 2.186805, 2.268697, 2.352071, 2.437370, 2.525903,
799 2.615415, 2.709082, 2.804198, 2.901704, 3.002606, 3.104412, 3.210406, 3.317733, 3.428386, 3.541634,
800 3.656634, 3.775988, 3.896306, 4.020480, 4.146814, 4.275356, 4.408257, 4.542282, 4.681174, 4.822524,
801 4.966424, 5.114948, 5.264973, 5.419906, 5.577056, 5.737688, 5.902347, 6.069138, 6.241065, 6.415155,
802 6.593317, 6.774853, 6.959322, 7.148845, 7.340334, 7.537156, 7.737358, 7.940882, 8.149932, 8.361576,
803 8.579150, 8.799591, 9.024378, 9.254584, 9.487362, 9.726535, 9.968784, 10.216089, 10.467716, 10.725293,
804 10.986, 11.25, 11.52,
805 };
806
807 /* The following times are set up by update() and refer
808 * to the same instant. The distinction between them
809 * is required by altaz().
810 */
811 static TLS double ss[5][8];
812 static TLS double cc[5][8];
813
814 static TLS double l; /* Moon's ecliptic longitude */
815 static TLS double B; /* Ecliptic latitude */
816
817 static TLS double moonpol[3];
818
819 /* Orbit calculation begins.
820 */
821 static TLS double SWELP;
822 static TLS double M;
823 static TLS double MP;
824 static TLS double D;
825 static TLS double NF;
826 static TLS double T;
827 static TLS double T2;
828
829 static TLS double T3;
830 static TLS double T4;
831 static TLS double f;
832 static TLS double g;
833 static TLS double Ve;
834 static TLS double Ea;
835 static TLS double Ma;
836 static TLS double Ju;
837 static TLS double Sa;
838 static TLS double cg;
839 static TLS double sg;
840 static TLS double l1;
841 static TLS double l2;
842 static TLS double l3;
843 static TLS double l4;
844
845 /* Calculate geometric coordinates of Moon
846 * without light time or nutation correction.
847 */
swi_moshmoon2(double J,double * pol)848 int swi_moshmoon2(double J, double *pol)
849 {
850 int i;
851 T = (J-J2000)/36525.0;
852 T2 = T*T;
853 mean_elements();
854 mean_elements_pl();
855 moon1();
856 moon2();
857 moon3();
858 moon4();
859 for( i=0; i<3; i++ )
860 pol[i] = moonpol[i];
861 return(0);
862 }
863
864 /* Moshier's moom
865 * tjd julian day
866 * xpm array of 6 doubles for moon's position and speed vectors
867 * serr pointer to error string
868 */
swi_moshmoon(double tjd,AS_BOOL do_save,double * xpmret,char * serr)869 int swi_moshmoon(double tjd, AS_BOOL do_save, double *xpmret, char *serr)
870 {
871 int i;
872 double a, b, x1[6], x2[6], t;
873 double xx[6], *xpm;
874 struct plan_data *pdp = &swed.pldat[SEI_MOON];
875 char s[AS_MAXCH];
876 if (do_save)
877 xpm = pdp->x;
878 else
879 xpm = xx;
880 /* allow 0.2 day tolerance so that true node interval fits in */
881 if (tjd < MOSHLUEPH_START - 0.2 || tjd > MOSHLUEPH_END + 0.2) {
882 if (serr != NULL) {
883 sprintf(s, "jd %f outside Moshier's Moon range %.2f .. %.2f ",
884 tjd, MOSHLUEPH_START, MOSHLUEPH_END);
885 if (strlen(serr) + strlen(s) < AS_MAXCH)
886 strcat(serr, s);
887 }
888 return(ERR);
889 }
890 /* if moon has already been computed */
891 if (tjd == pdp->teval && pdp->iephe == SEFLG_MOSEPH) {
892 if (xpmret != NULL)
893 for (i = 0; i <= 5; i++)
894 xpmret[i] = pdp->x[i];
895 return(OK);
896 }
897 /* else compute moon */
898 swi_moshmoon2(tjd, xpm);
899 if (do_save) {
900 pdp->teval = tjd;
901 pdp->xflgs = -1;
902 pdp->iephe = SEFLG_MOSEPH;
903 }
904 /* Moshier moon is referred to ecliptic of date. But we need
905 * equatorial positions for several reasons.
906 * e.g. computation of earth from emb and moon
907 * of heliocentric moon
908 * Besides, this helps to keep the program structure simpler
909 */
910 ecldat_equ2000(tjd, xpm);
911 /* speed */
912 /* from 2 other positions. */
913 /* one would be good enough for computation of osculating node,
914 * but not for osculating apogee */
915 t = tjd + MOON_SPEED_INTV;
916 swi_moshmoon2(t, x1);
917 ecldat_equ2000(t, x1);
918 t = tjd - MOON_SPEED_INTV;
919 swi_moshmoon2(t, x2);
920 ecldat_equ2000(t, x2);
921 for (i = 0; i <= 2; i++) {
922 #if 0
923 xpm[i+3] = (x1[i] - x2[i]) / MOON_SPEED_INTV / 2;
924 #else
925 b = (x1[i] - x2[i]) / 2;
926 a = (x1[i] + x2[i]) / 2 - xpm[i];
927 xpm[i+3] = (2 * a + b) / MOON_SPEED_INTV;
928 #endif
929 }
930 if (xpmret != NULL)
931 for (i = 0; i <= 5; i++)
932 xpmret[i] = xpm[i];
933 return(OK);
934 }
935
936 #ifdef MOSH_MOON_200
moon1()937 static void moon1()
938 {
939 double a;
940
941 sscc( 0, STR*D, 6 );
942 sscc( 1, STR*M, 4 );
943 sscc( 2, STR*MP, 4 );
944 sscc( 3, STR*NF, 4 );
945
946 moonpol[0] = 0.0;
947 moonpol[1] = 0.0;
948 moonpol[2] = 0.0;
949
950 /* terms in T^2, scale 1.0 = 10^-5" */
951 chewm( LRT2, NLRT2, 4, 2, moonpol );
952 chewm( BT2, NBT2, 4, 4, moonpol );
953
954 f = 18 * Ve - 16 * Ea;
955
956 g = STR*(f - MP ); /* 18V - 16E - l */
957 cg = cos(g);
958 sg = sin(g);
959 l = 6.367278 * cg + 12.747036 * sg; /* t^0 */
960 l1 = 23123.70 * cg - 10570.02 * sg; /* t^1 */
961 l2 = z[24] * cg + z[25] * sg; /* t^2 */
962 l3 = z[26] * cg + z[27] * sg; /* t^3 */
963 l4 = z[28] * cg + z[29] * sg; /* t^4 */
964 moonpol[2] += 5.01 * cg + 2.72 * sg;
965
966 g = STR * (10.*Ve - 3.*Ea - MP);
967 cg = cos(g);
968 sg = sin(g);
969 l += -0.253102 * cg + 0.503359 * sg;
970 l1 += 1258.46 * cg + 707.29 * sg;
971 l2 += z[30] * cg + z[31] * sg;
972 l3 += z[32] * cg + z[33] * sg;
973 l4 += z[34] * cg + z[35] * sg;
974
975 g = STR*(8.*Ve - 13.*Ea);
976 cg = cos(g);
977 sg = sin(g);
978 l += -0.187231 * cg - 0.127481 * sg;
979 l1 += -319.87 * cg - 18.34 * sg;
980 l2 += z[36] * cg + z[37] * sg;
981 l3 += z[38] * cg + z[39] * sg;
982 l4 += z[40] * cg + z[41] * sg;
983
984 a = 4.0*Ea - 8.0*Ma + 3.0*Ju;
985 g = STR * a;
986 cg = cos(g);
987 sg = sin(g);
988 l += -0.866287 * cg + 0.248192 * sg;
989 l1 += 41.87 * cg + 1053.97 * sg;
990 l2 += z[42] * cg + z[43] * sg;
991
992 g = STR*(a - MP);
993 cg = cos(g);
994 sg = sin(g);
995 l += -0.165009 * cg + 0.044176 * sg;
996 l1 += 4.67 * cg + 201.55 * sg;
997
998
999 g = STR*f; /* 18V - 16E */
1000 cg = cos(g);
1001 sg = sin(g);
1002 l += 0.330401 * cg + 0.661362 * sg;
1003 l1 += 1202.67 * cg - 555.59 * sg;
1004 l2 += z[44] * cg + z[45] * sg;
1005 l3 += z[46] * cg + z[47] * sg;
1006
1007 g = STR*(f - 2.0*MP ); /* 18V - 16E - 2l */
1008 cg = cos(g);
1009 sg = sin(g);
1010 l += 0.352185 * cg + 0.705041 * sg;
1011 l1 += 1283.59 * cg - 586.43 * sg;
1012 l2 += z[48] * cg + z[49] * sg;
1013 l3 += z[50] * cg + z[51] * sg;
1014
1015 g = STR * (2.0*Ju - 5.0*Sa);
1016 cg = cos(g);
1017 sg = sin(g);
1018 l += -0.034700 * cg + 0.160041 * sg;
1019 l2 += z[52] * cg + z[53] * sg;
1020
1021 g = STR * (SWELP - NF);
1022 cg = cos(g);
1023 sg = sin(g);
1024 l += 0.000116 * cg + 7.063040 * sg;
1025 l1 += 298.8 * sg;
1026 l2 += z[54] * cg + z[55] * sg;
1027
1028
1029 /* T^3 terms */
1030 sg = sin( STR * M );
1031 l3 += z[56] * sg;
1032 l4 += z[57] * sg;
1033
1034 g = STR * (2.0*D - M);
1035 sg = sin(g);
1036 cg = cos(g);
1037 l3 += z[58] * sg;
1038 l4 += z[59] * sg;
1039 moonpol[2] += -0.2655 * cg * T;
1040
1041 g = g - STR * MP;
1042 sg = sin(g);
1043 l3 += z[60] * sg;
1044 l4 += z[61] * sg;
1045
1046 g = STR * (M - MP);
1047 l3 += z[62] * sin( g );
1048 moonpol[2] += -0.1568 * cos( g ) * T;
1049
1050 g = STR * (M + MP);
1051 l3 += z[63] * sin( g );
1052 moonpol[2] += 0.1309 * cos( g ) * T;
1053
1054 g = STR * 2.0 * (D - M);
1055 sg = sin(g);
1056 l3 += z[64] * sg;
1057 l4 += z[65] * sg;
1058
1059 g = STR * 2.0 * M;
1060 sg = sin(g);
1061 l3 += z[66] * sg;
1062 l4 += z[67] * sg;
1063
1064 g = STR * (2.0*D - MP);
1065 sg = sin(g);
1066 l3 += z[68] * sg;
1067
1068 g = STR * (2.0*(D - M) - MP);
1069 sg = sin(g);
1070 l3 += z[69] * sg;
1071
1072 g = STR * (2.0*(D + M) - MP);
1073 sg = sin(g);
1074 cg = cos(g);
1075 l3 += z[70] * sg;
1076 moonpol[2] += 0.5568 * cg * T;
1077
1078 l2 += moonpol[0];
1079
1080 g = STR*(2.0*D - M - MP);
1081 moonpol[2] += -0.1910 * cos( g ) * T;
1082
1083
1084 moonpol[1] *= T;
1085 moonpol[2] *= T;
1086
1087 /* terms in T */
1088 moonpol[0] = 0.0;
1089 chewm( BT, NBT, 4, 4, moonpol );
1090 chewm( LRT, NLRT, 4, 1, moonpol );
1091 g = STR*(f - MP - NF - 2355767.6); /* 18V - 16E - l - F */
1092 moonpol[1] += -1127. * sin(g);
1093 g = STR*(f - MP + NF - 235353.6); /* 18V - 16E - l + F */
1094 moonpol[1] += -1123. * sin(g);
1095 g = STR*(Ea + D + 51987.6);
1096 moonpol[1] += 1303. * sin(g);
1097 g = STR*SWELP;
1098 moonpol[1] += 342. * sin(g);
1099
1100
1101 g = STR*(2.*Ve - 3.*Ea);
1102 cg = cos(g);
1103 sg = sin(g);
1104 l += -0.343550 * cg - 0.000276 * sg;
1105 l1 += 105.90 * cg + 336.53 * sg;
1106
1107 g = STR*(f - 2.*D); /* 18V - 16E - 2D */
1108 cg = cos(g);
1109 sg = sin(g);
1110 l += 0.074668 * cg + 0.149501 * sg;
1111 l1 += 271.77 * cg - 124.20 * sg;
1112
1113 g = STR*(f - 2.*D - MP);
1114 cg = cos(g);
1115 sg = sin(g);
1116 l += 0.073444 * cg + 0.147094 * sg;
1117 l1 += 265.24 * cg - 121.16 * sg;
1118
1119 g = STR*(f + 2.*D - MP);
1120 cg = cos(g);
1121 sg = sin(g);
1122 l += 0.072844 * cg + 0.145829 * sg;
1123 l1 += 265.18 * cg - 121.29 * sg;
1124
1125 g = STR*(f + 2.*(D - MP));
1126 cg = cos(g);
1127 sg = sin(g);
1128 l += 0.070201 * cg + 0.140542 * sg;
1129 l1 += 255.36 * cg - 116.79 * sg;
1130
1131 g = STR*(Ea + D - NF);
1132 cg = cos(g);
1133 sg = sin(g);
1134 l += 0.288209 * cg - 0.025901 * sg;
1135 l1 += -63.51 * cg - 240.14 * sg;
1136
1137 g = STR*(2.*Ea - 3.*Ju + 2.*D - MP);
1138 cg = cos(g);
1139 sg = sin(g);
1140 l += 0.077865 * cg + 0.438460 * sg;
1141 l1 += 210.57 * cg + 124.84 * sg;
1142
1143 g = STR*(Ea - 2.*Ma);
1144 cg = cos(g);
1145 sg = sin(g);
1146 l += -0.216579 * cg + 0.241702 * sg;
1147 l1 += 197.67 * cg + 125.23 * sg;
1148
1149 g = STR*(a + MP);
1150 cg = cos(g);
1151 sg = sin(g);
1152 l += -0.165009 * cg + 0.044176 * sg;
1153 l1 += 4.67 * cg + 201.55 * sg;
1154
1155 g = STR*(a + 2.*D - MP);
1156 cg = cos(g);
1157 sg = sin(g);
1158 l += -0.133533 * cg + 0.041116 * sg;
1159 l1 += 6.95 * cg + 187.07 * sg;
1160
1161 g = STR*(a - 2.*D + MP);
1162 cg = cos(g);
1163 sg = sin(g);
1164 l += -0.133430 * cg + 0.041079 * sg;
1165 l1 += 6.28 * cg + 169.08 * sg;
1166
1167 g = STR*(3.*Ve - 4.*Ea);
1168 cg = cos(g);
1169 sg = sin(g);
1170 l += -0.175074 * cg + 0.003035 * sg;
1171 l1 += 49.17 * cg + 150.57 * sg;
1172
1173 g = STR*(2.*(Ea + D - MP) - 3.*Ju + 213534.);
1174 l1 += 158.4 * sin(g);
1175 l1 += moonpol[0];
1176
1177 a = 0.1 * T; /* set amplitude scale of 1.0 = 10^-4 arcsec */
1178 moonpol[1] *= a;
1179 moonpol[2] *= a;
1180 }
1181 #else
moon1()1182 static void moon1()
1183 {
1184 double a;
1185 /* This code added by Bhanu Pinnamaneni, 17-aug-2009 */
1186 /* Note by Dieter: Bhanu noted that ss and cc are not sufficiently
1187 * initialised and random values are used for the calculation.
1188 * However, this may be only part of the bug.
1189 * The bug could be in sscc(). Or may be the bug is rather in
1190 * the 116th line of NLR, where the value "5" may be wrong.
1191 * Still, this will make a maximum difference of only 0.1", while the error
1192 * of the Moshier lunar ephemeris can reach 7". */
1193 int i, j;
1194 for (i = 0; i < 5; i++) {
1195 for (j = 0; j < 8; j++) {
1196 ss[i][j] = 0;
1197 cc[i][j] = 0;
1198 }
1199 }
1200 /* End of code addition */
1201 sscc( 0, STR*D, 6 );
1202 sscc( 1, STR*M, 4 );
1203 sscc( 2, STR*MP, 4 );
1204 sscc( 3, STR*NF, 4 );
1205 moonpol[0] = 0.0;
1206 moonpol[1] = 0.0;
1207 moonpol[2] = 0.0;
1208 /* terms in T^2, scale 1.0 = 10^-5" */
1209 chewm( LRT2, NLRT2, 4, 2, moonpol );
1210 chewm( BT2, NBT2, 4, 4, moonpol );
1211 f = 18 * Ve - 16 * Ea;
1212 g = STR*(f - MP ); /* 18V - 16E - l */
1213 cg = cos(g);
1214 sg = sin(g);
1215 l = 6.367278 * cg + 12.747036 * sg; /* t^0 */
1216 l1 = 23123.70 * cg - 10570.02 * sg; /* t^1 */
1217 l2 = z[12] * cg + z[13] * sg; /* t^2 */
1218 moonpol[2] += 5.01 * cg + 2.72 * sg;
1219 g = STR * (10.*Ve - 3.*Ea - MP);
1220 cg = cos(g);
1221 sg = sin(g);
1222 l += -0.253102 * cg + 0.503359 * sg;
1223 l1 += 1258.46 * cg + 707.29 * sg;
1224 l2 += z[14] * cg + z[15] * sg;
1225 g = STR*(8.*Ve - 13.*Ea);
1226 cg = cos(g);
1227 sg = sin(g);
1228 l += -0.187231 * cg - 0.127481 * sg;
1229 l1 += -319.87 * cg - 18.34 * sg;
1230 l2 += z[16] * cg + z[17] * sg;
1231 a = 4.0*Ea - 8.0*Ma + 3.0*Ju;
1232 g = STR * a;
1233 cg = cos(g);
1234 sg = sin(g);
1235 l += -0.866287 * cg + 0.248192 * sg;
1236 l1 += 41.87 * cg + 1053.97 * sg;
1237 l2 += z[18] * cg + z[19] * sg;
1238 g = STR*(a - MP);
1239 cg = cos(g);
1240 sg = sin(g);
1241 l += -0.165009 * cg + 0.044176 * sg;
1242 l1 += 4.67 * cg + 201.55 * sg;
1243 g = STR*f; /* 18V - 16E */
1244 cg = cos(g);
1245 sg = sin(g);
1246 l += 0.330401 * cg + 0.661362 * sg;
1247 l1 += 1202.67 * cg - 555.59 * sg;
1248 l2 += z[20] * cg + z[21] * sg;
1249 g = STR*(f - 2.0*MP ); /* 18V - 16E - 2l */
1250 cg = cos(g);
1251 sg = sin(g);
1252 l += 0.352185 * cg + 0.705041 * sg;
1253 l1 += 1283.59 * cg - 586.43 * sg;
1254 g = STR * (2.0*Ju - 5.0*Sa);
1255 cg = cos(g);
1256 sg = sin(g);
1257 l += -0.034700 * cg + 0.160041 * sg;
1258 l2 += z[22] * cg + z[23] * sg;
1259 g = STR * (SWELP - NF);
1260 cg = cos(g);
1261 sg = sin(g);
1262 l += 0.000116 * cg + 7.063040 * sg;
1263 l1 += 298.8 * sg;
1264 /* T^3 terms */
1265 sg = sin( STR * M );
1266 /* l3 += z[24] * sg; moshier! l3 not initialized! */
1267 l3 = z[24] * sg;
1268 l4 = 0;
1269 g = STR * (2.0*D - M);
1270 sg = sin(g);
1271 cg = cos(g);
1272 moonpol[2] += -0.2655 * cg * T;
1273 g = STR * (M - MP);
1274 moonpol[2] += -0.1568 * cos( g ) * T;
1275 g = STR * (M + MP);
1276 moonpol[2] += 0.1309 * cos( g ) * T;
1277 g = STR * (2.0*(D + M) - MP);
1278 sg = sin(g);
1279 cg = cos(g);
1280 moonpol[2] += 0.5568 * cg * T;
1281 l2 += moonpol[0];
1282 g = STR*(2.0*D - M - MP);
1283 moonpol[2] += -0.1910 * cos( g ) * T;
1284 moonpol[1] *= T;
1285 moonpol[2] *= T;
1286 /* terms in T */
1287 moonpol[0] = 0.0;
1288 chewm( BT, NBT, 4, 4, moonpol );
1289 chewm( LRT, NLRT, 4, 1, moonpol );
1290 g = STR*(f - MP - NF - 2355767.6); /* 18V - 16E - l - F */
1291 moonpol[1] += -1127. * sin(g);
1292 g = STR*(f - MP + NF - 235353.6); /* 18V - 16E - l + F */
1293 moonpol[1] += -1123. * sin(g);
1294 g = STR*(Ea + D + 51987.6);
1295 moonpol[1] += 1303. * sin(g);
1296 g = STR*SWELP;
1297 moonpol[1] += 342. * sin(g);
1298 g = STR*(2.*Ve - 3.*Ea);
1299 cg = cos(g);
1300 sg = sin(g);
1301 l += -0.343550 * cg - 0.000276 * sg;
1302 l1 += 105.90 * cg + 336.53 * sg;
1303 g = STR*(f - 2.*D); /* 18V - 16E - 2D */
1304 cg = cos(g);
1305 sg = sin(g);
1306 l += 0.074668 * cg + 0.149501 * sg;
1307 l1 += 271.77 * cg - 124.20 * sg;
1308 g = STR*(f - 2.*D - MP);
1309 cg = cos(g);
1310 sg = sin(g);
1311 l += 0.073444 * cg + 0.147094 * sg;
1312 l1 += 265.24 * cg - 121.16 * sg;
1313 g = STR*(f + 2.*D - MP);
1314 cg = cos(g);
1315 sg = sin(g);
1316 l += 0.072844 * cg + 0.145829 * sg;
1317 l1 += 265.18 * cg - 121.29 * sg;
1318 g = STR*(f + 2.*(D - MP));
1319 cg = cos(g);
1320 sg = sin(g);
1321 l += 0.070201 * cg + 0.140542 * sg;
1322 l1 += 255.36 * cg - 116.79 * sg;
1323 g = STR*(Ea + D - NF);
1324 cg = cos(g);
1325 sg = sin(g);
1326 l += 0.288209 * cg - 0.025901 * sg;
1327 l1 += -63.51 * cg - 240.14 * sg;
1328 g = STR*(2.*Ea - 3.*Ju + 2.*D - MP);
1329 cg = cos(g);
1330 sg = sin(g);
1331 l += 0.077865 * cg + 0.438460 * sg;
1332 l1 += 210.57 * cg + 124.84 * sg;
1333 g = STR*(Ea - 2.*Ma);
1334 cg = cos(g);
1335 sg = sin(g);
1336 l += -0.216579 * cg + 0.241702 * sg;
1337 l1 += 197.67 * cg + 125.23 * sg;
1338 g = STR*(a + MP);
1339 cg = cos(g);
1340 sg = sin(g);
1341 l += -0.165009 * cg + 0.044176 * sg;
1342 l1 += 4.67 * cg + 201.55 * sg;
1343 g = STR*(a + 2.*D - MP);
1344 cg = cos(g);
1345 sg = sin(g);
1346 l += -0.133533 * cg + 0.041116 * sg;
1347 l1 += 6.95 * cg + 187.07 * sg;
1348 g = STR*(a - 2.*D + MP);
1349 cg = cos(g);
1350 sg = sin(g);
1351 l += -0.133430 * cg + 0.041079 * sg;
1352 l1 += 6.28 * cg + 169.08 * sg;
1353 g = STR*(3.*Ve - 4.*Ea);
1354 cg = cos(g);
1355 sg = sin(g);
1356 l += -0.175074 * cg + 0.003035 * sg;
1357 l1 += 49.17 * cg + 150.57 * sg;
1358 g = STR*(2.*(Ea + D - MP) - 3.*Ju + 213534.);
1359 l1 += 158.4 * sin(g);
1360 l1 += moonpol[0];
1361 a = 0.1 * T; /* set amplitude scale of 1.0 = 10^-4 arcsec */
1362 moonpol[1] *= a;
1363 moonpol[2] *= a;
1364 }
1365 #endif /* MOSH_MOON_200 */
1366
moon2()1367 static void moon2()
1368 {
1369 /* terms in T^0 */
1370 g = STR*(2*(Ea-Ju+D)-MP+648431.172);
1371 l += 1.14307 * sin(g);
1372 g = STR*(Ve-Ea+648035.568);
1373 l += 0.82155 * sin(g);
1374 g = STR*(3*(Ve-Ea)+2*D-MP+647933.184);
1375 l += 0.64371 * sin(g);
1376 g = STR*(Ea-Ju+4424.04);
1377 l += 0.63880 * sin(g);
1378 g = STR*(SWELP + MP - NF + 4.68);
1379 l += 0.49331 * sin(g);
1380 g = STR*(SWELP - MP - NF + 4.68);
1381 l += 0.4914 * sin(g);
1382 g = STR*(SWELP+NF+2.52);
1383 l += 0.36061 * sin(g);
1384 g = STR*(2.*Ve - 2.*Ea + 736.2);
1385 l += 0.30154 * sin(g);
1386 g = STR*(2.*Ea - 3.*Ju + 2.*D - 2.*MP + 36138.2);
1387 l += 0.28282 * sin(g);
1388 g = STR*(2.*Ea - 2.*Ju + 2.*D - 2.*MP + 311.0);
1389 l += 0.24516 * sin(g);
1390 g = STR*(Ea - Ju - 2.*D + MP + 6275.88);
1391 l += 0.21117 * sin(g);
1392 g = STR*(2.*(Ea - Ma) - 846.36);
1393 l += 0.19444 * sin(g);
1394 g = STR*(2.*(Ea - Ju) + 1569.96);
1395 l -= 0.18457 * sin(g);
1396 g = STR*(2.*(Ea - Ju) - MP - 55.8);
1397 l += 0.18256 * sin(g);
1398 g = STR*(Ea - Ju - 2.*D + 6490.08);
1399 l += 0.16499 * sin(g);
1400 g = STR*(Ea - 2.*Ju - 212378.4);
1401 l += 0.16427 * sin(g);
1402 g = STR*(2.*(Ve - Ea - D) + MP + 1122.48);
1403 l += 0.16088 * sin(g);
1404 g = STR*(Ve - Ea - MP + 32.04);
1405 l -= 0.15350 * sin(g);
1406 g = STR*(Ea - Ju - MP + 4488.88);
1407 l += 0.14346 * sin(g);
1408 g = STR*(2.*(Ve - Ea + D) - MP - 8.64);
1409 l += 0.13594 * sin(g);
1410 g = STR*(2.*(Ve - Ea - D) + 1319.76);
1411 l += 0.13432 * sin(g);
1412 g = STR*(Ve - Ea - 2.*D + MP - 56.16);
1413 l -= 0.13122 * sin(g);
1414 g = STR*(Ve - Ea + MP + 54.36);
1415 l -= 0.12722 * sin(g);
1416 g = STR*(3.*(Ve - Ea) - MP + 433.8);
1417 l += 0.12539 * sin(g);
1418 g = STR*(Ea - Ju + MP + 4002.12);
1419 l += 0.10994 * sin(g);
1420 g = STR*(20.*Ve - 21.*Ea - 2.*D + MP - 317511.72);
1421 l += 0.10652 * sin(g);
1422 g = STR*(26.*Ve - 29.*Ea - MP + 270002.52);
1423 l += 0.10490 * sin(g);
1424 g = STR*(3.*Ve - 4.*Ea + D - MP - 322765.56);
1425 l += 0.10386 * sin(g);
1426 g = STR*(SWELP+648002.556);
1427 B = 8.04508 * sin(g);
1428 g = STR*(Ea+D+996048.252);
1429 B += 1.51021 * sin(g);
1430 g = STR*(f - MP + NF + 95554.332);
1431 B += 0.63037 * sin(g);
1432 g = STR*(f - MP - NF + 95553.792);
1433 B += 0.63014 * sin(g);
1434 g = STR*(SWELP - MP + 2.9);
1435 B += 0.45587 * sin(g);
1436 g = STR*(SWELP + MP + 2.5);
1437 B += -0.41573 * sin(g);
1438 g = STR*(SWELP - 2.0*NF + 3.2);
1439 B += 0.32623 * sin(g);
1440 g = STR*(SWELP - 2.0*D + 2.5);
1441 B += 0.29855 * sin(g);
1442 }
1443
moon3()1444 static void moon3()
1445 {
1446 /* terms in T^0 */
1447 moonpol[0] = 0.0;
1448 chewm( LR, NLR, 4, 1, moonpol );
1449 chewm( MB, NMB, 4, 3, moonpol );
1450 l += (((l4 * T + l3) * T + l2) * T + l1) * T * 1.0e-5;
1451 moonpol[0] = SWELP + l + 1.0e-4 * moonpol[0];
1452 moonpol[1] = 1.0e-4 * moonpol[1] + B;
1453 moonpol[2] = 1.0e-4 * moonpol[2] + 385000.52899; /* kilometers */
1454 }
1455
1456 /* Compute final ecliptic polar coordinates
1457 */
moon4()1458 static void moon4()
1459 {
1460 moonpol[2] /= AUNIT / 1000;
1461 moonpol[0] = STR * mods3600( moonpol[0] );
1462 moonpol[1] = STR * moonpol[1];
1463 B = moonpol[1];
1464 }
1465
1466 #define CORR_MNODE_JD_T0GREG -3063616.5 /* 1 jan -13100 greg. */
1467 #define CORR_MNODE_JD_T1GREG 844477.5 /* 1 jan -2400 jul. */
1468 #define CORR_MNODE_JD_T2GREG 2780263.5 /* 1 jan 2900 jul. */
1469 #define CORR_MNODE_JD_T3GREG 7930182.5 /* 1 jan 17000 greg. */
corr_mean_node(double J)1470 static double corr_mean_node(double J)
1471 {
1472 double J0, dJ, dayscty, dcor, dcor0, dcor1, dfrac;
1473 int i;
1474 J0 = CORR_MNODE_JD_T0GREG; /* 1-jan--13000 greg */
1475 dayscty = 36524.25; /* days per Gregorian century */
1476 if (J < JPL_DE431_START) return 0;
1477 if (J > JPL_DE431_END) return 0;
1478 /*if (J > CORR_MNODE_JD_T1GREG && J < CORR_MNODE_JD_T2GREG) return 0;*/
1479 dJ = J - J0;
1480 i = (int) floor(dJ / dayscty); /* centuries = index of lower correction value */
1481 dfrac = (dJ - i * dayscty) / dayscty;
1482 dcor0 = mean_node_corr[i];
1483 dcor1 = mean_node_corr[i + 1];
1484 dcor = dcor0 + dfrac * (dcor1 - dcor0);
1485 return dcor;
1486 }
1487
1488 /* mean lunar node
1489 * J julian day
1490 * pol return array for position and velocity
1491 * (polar coordinates of ecliptic of date)
1492 */
swi_mean_node(double J,double * pol,char * serr)1493 int swi_mean_node(double J, double *pol, char *serr)
1494 {
1495 #if 0
1496 double a, b, c;
1497 #endif
1498 char s[AS_MAXCH];
1499 double dcor;
1500 T = (J-J2000)/36525.0;
1501 T2 = T*T;
1502 T3 = T*T2;
1503 T4 = T2*T2;
1504 /* with elements from swi_moshmoon2(), which are fitted to jpl-ephemeris */
1505 if (J < MOSHNDEPH_START || J > MOSHNDEPH_END) {
1506 if (serr != NULL) {
1507 sprintf(s, "jd %f outside mean node range %.2f .. %.2f ",
1508 J, MOSHNDEPH_START, MOSHNDEPH_END);
1509 if (strlen(serr) + strlen(s) < AS_MAXCH)
1510 strcat(serr, s);
1511 }
1512 return ERR;
1513 }
1514 mean_elements();
1515 dcor = corr_mean_node(J) * 3600;
1516 /* longitude */
1517 pol[0] = swi_mod2PI((SWELP - NF - dcor) * STR);
1518 /* latitude */
1519 pol[1] = 0.0;
1520 /* distance */
1521 pol[2] = MOON_MEAN_DIST / AUNIT; /* or should it be derived from mean
1522 * orbital ellipse? */
1523 #if 0
1524 a = pol[0];
1525 /* Chapront, according to Meeus, German, p. 339 */
1526 pol[0] = 125.0445550 - 1934.1361849 * T + 0.0020762 * T2 +
1527 T3 / 467410 - T4 / 60616000;
1528 pol[0] = swi_mod2PI(pol[0] * DEGTORAD);
1529 c = pol[0];
1530 printf ("mean node\n");
1531 printf ("moshier de404 - chapront %f\"\n", (a-c) * RADTODEG * 3600);
1532 #endif
1533 return OK;
1534 }
1535
1536 #define CORR_MAPOG_JD_T0GREG -3063616.5 /* 1 jan -13100 greg. */
1537 #define CORR_MAPOG_JD_T1GREG 1209720.5 /* 1 jan -1400 greg. */
1538 #define CORR_MAPOG_JD_T2GREG 2780263.5 /* 1 jan 2900 greg. */
1539 #define CORR_MAPOG_JD_T3GREG 7930182.5 /* 1 jan 17000 greg. */
corr_mean_apog(double J)1540 static double corr_mean_apog(double J)
1541 {
1542 double J0, dJ, dayscty, dcor, dcor0, dcor1, dfrac;
1543 int i;
1544 J0 = CORR_MAPOG_JD_T0GREG; /* 1-jan--13000 greg */
1545 dayscty = 36524.25; /* days per Gregorian century */
1546 if (J < JPL_DE431_START) return 0;
1547 if (J > JPL_DE431_END) return 0;
1548 /*if (J > CORR_MAPOG_JD_T1GREG && J < CORR_MAPOG_JD_T2GREG) return 0;*/
1549 dJ = J - J0;
1550 i = (int) floor(dJ / dayscty); /* centuries = index of lower correction value */
1551 dfrac = (dJ - i * dayscty) / dayscty;
1552 dcor0 = mean_apsis_corr[i];
1553 dcor1 = mean_apsis_corr[i + 1];
1554 dcor = dcor0 + dfrac * (dcor1 - dcor0);
1555 return dcor;
1556 }
1557
1558 /* mean lunar apogee ('dark moon', 'lilith')
1559 * J julian day
1560 * pol return array for position
1561 * (polar coordinates of ecliptic of date)
1562 * serr error return string
1563 */
swi_mean_apog(double J,double * pol,char * serr)1564 int swi_mean_apog(double J, double *pol, char *serr)
1565 {
1566 #if 0
1567 int i;
1568 double a, b;
1569 double x[3];
1570 #endif
1571 double node, dcor;
1572 char s[AS_MAXCH];
1573 T = (J-J2000)/36525.0;
1574 T2 = T*T;
1575 T3 = T*T2;
1576 T4 = T2*T2;
1577 /* with elements from swi_moshmoon2(), which are fitted to jpl-ephemeris */
1578 if (J < MOSHNDEPH_START || J > MOSHNDEPH_END) {
1579 if (serr != NULL) {
1580 sprintf(s, "jd %f outside mean apogee range %.2f .. %.2f ",
1581 J, MOSHNDEPH_START, MOSHNDEPH_END);
1582 if (strlen(serr) + strlen(s) < AS_MAXCH)
1583 strcat(serr, s);
1584 }
1585 return(ERR);
1586 }
1587 mean_elements();
1588 pol[0] = swi_mod2PI((SWELP - MP) * STR + PI);
1589 pol[1] = 0;
1590 pol[2] = MOON_MEAN_DIST * (1 + MOON_MEAN_ECC) / AUNIT; /* apogee */
1591 /* Lilith or Dark Moon is either the empty focal point of the mean
1592 * lunar ellipse or, for some people, its apogee ("aphelion").
1593 * This is 180 degrees from the perigee.
1594 *
1595 * Since the lunar orbit is not in the ecliptic, the apogee must be
1596 * projected onto the ecliptic.
1597 * Joelle de Gravelaine has in her book "Lilith der schwarze Mond"
1598 * (Astrodata, 1990) an ephemeris which gives noon (12.00) positions
1599 * but does not project them onto the ecliptic.
1600 * This results in a mistake of several arc minutes.
1601 *
1602 * There is also another problem. The other focal point doesn't
1603 * coincide with the geocenter but with the barycenter of the
1604 * earth-moon-system. The difference is about 4700 km. If one
1605 * took this into account, it would result in an oscillation
1606 * of the Black Moon. If defined as the apogee, this oscillation
1607 * would be about +/- 40 arcmin.
1608 * If defined as the second focus, the effect is very large:
1609 * +/- 6 deg!
1610 * We neglect this influence.
1611 */
1612 dcor = corr_mean_apog(J) * DEGTORAD;
1613 pol[0] = swi_mod2PI(pol[0] - dcor);
1614 /* apogee is now projected onto ecliptic */
1615 node = (SWELP - NF) * STR;
1616 dcor = corr_mean_node(J) * DEGTORAD;
1617 node = swi_mod2PI(node - dcor);
1618 pol[0] = swi_mod2PI(pol[0] - node);
1619 swi_polcart(pol, pol);
1620 swi_coortrf(pol, pol, -MOON_MEAN_INCL * DEGTORAD);
1621 swi_cartpol(pol, pol);
1622 pol[0] = swi_mod2PI(pol[0] + node);
1623 return OK;
1624 }
1625
1626 /* Program to step through the perturbation table
1627 */
chewm(const short * pt,int nlines,int nangles,int typflg,double * ans)1628 static void chewm(const short *pt, int nlines, int nangles, int typflg, double *ans )
1629 {
1630 int i, j, k, k1, m;
1631 double cu, su, cv, sv, ff;
1632 for( i=0; i<nlines; i++ ) {
1633 k1 = 0;
1634 sv = 0.0;
1635 cv = 0.0;
1636 for( m=0; m<nangles; m++ ) {
1637 j = *pt++; /* multiple angle factor */
1638 if( j ) {
1639 k = j;
1640 if( j < 0 ) k = -k; /* make angle factor > 0 */
1641 /* sin, cos (k*angle) from lookup table */
1642 su = ss[m][k-1];
1643 cu = cc[m][k-1];
1644 if( j < 0 ) su = -su; /* negative angle factor */
1645 if( k1 == 0 ) {
1646 /* Set sin, cos of first angle. */
1647 sv = su;
1648 cv = cu;
1649 k1 = 1;
1650 }
1651 else {
1652 /* Combine angles by trigonometry. */
1653 ff = su*cv + cu*sv;
1654 cv = cu*cv - su*sv;
1655 sv = ff;
1656 }
1657 }
1658 }
1659 /* Accumulate
1660 */
1661 switch( typflg ) {
1662 /* large longitude and radius */
1663 case 1:
1664 j = *pt++;
1665 k = *pt++;
1666 ans[0] += (10000.0 * j + k) * sv;
1667 j = *pt++;
1668 k = *pt++;
1669 if( k ) ans[2] += (10000.0 * j + k) * cv;
1670 break;
1671 /* longitude and radius */
1672 case 2:
1673 j = *pt++;
1674 k = *pt++;
1675 ans[0] += j * sv;
1676 ans[2] += k * cv;
1677 break;
1678 /* large latitude */
1679 case 3:
1680 j = *pt++;
1681 k = *pt++;
1682 ans[1] += ( 10000.0*j + k)*sv;
1683 break;
1684 /* latitude */
1685 case 4:
1686 j = *pt++;
1687 ans[1] += j * sv;
1688 break;
1689 }
1690 }
1691 }
1692
1693 /* Prepare lookup table of sin and cos ( i*Lj )
1694 * for required multiple angles
1695 */
sscc(int k,double arg,int n)1696 static void sscc(int k, double arg, int n )
1697 {
1698 double cu, su, cv, sv, s;
1699 int i;
1700 su = sin(arg);
1701 cu = cos(arg);
1702 ss[k][0] = su; /* sin(L) */
1703 cc[k][0] = cu; /* cos(L) */
1704 sv = 2.0*su*cu;
1705 cv = cu*cu - su*su;
1706 ss[k][1] = sv; /* sin(2L) */
1707 cc[k][1] = cv;
1708 for( i=2; i<n; i++ ) {
1709 s = su*cv + cu*sv;
1710 cv = cu*cv - su*sv;
1711 sv = s;
1712 ss[k][i] = sv; /* sin( i+1 L ) */
1713 cc[k][i] = cv;
1714 }
1715 }
1716
1717 /* converts from polar coordinates of ecliptic of date
1718 * to cartesian coordinates of equator 2000
1719 * tjd date
1720 * x array of position
1721 */
ecldat_equ2000(double tjd,double * xpm)1722 static void ecldat_equ2000(double tjd, double *xpm) {
1723 /* cartesian */
1724 swi_polcart(xpm, xpm);
1725 /* equatorial */
1726 swi_coortrf2(xpm, xpm, -swed.oec.seps, swed.oec.ceps);
1727 /* j2000 */
1728 swi_precess(xpm, tjd, 0, J_TO_J2000);/**/
1729 }
1730
1731 /* Reduce arc seconds modulo 360 degrees
1732 * answer in arc seconds
1733 */
mods3600(double x)1734 static double mods3600(double x)
1735 {
1736 double lx;
1737 lx = x;
1738 lx = lx - 1296000.0 * floor( lx/1296000.0 );
1739 return( lx );
1740 }
1741
swi_mean_lunar_elements(double tjd,double * node,double * dnode,double * peri,double * dperi)1742 void swi_mean_lunar_elements(double tjd,
1743 double *node, double *dnode,
1744 double *peri, double *dperi)
1745 {
1746 double dcor;
1747 T = (tjd - J2000) / 36525.0;
1748 T2 = T*T;
1749 mean_elements();
1750 *node = swe_degnorm((SWELP - NF) * STR * RADTODEG);
1751 *peri = swe_degnorm((SWELP - MP) * STR * RADTODEG);
1752 T -= 1.0 / 36525;
1753 mean_elements();
1754 *dnode = swe_degnorm(*node - (SWELP - NF) * STR * RADTODEG);
1755 *dnode -= 360;
1756 *dperi = swe_degnorm(*peri - (SWELP - MP) * STR * RADTODEG);
1757 dcor = corr_mean_node(tjd);
1758 *node = swe_degnorm(*node - dcor);
1759 dcor = corr_mean_apog(tjd);
1760 *peri = swe_degnorm(*peri - dcor);
1761 }
1762
mean_elements()1763 static void mean_elements()
1764 {
1765 double fracT = fmod(T, 1);
1766 /* Mean anomaly of sun = l' (J. Laskar) */
1767 /*M = mods3600(129596581.038354 * T + 1287104.76154);*/
1768 M = mods3600(129600000.0 * fracT - 3418.961646 * T + 1287104.76154);
1769 M += ((((((((
1770 1.62e-20 * T
1771 - 1.0390e-17 ) * T
1772 - 3.83508e-15 ) * T
1773 + 4.237343e-13 ) * T
1774 + 8.8555011e-11 ) * T
1775 - 4.77258489e-8 ) * T
1776 - 1.1297037031e-5 ) * T
1777 + 1.4732069041e-4 ) * T
1778 - 0.552891801772 ) * T2;
1779 #ifdef MOSH_MOON_200
1780 /* Mean distance of moon from its ascending node = F */
1781 NF = mods3600( 1739527263.0983 * T + 335779.55755 );
1782 /* Mean anomaly of moon = l */
1783 MP = mods3600( 1717915923.4728 * T + 485868.28096 );
1784 /* Mean elongation of moon = D */
1785 D = mods3600( 1602961601.4603 * T + 1072260.73512 );
1786 /* Mean longitude of moon */
1787 SWELP = mods3600( 1732564372.83264 * T + 785939.95571 );
1788 /* Higher degree secular terms found by least squares fit */
1789 NF += (((((z[5] *T+z[4] )*T + z[3] )*T + z[2] )*T + z[1] )*T + z[0] )*T2;
1790 MP += (((((z[11]*T+z[10])*T + z[9] )*T + z[8] )*T + z[7] )*T + z[6] )*T2;
1791 D += (((((z[17]*T+z[16])*T + z[15])*T + z[14])*T + z[13])*T + z[12])*T2;
1792 SWELP += (((((z[23]*T+z[22])*T + z[21])*T + z[20])*T + z[19])*T + z[18])*T2;
1793 #else
1794 /* Mean distance of moon from its ascending node = F */
1795 /*NF = mods3600((1739527263.0983 - 2.079419901760e-01) * T + 335779.55755);*/
1796 NF = mods3600(1739232000.0 * fracT + 295263.0983 * T - 2.079419901760e-01 * T + 335779.55755);
1797 /* Mean anomaly of moon = l */
1798 /*MP = mods3600((1717915923.4728 - 2.035946368532e-01) * T + 485868.28096);*/
1799 MP = mods3600(1717200000.0 * fracT + 715923.4728 * T - 2.035946368532e-01 * T + 485868.28096);
1800 /* Mean elongation of moon = D */
1801 /*D = mods3600((1602961601.4603 + 3.962893294503e-01) * T + 1072260.73512);*/
1802 D = mods3600(1601856000.0 * fracT + 1105601.4603 * T + 3.962893294503e-01 * T + 1072260.73512);
1803 /* Mean longitude of moon, referred to the mean ecliptic and equinox of date */
1804 /*SWELP = mods3600((1732564372.83264 - 6.784914260953e-01) * T + 785939.95571);*/
1805 SWELP = mods3600(1731456000.0 * fracT + 1108372.83264 * T - 6.784914260953e-01 * T + 785939.95571);
1806 /* Higher degree secular terms found by least squares fit */
1807 NF += ((z[2]*T + z[1])*T + z[0])*T2;
1808 MP += ((z[5]*T + z[4])*T + z[3])*T2;
1809 D += ((z[8]*T + z[7])*T + z[6])*T2;
1810 SWELP += ((z[11]*T + z[10])*T + z[9])*T2;
1811 #endif /* ! MOSH_MOON_200 */
1812 /* sensitivity of mean elements
1813 * delta argument = scale factor times delta amplitude (arcsec)
1814 * cos l 9.0019 = mean eccentricity
1815 * cos 2D 43.6
1816 * cos F 11.2 (latitude term)
1817 */
1818 }
1819
mean_elements_pl()1820 void mean_elements_pl()
1821 {
1822 /* Mean longitudes of planets (Laskar, Bretagnon) */
1823 Ve = mods3600( 210664136.4335482 * T + 655127.283046 );
1824 Ve += ((((((((
1825 -9.36e-023 * T
1826 - 1.95e-20 ) * T
1827 + 6.097e-18 ) * T
1828 + 4.43201e-15 ) * T
1829 + 2.509418e-13 ) * T
1830 - 3.0622898e-10 ) * T
1831 - 2.26602516e-9 ) * T
1832 - 1.4244812531e-5 ) * T
1833 + 0.005871373088 ) * T2;
1834 Ea = mods3600( 129597742.26669231 * T + 361679.214649 );
1835 Ea += (((((((( -1.16e-22 * T
1836 + 2.976e-19 ) * T
1837 + 2.8460e-17 ) * T
1838 - 1.08402e-14 ) * T
1839 - 1.226182e-12 ) * T
1840 + 1.7228268e-10 ) * T
1841 + 1.515912254e-7 ) * T
1842 + 8.863982531e-6 ) * T
1843 - 2.0199859001e-2 ) * T2;
1844 Ma = mods3600( 68905077.59284 * T + 1279559.78866 );
1845 Ma += (-1.043e-5*T + 9.38012e-3)*T2;
1846 Ju = mods3600( 10925660.428608 * T + 123665.342120 );
1847 Ju += (1.543273e-5*T - 3.06037836351e-1)*T2;
1848 Sa = mods3600( 4399609.65932 * T + 180278.89694 );
1849 Sa += (( 4.475946e-8*T - 6.874806E-5 ) * T + 7.56161437443E-1)*T2;
1850 }
1851
1852 /* Calculate geometric coordinates of true interpolated Moon apsides
1853 */
swi_intp_apsides(double J,double * pol,int ipli)1854 int swi_intp_apsides(double J, double *pol, int ipli)
1855 {
1856 double dd;
1857 double rsv[3];
1858 double sNF, sD, sLP, sMP, sM, sVe, sEa, sMa, sJu, sSa, fM, fVe, fEa, fMa, fJu, fSa, cMP, zMP, fNF, fD, fLP;
1859 double dMP, mLP, mNF, mD, mMP;
1860 int i, ii, iii, niter = 4; /* niter: silence compiler warning */
1861 ii=1;
1862 zMP=27.55454988;
1863 fNF = 27.212220817/zMP;/**/
1864 fD = 29.530588835/zMP;/**/
1865 fLP = 27.321582/zMP;/**/
1866 fM = 365.2596359/zMP;
1867 fVe = 224.7008001/zMP;
1868 fEa = 365.2563629/zMP;
1869 fMa = 686.9798519/zMP;
1870 fJu = 4332.589348/zMP;
1871 fSa = 10759.22722/zMP;
1872 T = (J-J2000)/36525.0;
1873 T2 = T*T;
1874 T4 = T2*T2;
1875 mean_elements();
1876 mean_elements_pl();
1877 sNF = NF;
1878 sD = D;
1879 sLP = SWELP;
1880 sMP = MP;
1881 sM = M ;
1882 sVe = Ve;
1883 sEa = Ea;
1884 sMa = Ma;
1885 sJu = Ju;
1886 sSa = Sa;
1887 sNF = mods3600(NF);
1888 sD = mods3600(D);
1889 sLP = mods3600(SWELP);
1890 sMP = mods3600(MP);
1891 if (ipli == SEI_INTP_PERG) {MP = 0.0; niter = 5;}
1892 if (ipli == SEI_INTP_APOG) {MP = 648000.0; niter = 4;}
1893 cMP = 0;
1894 dd = 18000.0;
1895 for (iii= 0; iii<=niter; iii++) {/**/
1896 dMP = sMP - MP;
1897 mLP = sLP - dMP;
1898 mNF = sNF - dMP;
1899 mD = sD - dMP;
1900 mMP = sMP - dMP;
1901 for (ii = 0; ii <=2; ii++) {/**/
1902 MP = mMP + (ii-1)*dd; /**/
1903 NF = mNF + (ii-1)*dd/fNF;
1904 D = mD + (ii-1)*dd/fD;
1905 SWELP = mLP + (ii-1)*dd/fLP;
1906 M = sM + (ii-1)*dd/fM ;
1907 Ve = sVe + (ii-1)*dd/fVe;
1908 Ea = sEa + (ii-1)*dd/fEa;
1909 Ma = sMa + (ii-1)*dd/fMa;
1910 Ju = sJu + (ii-1)*dd/fJu;
1911 Sa = sSa + (ii-1)*dd/fSa;
1912 moon1();
1913 moon2();
1914 moon3();
1915 moon4();
1916 if (ii==1) {
1917 for( i=0; i<3; i++ ) pol[i] = moonpol[i];
1918 }
1919 rsv[ii] = moonpol[2];
1920 }
1921 cMP = (1.5*rsv[0] - 2*rsv[1] + 0.5*rsv[2]) / (rsv[0] + rsv[2] - 2*rsv[1]);/**/
1922 cMP *= dd;
1923 cMP = cMP - dd;
1924 mMP += cMP;
1925 MP = mMP;
1926 dd /= 10;
1927 }
1928 return(0);
1929 }
1930
1931