1 // SPDX-License-Identifier: LGPL-2.1-or-later
2 //
3 // SPDX-FileCopyrightText: 2014 Gerhard Holtkamp
4 //
5
6 /* =========================================================================
7 Procedures needed for standard astronomy programs.
8 The majority of procedures in this unit are taken from Montenbruck,
9 Pfleger "Astronomie mit dem Personal Computer", Springer Verlag, 1989
10 as well as from the "Explanatory Supplement to the Astronomical Almanac"
11 University Science Books, Mill Valley, California, 1992
12 and modified correspondingly.
13
14 License: GNU LGPL Version 2+
15 Copyright : Gerhard HOLTKAMP 15-APR-2015
16 ========================================================================= */
17
18 #include <cmath>
19 using namespace std;
20
21 #include "astrolib.h"
22
frac(double f)23 double frac (double f)
24 { return fmod(f,1.0); }
25
26 /*--------------- Function ddd ------------------------------------------*/
27
ddd(int d,int m,double s)28 double ddd (int d, int m, double s)
29 {
30 /*
31 Conversion from Degrees, Minutes, Seconds into decimal degrees
32
33 INPUT :
34 d : degrees
35 m : minutes
36 s : seconds (and fractions)
37
38 OUTPUT :
39 dd : decimal degrees
40 */
41 double dd;
42
43 if ( (d < 0) || (m < 0) || (s < 0)) dd = -1.0;
44 else dd = 1.0;
45 dd = dd * (fabs(double(d)) + fabs(double(m))/60.0 + fabs(s)/3600.0);
46
47 return dd;
48 }
49
50 /*--------------- Function dms ------------------------------------------*/
51
dms(double dd,int & d,int & m,double & s)52 void dms (double dd, int &d, int &m, double &s)
53 /*
54 Conversion from decimal degrees into Degrees, Minutes, Seconds
55
56 INPUT :
57 dd : decimal degrees
58
59 OUTPUT :
60 d : degrees
61 m : minutes
62 s : seconds (and fractions)
63 */
64 {
65 double d1;
66
67 d1 = fabs(dd);
68 // d = int(floor(d1));
69 d = int(d1);
70 d1 = (d1 - d) * 60.0;
71 // m = int(floor(d1));
72 m = int(d1);
73 s = (d1 - m) * 60.0;
74 if (dd < 0)
75 {
76 if (d != 0) d = -d;
77 else
78 {
79 if (m != 0) m = -m;
80 else s = -s;
81 };
82 };
83 }
84
85 /*------------ FUNCTION mjd ----------------------------------------------*/
86
mjd(int day,int month,int year,double hour)87 double mjd (int day, int month, int year, double hour)
88 /*
89 Modified Julian Date ( MJD = Julian Date - 2400000.5)
90 valid for every date
91 Julian Calendar up to 4-OCT-1582,
92 Gregorian Calendar from 15-OCT-1582.
93 */
94 {
95 double a;
96 long int b, c;
97
98 a = 10000.0 * year + 100.0 * month + day;
99 if (month <= 2)
100 {
101 month = month + 12;
102 year = year - 1;
103 };
104 if (a <= 15821004.1)
105 {
106 b = ((year+4716)/4) - 1181;
107 if (year < -4716)
108 {
109 c = year + 4717;
110 c = -c;
111 c = c / 4;
112 c = -c;
113 b = c - 1182;
114 };
115 }
116 else b = (year/400) - (year/100) + (year/4);
117 // { b = -2 + floor((year+4716)/4) - 1179;}
118 // else {b = floor(year/400) - floor(year/100) + floor(year/4);};
119
120 a = 365.0 * year - 679004.0;
121 a = a + b + int(30.6001 * (month + 1)) + day + hour / 24.0;
122
123 return a;
124 }
125
126 /*---------------- Function julcent ---------------------------------------*/
127
julcent(double mjuld)128 double julcent (double mjuld)
129 {
130 /*
131 Julian Centuries since J2000.0 of Modified Julian Date mjuld.
132 */
133 return (mjuld - 51544.5) / 36525.0;
134 }
135
136 /*---------------- Function caldat -----------------------------------------*/
137
caldat(double mjd,int & day,int & month,int & year,double & hour)138 void caldat (double mjd, int &day, int &month, int &year, double &hour)
139 /*
140 Calculate the calendar date from the Modified Julian Date
141
142 INPUT :
143 mjd : Modified Julian Date (Julian Date - 2400000.5)
144
145 OUTPUT :
146 day, month, year : corresponding date
147 hour : corresponding hours of the above date
148 */
149 {
150 long int b, c, d, e, f, jd0;
151
152 jd0 = long(mjd + 2400001.0);
153 if (jd0 < 2299161) c = jd0 + 1524; /* Julian calendar */
154 else
155 { /* Gregorian calendar */
156 b = long (( jd0 - 1867216.25) / 36524.25);
157 c = jd0 + b - (b/4) + 1525;
158 };
159
160 if (mjd < -2400001.0) // special case for year < -4712
161 {
162 if (mjd == floor(mjd)) jd0 = jd0 + 1;
163 c = long((-jd0 - 0.1)/ 365.25);
164 c = c + 1;
165 year = -4712 - c;
166 d = c / 4;
167 d = c * 365 + d; // number of days from JAN 1, -4712
168 f = d + jd0; // day of the year
169 if ((c % 4) == 0) e = 61;
170 else e = 60;
171 if (f == 0)
172 {
173 year = year - 1;
174 month = 12;
175 day = 31;
176 f = 500; // set as a marker
177 };
178 if (f < e)
179 {
180 if (f < 32)
181 {
182 month = 1;
183 day = f;
184 }
185 else
186 {
187 month = 2;
188 day = f - 31;
189 };
190 }
191 else
192 {
193 if (f < 500)
194 {
195 f = f - e;
196 month = long((f + 123.0) / 30.6001);
197 day = f - long(month * 30.6001) + 123;
198 month = month - 1;
199 };
200 };
201 }
202 else // normal case
203 {
204 d = long ((c - 122.1) / 365.25);
205 e = 365 * d + (d/4);
206 f = long ((c - e) / 30.6001);
207 day = c - e - long(30.6001 * f);
208 month = f - 1 - 12 * (f / 14);
209 year = d - 4715 - ((7 + month) / 10);
210 };
211
212 hour = 24.0 * (mjd - floor(mjd));
213 }
214
215 /*---------------- Function DefTdUt --------------------------------------*/
DefTdUt(int yi)216 double DefTdUt (int yi)
217 {
218 /*
219 Get a suitable default for value of TDT - UT in year yi in seconds.
220 */
221
222 const int td[9] = {55,55,56,56,57,58,58,59,60};
223 double t, result;
224 int yy, yr;
225
226 yr = yi;
227
228 if (yr > 1899)
229 {
230 if (yr >= 2000) // this corrects for observations until 2013
231 {
232 if (yr > 2006) yr -= 1;
233 if (yr > 2006) yr -= 1; // no leap second at the end of 2007
234 if (yr > 2007) yr -= 1; // no leap second at the end of 2009
235 if (yr > 2007) yr -= 1; // no leap second at the end of 2010
236 if (yr > 2008) yr -= 1; // no leap second at the end of 2012
237 yr -= 6;
238 if (yr < 1999) yr = 1999;
239 };
240 t = yr - 2000;
241 t = t / 100.0;
242 result = (((-339.84*t - 516.12)*t -160.22)*t + 92.23)*t + 71.28;
243 if (yr > 2013)
244 {
245 t = yr - 2013;
246 t = t / 100.0;
247 result = (27.5*t + 75.0)*t + 73.0;
248 }
249 else if (yr > 1994)
250 {
251 result -= 6.28;
252 }
253 else if (yr > 1985)
254 {
255 yy = yr - 1986;
256 result = td[yy];
257 };
258 };
259
260 if (yr < 1900)
261 {
262 t = yr - 1800;
263 t = t / 100;
264 if (yr > 1649)
265 {
266 t = t - 0.19;
267 result = 5.156 + 13.3066 * t * t;
268 if (yr > 1864) result = -0.6*(yr - 1865) + 6;
269 if (yr > 1884) result = -0.2*(yr - 1885) - 6;
270 }
271 else
272 {
273 if (yr > 947) result = 25.5 * t*t;
274 else result = 1360 + (44.3*t + 320.0) * t;
275 };
276 };
277
278 // round to full seconds
279
280 if (result < 0)
281 {
282 t = -result;
283 t += 0.5;
284 result = - floor(t);
285 }
286 else result = floor(result + 0.5);
287
288 result += 0.184; // because of TT = TAI + 32.184;
289
290 return result;
291 }
292
293 /*---------------- Function lsidtim -------------------------------------*/
294
lsidtim(double jd,double lambda,double ep2)295 double lsidtim (double jd, double lambda, double ep2)
296 {
297 /* Calculate the Apparent Local Mean Sidereal Time (in decimal hours) for
298 the Modified Julian Date jd and geographic longitude lambda (in degrees)
299 and with the correction (AST - MST) eps (given in seconds)
300 */
301 double lmst;
302 double t, ut, gmst;
303 int mjd0;
304
305 mjd0 = int(jd);
306 ut = (jd - mjd0)*24.0;
307 t = (mjd0 - 51544.5) / 36525.0;
308 gmst = 6.697374558 + 1.0027379093*ut
309 + (8640184.812866 + (0.093104 - 6.2e-6*t)*t)*t/3600.0;
310 lmst = 24.0 * frac((gmst + lambda/15.0) / 24.0);
311 lmst = lmst + ep2 / 3600.0;
312
313 return lmst;
314 }
315
316 /*---------------- Function eps -----------------------------------------*/
317
eps(double t)318 double eps (double t) // obliquity of ecliptic
319 {
320 /*
321 Obliquety of ecliptic eps (in radians)
322 at time t (in Julian Centuries from J2000.0)
323 */
324 double tp;
325
326 tp = 23.43929111 - (46.815+(0.00059-0.001813*t)*t)*t/3600.0;
327 tp = 1.74532925199e-2 * tp;
328
329 return tp;
330 }
331
332 /*---------------- Function eclequ -----------------------------------------*/
333
eclequ(double t,Vec3 & r1)334 Vec3 eclequ (double t, Vec3& r1) // ecliptic -> equatorial
335 {
336 /*
337 Convert position vector r1 from ecliptic into equatorial coordinates
338 at t (in Julian Centuries from J2000.0 )
339 */
340
341 Mat3 m;
342 Vec3 r2;
343
344 m = xrot (-eps(t));
345 r2 = mxvct (m, r1);
346 return r2;
347 }
348
349 /*---------------- Function equecl -----------------------------------------*/
350
equecl(double t,Vec3 & r1)351 Vec3 equecl (double t, Vec3& r1) // equatorial -> ecliptic
352 {
353 /*
354 Convert position vector r1 from equatorial into ecliptic coordinates
355 at t (in Julian Centuries from J2000.0)
356 */
357
358 Mat3 m;
359 Vec3 r2;
360
361 m = xrot (eps(t));
362 r2 = mxvct (m, r1);
363 return r2;
364 }
365
366 /*---------------- Function pmatecl -----------------------------------------*/
367
pmatecl(double t1,double t2)368 Mat3 pmatecl (double t1, double t2) // ecl. precession
369 {
370 /*
371 Calculate ecliptic precession matrix a from t1 to t2
372 (times in Julian Centuries from J2000.0)
373 */
374
375 const double secrad = 4.8481368111e-6; // arcsec -> radians
376
377 double ppi, pii, pa, dt;
378 Mat3 m1, m2, a;
379
380 dt = t2 - t1;
381 ppi = 174.876383889 + (((3289.4789+0.60622*t1)*t1) +
382 ((-869.8089-0.50491*t1) + 0.03536*dt)*dt) / 3600.0;
383 ppi = ppi * 1.74532925199e-2;
384 pii = ((47.0029-(0.06603-0.000598*t1)*t1) +
385 ((-0.03302+0.000598*t1)+0.00006*dt)*dt)*dt * secrad;
386 pa = ((5029.0966+(2.22226-0.000042*t1)*t1) +
387 ((1.11113-0.000042*t1)-0.000006*dt)*dt)*dt * secrad;
388
389 pa = ppi + pa;
390 m1 = zrot (-pa);
391 m2 = xrot (pii);
392 m1 = m1 * m2;
393 m2 = zrot (ppi);
394 a = m1 * m2;
395
396 return a;
397 }
398
399 /*---------------- Function pmatequ --------------------------------------*/
400
pmatequ(double t1,double t2)401 Mat3 pmatequ (double t1, double t2) // equ. precession
402 {
403 /*
404 Calculate equatorial precession matrix a from t1 to t2
405 (times in Julian Centuries from J2000.0)
406 */
407 const double secrad = 4.8481368111e-6; // arcsec -> radians
408
409 double dt, zeta, z, theta;
410 Mat3 m1, m2, a;
411
412 dt = t2 - t1;
413 zeta = ((2306.2181+(1.39656-0.000139*t1)*t1) +
414 ((0.30188-0.000345*t1)+0.017998*dt)*dt)*dt * secrad;
415 z = zeta + ((0.7928+0.000411*t1)+0.000205*dt)*dt*dt * secrad;
416 theta = ((2004.3109-(0.8533+0.000217*t1)*t1) -
417 ((0.42665+0.000217*t1)+0.041833*dt)*dt)*dt * secrad;
418
419 m1 = zrot (-z);
420 m2 = yrot (theta);
421 m1 = m1 * m2;
422 m2 = zrot (-zeta);
423 a = m1 * m2;
424
425 return a;
426 }
427
428 /*---------------- Function nutmat --------------------------------------*/
429
nutmat(double t,double & ep2,bool hpr)430 Mat3 nutmat (double t, double& ep2, bool hpr)
431 {
432 /*
433 Calculate nutation matrix a from mean to true equatorial coordinates
434 at time t (in Julian Centuries from J2000.0)
435 Also calculates the correction ep2 for apparent sidereal time in sec
436
437 if hpr is true a high precision is used, otherwise a low precision
438 (only the first 50 terms of the nutation theory are used)
439 */
440 const int ntb1 = 15;
441 const int tb1[ntb1][5] =
442 {
443 { 0, 0, 0, 0, 1}, // 1
444 { 0, 0, 0, 0, 2}, // 2
445 { 0, 0, 2,-2, 2}, // 9
446 { 0, 1, 0, 0, 0}, // 10
447 { 0, 1, 2,-2, 2}, // 11
448 { 0,-1, 2,-2, 2}, // 12
449 { 0, 0, 2,-2, 1}, // 13
450 { 0, 2, 0, 0, 0}, // 16
451 { 0, 2, 2,-2, 2}, // 18
452 { 0, 0, 2, 0, 2}, // 31
453 { 1, 0, 0, 0, 0}, // 32
454 { 0, 0, 2, 0, 1}, // 33
455 { 1, 0, 2, 0, 2}, // 34
456 { 1, 0, 0, 0, 1}, // 38
457 { -1, 0, 0, 0, 1}, // 39
458 };
459
460 const int ntb2 = 35;
461 const int tb2[ntb2][5] =
462 {
463 { -2, 0, 2, 0, 1}, // 3
464 { 2, 0,-2, 0, 0}, // 4
465 { -2, 0, 2, 0, 2}, // 5
466 { 1,-1, 0,-1, 0}, // 6
467 { 0,-2, 2,-2, 1}, // 7
468 { 2, 0,-2, 0, 1}, // 8
469 { 2, 0, 0,-2, 0}, // 14
470 { 0, 0, 2,-2, 0}, // 15
471 { 0, 1, 0, 0, 1}, // 17
472 { 0,-1, 0, 0, 1}, // 19
473 { -2, 0, 0, 2, 1}, // 20
474 { 0,-1, 2,-2, 1}, // 21
475 { 2, 0, 0,-2, 1}, // 22
476 { 0, 1, 2,-2, 1}, // 23
477 { 1, 0, 0,-1, 0}, // 24
478 { 2, 1, 0,-2, 0}, // 25
479 { 0, 0,-2, 2, 1}, // 26
480 { 0, 1,-2, 2, 0}, // 27
481 { 0, 1, 0, 0, 2}, // 28
482 { -1, 0, 0, 1, 1}, // 29
483 { 0, 1, 2,-2, 0}, // 30
484 { 1, 0, 0,-2, 0}, // 35
485 { -1, 0, 2, 0, 2}, // 36
486 { 0, 0, 0, 2, 0}, // 37
487 { -1, 0, 2, 2, 2}, // 40
488 { 1, 0, 2, 0, 1}, // 41
489 { 0, 0, 2, 2, 2}, // 42
490 { 2, 0, 0, 0, 0}, // 43
491 { 1, 0, 2,-2, 2}, // 44
492 { 2, 0, 2, 0, 2}, // 45
493 { 0, 0, 2, 0, 0}, // 46
494 { -1, 0, 2, 0, 1}, // 47
495 { -1, 0, 0, 2, 1}, // 48
496 { 1, 0, 0,-2, 1}, // 49
497 { -1, 0, 2, 2, 1}, // 50
498 };
499
500 const double tb3[ntb1][4] =
501 {
502 {-171996.0,-174.2, 92025.0, 8.9 }, // 1
503 { 2062.0, 0.2, -895.0, 0.5 }, // 2
504 { -13187.0, -1.6, 5736.0, -3.1 }, // 9
505 { 1426.0, -3.4, 54.0, -0.1 }, // 10
506 { -517.0, 1.2, 224.0, -0.6 }, // 11
507 { 217.0, -0.5, -95.0, 0.3 }, // 12
508 { 129.0, 0.1, -70.0, 0.0 }, // 13
509 { 17.0, -0.1, 0.0, 0.0 }, // 16
510 { -16.0, 0.1, 7.0, 0.0 }, // 18
511 { -2274.0, -0.2, 977.0, -0.5 }, // 31
512 { 712.0, 0.1, -7.0, 0.0 }, // 32
513 { -386.0, -0.4, 200.0, 0.0 }, // 33
514 { -301.0, 0.0, 129.0, -0.1 }, // 34
515 { 63.0, 0.1, -33.0, 0.0 }, // 38
516 { -58.0, -0.1, 32.0, 0.0 }, // 39
517 };
518
519 const double tb4[ntb2][2] =
520 {
521 { 46.0, -24.0}, // 3
522 { 11.0, 0.0 }, // 4
523 { -3.0, 1.0 }, // 5
524 { -3.0, 0.0 }, // 6
525 { -2.0, 1.0 }, // 7
526 { 1.0, 0.0 }, // 8
527 { 48.0, 1.0 }, // 14
528 {-22.0, 0.0 }, // 15
529 {-15.0, 9.0 }, // 17
530 {-12.0, 6.0 }, // 19
531 { -6.0, 3.0 }, // 20
532 { -5.0, 3.0 }, // 21
533 { 4.0, -2.0 }, // 22
534 { 4.0, -2.0 }, // 23
535 { -4.0, 0.0 }, // 24
536 { 1.0, 0.0 }, // 25
537 { 1.0, 0.0 }, // 26
538 { -1.0, 0.0 }, // 27
539 { 1.0, 0.0 }, // 28
540 { 1.0, 0.0 }, // 29
541 { -1.0, 0.0 }, // 30
542 {-158.0,-1.0 }, // 35
543 {123.0,-53.0 }, // 36
544 { 63.0, -2.0 }, // 37
545 {-59.0, 26.0 }, // 40
546 {-51.0, 27.0 }, // 41
547 {-38.0, 16.0 }, // 42
548 { 29.0, -1.0 }, // 43
549 { 29.0,-12.0 }, // 44
550 {-31.0, 13.0 }, // 45
551 { 26.0, -1.0 }, // 46
552 { 21.0,-10.0 }, // 47
553 { 16.0, -8.0 }, // 48
554 {-13.0, 7.0 }, // 49
555 {-10.0, 5.0 }, // 50
556 };
557
558 const double secrad = 4.8481368111e-6; // arcsec -> radians
559 const double p2 = 2.0 * M_PI;
560
561 double ls, lm, d, f, n, dpsi, deps, ep0;
562 int j;
563 Mat3 m1, m2, a;
564
565 if (hpr)
566 {
567 lm = 2.355548393544 +
568 (8328.691422883903 + (0.000151795164 + 0.000000310281*t)*t)*t;
569 ls = 6.240035939326 +
570 (628.301956024185 + (-0.000002797375 - 0.000000058178*t)*t)*t;
571 f = 1.627901933972 +
572 (8433.466158318464 + (-0.000064271750 + 0.000000053330*t)*t)*t;
573 d = 5.198469513580 +
574 (7771.377146170650 + (-0.000033408511 + 0.000000092115*t)*t)*t;
575 n = 2.182438624361 +
576 (-33.757045933754 + (0.000036142860 + 0.000000038785*t)*t)*t;
577
578 lm = fmod(lm,p2);
579 ls = fmod(ls,p2);
580 f = fmod(f,p2);
581 d = fmod(d,p2);
582 n = fmod(n,p2);
583
584 dpsi = 0.0;
585 deps = 0.0;
586 for(j=0; j<ntb1; j++)
587 {
588 ep0 = tb1[j][0]*lm + tb1[j][1]*ls + tb1[j][2]*f + tb1[j][3]*d + tb1[j][4]*n;
589 dpsi = dpsi + (tb3[j][0]+tb3[j][1]*t) * sin(ep0);
590 deps = deps + (tb3[j][2]+tb3[j][3]*t) * cos(ep0);
591 };
592 for(j=0; j<ntb2; j++)
593 {
594 ep0 = tb2[j][0]*lm + tb2[j][1]*ls + tb2[j][2]*f + tb2[j][3]*d + tb2[j][4]*n;
595 dpsi = dpsi + tb4[j][0] * sin(ep0);
596 deps = deps + tb4[j][1] * cos(ep0);
597 };
598 dpsi = 1.0e-4 * dpsi * secrad;
599 deps = 1.0e-4 * deps * secrad;
600 }
601
602 else // low precision
603 {
604 ls = p2 * frac (0.993133+99.997306*t); // mean anomaly sun
605 d = p2 * frac (0.827362+1236.853087*t); // diff long. moon-sun
606 f = p2 * frac (0.259089+1342.227826*t); // dist. node
607 n = p2 * frac (0.347346 - 5.372447*t); // long. node
608
609 dpsi = (-17.2*sin(n) - 1.319*sin(2*(f-d+n)) - 0.227*sin(2*(f+n))
610 + 0.206*sin(2*n) + 0.143*sin(ls)) * secrad;
611 deps = (+9.203*cos(n) + 0.574*cos(2*(f-d+n)) + 0.098*cos(2*(f+n))
612 -0.09*cos(2*n) ) * secrad;
613 };
614
615 ep0 = eps (t);
616 ep2 = ep0 + deps;
617 m1 = xrot (ep0);
618 m2 = zrot (-dpsi);
619 m1 *= m2;
620 m2 = xrot (-ep2);
621 a = m2 * m1;
622 ep2 = dpsi * cos (ep2);
623 ep2 *= 13750.9870831; // convert radians into time-seconds
624
625 return a;
626 }
627
628 /*---------------- Function nutecl --------------------------------------*/
629
nutecl(double t,double & ep2)630 Mat3 nutecl (double t, double& ep2) // nutation matrix (ecliptic)
631 {
632 /*
633 Calculate nutation matrix a from mean to true ecliptic coordinates
634 at time t (in Julian Centuries from J2000.0)
635 Also calculates the correction ep2 for apparent sidereal time in sec
636 */
637
638 const double secrad = 4.8481368111e-6; // arcsec -> radians
639 const double p2 = 2.0 * M_PI;
640
641 double ls, d, f, n, dpsi, deps, ep0;
642 Mat3 m1, m2, a;
643
644 ls = p2 * frac (0.993133+99.997306*t); // mean anomaly sun
645 d = p2 * frac (0.827362+1236.853087*t); // diff long. moon-sun
646 f = p2 * frac (0.259089+1342.227826*t); // dist. node
647 n = p2 * frac (0.347346 - 5.372447*t); // long. node
648
649 dpsi = (-17.2*sin(n) - 1.319*sin(2*(f-d+n)) - 0.227*sin(2*(f+n))
650 + 0.206*sin(2*n) + 0.143*sin(ls)) * secrad;
651
652 deps = (+9.203*cos(n) + 0.574*cos(2*(f-d+n)) + 0.098*cos(2*(f+n))
653 -0.09*cos(2*n) ) * secrad;
654
655 ep0 = eps (t);
656 ep2 = ep0 + deps;
657 m1 = xrot (-deps);
658 m2 = zrot (-dpsi);
659 a = m1 * m2;
660 ep2 = dpsi * cos (ep2);
661 ep2 *= 13750.9870831; // convert radians into time-seconds
662
663 return a;
664 }
665
666 /*---------------- Function pmatequ --------------------------------------*/
667
PoleMx(double xp,double yp)668 Mat3 PoleMx (double xp, double yp)
669 {
670 /* Returns Polar Motion matrix.
671 xp and yp are the coordinates of the Celestial Ephemeris Pole
672 with respect to the IERS Reference Pole in arcsec as published
673 in the IERS Bulletin B.
674 */
675 const double arctrd = M_PI / (180.0*3600.0);
676 Mat3 res;
677 double xr, yr;
678
679 xr = xp * arctrd;
680 yr = yp * arctrd;
681 res.assign(1.0,0.0,xr,0.0,1.0,-yr,-xr,yr,1.0);
682
683 return res;
684 }
685
686 /*---------------- Function aberrat --------------------------------------*/
687
aberrat(double t,Vec3 & ve)688 Vec3 aberrat (double t, Vec3& ve) // aberration
689 {
690 /*
691 Correct position vector ve for aberration into new position va
692 at time t (in Julian Centuries from J2000.0)
693 */
694 const double p2 = 2.0 * M_PI;
695
696 double l, cl, d0;
697 Vec3 va;
698
699 d0 = abs(ve);
700 l = p2 * frac(0.27908+100.00214*t);
701 cl = cos(l)*d0;
702 va[0] = ve[0] - 9.934e-5 * sin(l) * d0;
703 va[1] = ve[1] + 9.125e-5 * cl;
704 va[2] = ve[2] + 3.927e-5 * cl;
705
706 return va;
707 }
708
709 /*--------------------- Function GeoPos ----------------------------------*/
710
GeoPos(double jd,double ep2,double lat,double lng,double ht)711 Vec3 GeoPos (double jd, double ep2, double lat, double lng, double ht)
712 {
713 /* Return the geocentric vector (in the equatorial system) of the
714 geographic position given by the latitude lat and longitude lng
715 (in radians) and height ht (in m) at the MJD-time jd (UT).
716 ep2 : correction for apparent sidereal time in sec
717
718 The length unit of the vector is in terms of the equatorial
719 Earth radius (6378.137 km)
720 */
721 double const e2 = 6.69438499959e-3;
722 double np, h, sp;
723
724 Vec3 r;
725
726 sp = sin(lat);
727 h = ht / 6378.137e3;
728 np = 1.0 / (sqrt(1.0 - e2*sp*sp));
729 r[2] = ((1.0 - e2)*np + h)*sp;
730
731 sp = (np + h) * cos(lat);
732 np = lsidtim(jd, (lng*180.0/M_PI), ep2) * M_PI/12.0;
733
734 r[0] = sp * cos(np);
735 r[1] = sp * sin(np);
736
737 return r;
738 }
739
GeoPos(double jd,double ep2,double lat,double lng,double ht,double xp,double yp)740 Vec3 GeoPos (double jd, double ep2, double lat, double lng, double ht,
741 double xp, double yp)
742 {
743 /* Return the geocentric vector (in the equatorial system) of the
744 geographic position given by the latitude lat and longitude lng
745 (in radians) and height ht (in m) at the MJD-time jd (UT).
746
747 ep2 : correction for apparent sidereal time in sec
748 xp, yp: coordinates of polar motion in arcsec.
749
750 The length unit of the vector is in terms of the equatorial
751 Earth radius (6378.137 km)
752 */
753 double const e2 = 6.69438499959e-3;
754 double np, h, sp;
755 Mat3 prx;
756 Vec3 r;
757
758 sp = sin(lat);
759 h = ht / 6378.137e3;
760 np = 1.0 / (sqrt(1.0 - e2*sp*sp));
761 r[2] = ((1.0 - e2)*np + h)*sp;
762 sp = (np + h) * cos(lat);
763 r[0] = sp * cos(lng);
764 r[1] = sp * sin(lng);
765
766 // correct for polar motion
767 if ((xp != 0) || (yp != 0))
768 {
769 prx = mxtrn(PoleMx(xp, yp));
770 r = mxvct(prx, r);
771 };
772
773 // convert into equatorial
774 np = lsidtim(jd, 0.0, ep2) * M_PI/12.0;
775 prx = zrot(-np);
776 r = mxvct(prx, r);
777
778 return r;
779 }
780
781 /*--------------------- Function EquHor ----------------------------------*/
782
EquHor(double jd,double ep2,double lat,double lng,Vec3 r)783 Vec3 EquHor (double jd, double ep2, double lat, double lng, Vec3 r)
784 {
785 /* convert vector r from the equatorial system into the horizontal system.
786 jd = MJD-time (UT)
787 ep2 : correction for apparent sidereal time in sec
788 lat, lng : geographic latitude and longitude (in radians)
789 */
790 double lst;
791 Vec3 s;
792 Mat3 mx;
793
794 lst = lsidtim(jd, (lng*180.0/M_PI), ep2) * M_PI/12.0;
795 mx = zrot(lst);
796 s = mxvct(mx, r);
797 mx = yrot(M_PI/2.0 - lat);
798 s = mxvct(mx, s);
799
800 return s;
801 }
802
803 /*--------------------- Function HorEqu ----------------------------------*/
804
805
HorEqu(double jd,double ep2,double lat,double lng,Vec3 r)806 Vec3 HorEqu (double jd, double ep2, double lat, double lng, Vec3 r)
807 {
808 /* convert vector r from the horizontal system into the equatorial system.
809 jd = MJD-time (UT)
810 ep2 : correction for apparent sidereal time in sec
811 lat, lng : geographic latitude and longitude (in radians)
812 */
813 double lst;
814 Vec3 s;
815 Mat3 mx;
816
817 mx = yrot(lat - M_PI/2.0);
818 s = mxvct(mx, r);
819 lst = -lsidtim(jd, (lng*180.0/M_PI), ep2) * M_PI/12.0;
820 mx = zrot(lst);
821 s = mxvct(mx, s);
822
823 return s;
824 }
825
826 /*--------------------- Function AppPos ----------------------------------*/
827
AppPos(double jd,double ep2,double lat,double lng,double ht,int solsys,const Vec3 & r,double & azim,double & elev,double & dist)828 void AppPos (double jd, double ep2, double lat, double lng, double ht,
829 int solsys, const Vec3& r, double& azim, double& elev, double& dist)
830 {
831 /* get apparent position in the horizontal system
832 jd = MJD-time (UT)
833 ep2 : correction for apparent sidereal time in sec
834 lat, lng : geographic latitude and longitude (in radians)
835 ht : height above normal in meters.
836 solsys : = 1 if object is in solar system and parallax has to be
837 taken into account, 0 otherwise.
838 r = vector of celestial object. The unit of length of this vector
839 has to be in terms of the equatorial Earth radius (6378.14 km)
840 if solsys = 1, otherwise it's arbitrary.
841 azim : azimuth in radians (0 is to the North).
842 elev : elevation in radians
843 dist : distance (if solsys = 1; otherwise abs(r))
844 */
845 Vec3 s;
846
847 // correct for topocentric position (parallax)
848 if (solsys) s = r - GeoPos(jd, ep2, lat, lng, ht);
849 else s = r;
850
851 s = EquHor(jd, ep2, lat, lng, s);
852 s = carpol(s);
853 dist = s[0];
854 elev = s[2];
855 azim = M_PI - s[1];
856 }
857
858 /*--------------------- Function AppRADec ----------------------------------*/
859
AppRADec(double jd,double ep2,double lat,double lng,double azim,double elev,double & ra,double & dec)860 void AppRADec (double jd, double ep2, double lat, double lng,
861 double azim, double elev, double& ra, double& dec)
862 {
863 /* get apparent position in the horizontal system
864 jd = MJD-time (UT)
865 ep2 : correction for apparent sidereal time in sec
866 lat, lng : geographic latitude and longitude (in radians)
867 azim : azimuth in radians (0 is to the North).
868 elev : elevation in radians
869 ra : Right Ascension (in radians)
870 dec : Declination (in radians)
871 */
872 Vec3 s;
873
874 s[0] = 1.0;
875 s[1] = M_PI - azim;
876 s[2] = elev;
877 s = polcar(s);
878 s = HorEqu(jd, ep2, lat, lng, s);
879 s = carpol(s);
880 dec = s[2];
881 ra = s[1];
882 }
883
884 /*------------------------- Refract ---------------------------------*/
885
Refract(double h,double p,double t)886 double Refract (double h, double p, double t)
887 {
888 /* Calculate atmospheric refraction with low precision
889 h = height (in radians) of object
890 p = presure (in millibars)
891 t = temperature (in degrees Celsius)
892
893 RETURN: refraction angle in radians
894 */
895 double const raddeg = 180.0 / M_PI;
896 double r;
897
898 r = h * raddeg;
899 r = (r + 7.31 / (r + 4.4)) / raddeg;
900 r = (0.28*p/(t+273.0)) * 0.0167 / tan(r);
901 r = r / raddeg;
902
903 return r;
904 }
905
906 /*---------------- Function eccanom --------------------------------------*/
907
eccanom(double man,double ecc)908 double eccanom (double man, double ecc)
909 {
910 /*
911 Solve Kepler equation for eccentric anomaly of elliptical orbits.
912 man : Mean Anomaly (in radians)
913 ecc : eccentricity
914 */
915 const double p2 = 2.0*M_PI;
916 const double eps = 1E-11;
917 const int maxit = 15;
918
919 double m, e, f;
920 int i, mi;
921
922 m=man/p2;
923 mi = int(m);
924 m=p2*(m-mi);
925 if (m < 0) m=m+p2;
926 if (ecc < 0.8) e=m;
927 else e=M_PI;
928 f = e - ecc*sin(e) - m;
929 i=0;
930
931 while ((fabs(f) > eps) && (i < maxit))
932 {
933 e = e - f / (1.0 - ecc*cos(e));
934 f = e - ecc*sin(e) - m;
935 i = i + 1;
936 }
937
938 return e;
939 }
940
941
942 /*---------------- Function hypanom --------------------------------------*/
943
hypanom(double mh,double ecc)944 double hypanom (double mh, double ecc)
945 {
946 /*
947 Solve Kepler equation for eccentric anomaly of hyperbolic orbits.
948 mh : Mean Anomaly
949 ecc : eccentricity
950 */
951 const double eps = 1E-11;
952 const int maxit = 15;
953
954 double h, f;
955 int i;
956
957 h = log (2.0*fabs(mh)/ecc+1.8);
958 if (mh < 0.0) h = -h;
959 f = ecc * sinh(h) - h - mh;
960 i = 0;
961
962 while ((fabs(f) > eps*(1.0+fabs(h+mh))) && (i < maxit))
963 {
964 h = h - f / (ecc*cosh(h) - 1.0);
965 f = ecc*sinh(h) - h - mh;
966 i = i + 1;
967 }
968
969 return h;
970 }
971
972
973 /*------------------ Function ellip ------------------------------------*/
974
ellip(double gm,double t0,double t,double a,double ecc,double m0,Vec3 & r1,Vec3 & v1)975 void ellip (double gm, double t0, double t, double a, double ecc,
976 double m0, Vec3& r1, Vec3& v1)
977 {
978 /*
979 Calculate position r1 and velocity v1 of for elliptic orbit at time t.
980 gm : Gravitational Constant
981 t0 : epoch
982 a : semim-ajor axis
983 ecc : eccentricity
984 m0 : Mean Anomaly at epoch (in radians).
985 The units must be consistent with each other.
986 */
987 double m, e, fac, k, c, s;
988
989 if (fabs(a) < 1e-60) a = 1e-60;
990 k = gm / a;
991 if (k >= 0) k = sqrt(k);
992 else k = 0; // just in case
993
994 m = k * (t - t0) / a + m0;
995 e = eccanom(m, ecc);
996 fac = sqrt(1.0 - ecc*ecc);
997 c = cos(e);
998 s = sin(e);
999 r1.assign (a*(c - ecc), a*fac*s, 0.0);
1000 m = 1.0 - ecc*c;
1001 v1.assign (-k*s/m, k*fac*c/m, 0.0);
1002 }
1003
1004 /*---------------- Function hyperb --------------------------------------*/
1005
hyperb(double gm,double t0,double t,double a,double ecc,Vec3 & r1,Vec3 & v1)1006 void hyperb (double gm, double t0, double t, double a, double ecc,
1007 Vec3& r1, Vec3& v1)
1008 {
1009 /*
1010 Calculate position r1 and velocity v1 of for hyperbolic orbit at time t.
1011 gm : Gravitational Constant
1012 t0 : time of perihelion passage
1013 a : semi-major axis
1014 ecc : eccentricity
1015 The units must be consistent with each other.
1016 */
1017 double mh, h, fac, k, c, s;
1018
1019 a = fabs(a);
1020 if (a < 1e-60) a = 1e-60;
1021 k = gm / a;
1022 if (k >= 0) k = sqrt(k);
1023 else k = 0; // just in case
1024
1025 mh = k * (t - t0) / a;
1026 h = hypanom (mh, ecc);
1027 fac = sqrt(ecc*ecc-1.0);
1028 c = cosh(h);
1029 s = sinh(h);
1030 r1.assign (a*(ecc-c), a*fac*s, 0.0);
1031 mh = ecc*c - 1.0;
1032 v1.assign (-k*s/mh, k*fac*c/mh, 0.0);
1033 }
1034
1035 /*---------------- Function parab --------------------------------------*/
1036
stumpff(double e2,double & c1,double & c2,double & c3)1037 void stumpff (double e2, double& c1, double& c2, double& c3)
1038 {
1039 /*
1040 Calculation of Stumpff functions c1=sin(e)/c,
1041 c2=(1-cos(e))/e^2 and c3=(e-sin(e))/e^3 for the
1042 argument e2=e^2 (e: eccentric anomaly)
1043 */
1044 const double eps=1E-12;
1045 double n, add;
1046
1047 c1=0.0; c2=0.0; c3=0.0; add=1.0; n=1.0;
1048 do
1049 {
1050 c1=c1+add; add=add/(2.0*n);
1051 c2=c2+add; add=add/(2.0*n+1);
1052 c3=c3+add; add=-e2*add;
1053 n=n+1.0;
1054 } while (abs(add)>eps);
1055 }
1056
parab(double gm,double t0,double t,double q,double ecc,Vec3 & r1,Vec3 & v1)1057 void parab (double gm, double t0, double t, double q, double ecc,
1058 Vec3& r1, Vec3& v1)
1059 {
1060 /*
1061 Calculate position r1 and velocity v1 of for parabolic
1062 (or near parabolic) orbit at time t using Stumpff's method.
1063 gm : Gravitational Constant
1064 t0 : time of perihelion passage
1065 q : perihelion distance
1066 ecc : eccentricity
1067 The units must be consistent with each other.
1068 */
1069 const double eps = 1E-9;
1070 const int maxit = 15;
1071
1072 double e2, e20, fac, c1, c2, c3, k, tau, a, u, u2;
1073 double x, y;
1074 int i, fle;
1075
1076 ecc = fabs(ecc);
1077 e2=0.0; fac=0.5*ecc; i=0;
1078
1079 q = fabs(q);
1080 if (q < 1e-40) q = 1e-40;
1081 k = gm / (q*(1+ecc));
1082
1083 if (k >= 0) k = sqrt(k);
1084 else k = 0; // just in case
1085
1086 tau = gm / (q*q*q);
1087 if (tau >= 0) tau= 1.5*sqrt(tau)*(t-t0);
1088 else tau = 0;
1089
1090 do
1091 {
1092 i = i + 1;
1093 e20 = e2;
1094 if (fac < 0.0) a = -1.0;
1095 else a = tau * sqrt(fac);
1096 a = sqrt(a*a+1.0)+a;
1097 if (a > 0.0) a = exp(log(a)/3.0);
1098 if (a == 0.0) u = 0.0;
1099 else u = a - 1.0/a;
1100 u2 = u*u;
1101 if (fac != 0) e2 = u2*(1.0-ecc)/fac;
1102 else e2 = 1;
1103 stumpff (e2, c1, c2, c3);
1104 fac = 3.0*ecc*c3;
1105 fle = (abs(e2-e20) < eps) || (i > maxit);
1106 } while (!fle);
1107
1108 if (fac != 0)
1109 {
1110 tau = q*(1.0+u2*c2*ecc/fac);
1111 x = q*(1.0-u2*c2/fac);
1112 y = (1.0+ecc)/fac;
1113 if (y >= 0) y = q*sqrt(y)*u*c1;
1114 else y = 0;
1115 r1.assign (x, y, 0.0);
1116 v1.assign (-k*y/tau, k*(x/tau+ecc), 0.0);
1117 }
1118 else // just in case
1119 {
1120 r1.assign (0, 0, 0);
1121 v1.assign (0, 0, 0);
1122 };
1123 }
1124
1125 /*---------------- Function kepler --------------------------------------*/
1126
kepler(double gm,double t0,double t,double m0,double a,double ecc,double ran,double aper,double inc,Vec3 & r1,Vec3 & v1)1127 void kepler (double gm, double t0, double t, double m0, double a, double ecc,
1128 double ran, double aper, double inc, Vec3& r1, Vec3& v1)
1129 {
1130 /*
1131 Calculate position r1 and velocity v1 of body at time t as an
1132 undisturbed 2-body Kepler problem.
1133
1134 gm : Gravitational Constant
1135 t0 : time of perihelion passage (hyperbolic or near-parabolic)
1136 epoch of m0 for elliptical orbits.
1137 m0 : Mean Anomaly at epoch for elliptical orbits in degrees in the
1138 range 0 to 360. If m0 < 0 the calculations will be done as
1139 near-parabolic. Set m0 = 0 for hyperbolic orbits.
1140 a : semi-major axis for elliptical (positive) or hyperbolic orbits
1141 (negative). If m0 < 0, a signifies perihelion distance q instead
1142 ecc : eccentricity
1143 ran : Right Ascension of Ascending Node (in degrees)
1144 aper : Argument of Perihelion (in degrees)
1145 inc : Inclination (in degrees)
1146
1147 The units must be consistent with each other.
1148
1149 NOTE : If ecc = 1 the orbit will always be calculated as a parabolic one.
1150 If ecc < 1 (elliptical orbits) two possibilities exist:
1151 If m0 >= 0 the calculation will be done in the classical
1152 way with m0 signifying the mean anomaly at time t0 and
1153 a the semi-major axis.
1154 If m0 < 0 the orbit will be calculated as a near-parabolic
1155 orbit with a signifying the perihelion distance at time t0.
1156 (This latter method will be useful for high eccentricity
1157 comet orbits for which typically q is given instead of a)
1158 If ecc > 1 (hyperbolic orbits) two possibilities exist:
1159 If m0 >= 0 (the value is ignored for hyperbolic orbits)
1160 a classical hyperbolic calculation is performed with a
1161 signifying the semi-major axis (negative for hyperbolas).
1162 If m0 < 0 the orbit will be calculated as a near-parabolic
1163 orbit with a signifying the perihelion distance. (This can
1164 be useful for comet orbits which would rarely have an
1165 eccentricity >> 1)
1166 In either case t0 signifies the time of perihelion passage.
1167 */
1168 const double dgrd = M_PI / 180.0;
1169 enum {ell, par, hyp} kepc;
1170 Mat3 m1,m2;
1171
1172 kepc = ell; // just to keep the compiler happy
1173
1174 // convert into radians
1175 m0 *= dgrd;
1176 ran *= dgrd;
1177 aper *= dgrd;
1178 inc *= dgrd;
1179
1180 // calculate the position and velocity within the orbit plane
1181 if (ecc == 1.0) kepc = par;
1182 if (ecc < 1.0)
1183 {
1184 if (m0 < 0.0) kepc = par;
1185 else kepc = ell;
1186 }
1187 if (ecc > 1.0)
1188 {
1189 if (m0 < 0.0) kepc = par;
1190 else kepc = hyp;
1191 }
1192
1193 switch (kepc)
1194 {
1195 case ell : ellip (gm, t0, t, a, ecc, m0, r1, v1); break;
1196 case par : parab (gm, t0, t, a, ecc, r1, v1); break;
1197 case hyp : hyperb (gm, t0, t, a, ecc, r1, v1); break;
1198 }
1199
1200 // convert into reference plane
1201 m1 = zrot (-aper);
1202 m2 = xrot (-inc);
1203 m1 *= m2;
1204 m2 = zrot (-ran);
1205 m2 = m2 * m1;
1206
1207 r1 = mxvct (m2, r1);
1208 v1 = mxvct (m2, v1);
1209 }
1210
1211 /*---------------- Function oscelm --------------------------------------*/
1212
oscelm(double gm,double t,Vec3 & r1,Vec3 & v1,double & t0,double & m0,double & a,double & ecc,double & ran,double & aper,double & inc)1213 void oscelm (double gm, double t, Vec3& r1, Vec3& v1,
1214 double& t0, double& m0, double& a, double& ecc,
1215 double& ran, double& aper, double& inc)
1216 {
1217 /*
1218 Get the osculating Kepler elements at epoch t from position r1 and
1219 velocity v1 of body.
1220
1221 gm : Gravitational Constant. If set negative, its absolute value will
1222 be taken as gm and the perihelion distance will be returned
1223 instead of the semimajor axis.
1224 t0 : time of perihelion passage
1225 m0 : Mean Anomaly at epoch for elliptical orbits in degrees in the
1226 range 0 to 360. Set m0 = 0 for hyperbolic orbits.
1227 m0 will be set to -1 if perihelion distance is to be returned
1228 instead of a.
1229 a : semi-major axis for elliptical (positive) or hyperbolic orbits
1230 (negative). If m0 < 0, a signifies perihelion distance q instead.
1231 If gm is negative, the perihelion distance q will be returned
1232 ecc : eccentricity
1233 ran : Right Ascension of Ascending Node (in degrees)
1234 aper : Argument of Perihelion (in degrees)
1235 inc : Inclination (in degrees)
1236
1237 The units must be consistent with each other.
1238 */
1239 const double rddg = 180.0/M_PI;
1240 const double p2 = 2.0*M_PI;
1241
1242 Vec3 c;
1243 double cabs, u, cu, su, p, r;
1244 int pflg; // 1, if perihelion distance to be returned
1245
1246 pflg = 0;
1247 if (gm < 0.0)
1248 {
1249 gm = -gm;
1250 pflg = 1;
1251 };
1252
1253 if (gm < 1e-60) gm = 1e-60;
1254 c = r1 * v1;
1255 cabs = abs(c);
1256 if (fabs(cabs) < 1e-40) cabs = 1e-40;
1257 ran = atan20 (c[0], -c[1]);
1258 inc = c[2]/cabs;
1259 if (fabs(inc) <= 1.0) inc = acos (inc);
1260 else inc = 0;
1261
1262 r = abs(r1);
1263 if (fabs(r) < 1e-40) r = 1e-40;
1264 su = sin(inc);
1265 if (su != 0.0) su = r1[2]/su;
1266 cu = r1[0]*cos(ran)+r1[1]*sin(ran);
1267 u = atan20 (su, cu); // argument of latitude
1268
1269 a = abs(v1);
1270 a = 2.0/r - a*a/gm;
1271 if (fabs(a) < 1.0E-30) ecc = 1.0; // parabolic
1272 else
1273 {
1274 a = 1.0/a; // semimajor axis
1275 ecc = 0;
1276 };
1277
1278 p = cabs*cabs/gm; // semilatum rectum
1279
1280 if (ecc == 1.0)
1281 {
1282 p = p / 2.0; // perihelion distance
1283 a = 2*p;
1284 }
1285 else
1286 {
1287 ecc = 1.0 - p/a;
1288 if (ecc >= 0) ecc = sqrt(ecc);
1289 else ecc = 0;
1290 p = p / (1.0 + ecc); // perihelion distance
1291 }
1292
1293 if (fabs(a) > 1e-60) cu = (1.0 - r/a);
1294 else cu = 0; // just in case
1295 su = dot(r1,v1) / sqrt(fabs(a)*gm);
1296 if (ecc < 1.0)
1297 {
1298 m0 = atan20(su,cu); // eccentric anomaly
1299 su = sin(m0);
1300 cu = cos(m0);
1301 aper = 1.0-ecc*ecc;
1302 if (aper >= 0) aper = atan20(sqrt(aper)*su,(cu-ecc)); // true anomaly
1303 m0 = m0 - ecc*su; // mean anomaly
1304 }
1305 else if (ecc > 1.0)
1306 {
1307 su = su / ecc;
1308 m0 = su + sqrt(su*su + 1.0);
1309 if (m0 >= 0) m0 = log(m0); // = asinh (su); hyperbolic anomaly
1310 aper = (ecc+1.0)/(ecc-1.0);
1311 if (aper >= 0) aper = 2.0 * atan(sqrt(aper)*tanh(m0/2.0));
1312 m0 = ecc*su-m0;
1313 }
1314
1315 if (ecc != 1.0)
1316 {
1317 aper = u - aper; // argument of perihelion
1318 u = fabs(a);
1319 t0 = u/gm;
1320 if (t0 >= 0) t0 = t - u*sqrt(t0)*m0; // time of perihelion transit
1321 else t0 = t;
1322 }
1323 else // parabolic
1324 {
1325 pflg = 1;
1326
1327 aper = 2.0 * atan(su); // true anomaly
1328 aper = u - aper; // argument of perihelion
1329 t0 = 2.0*p*p*p/gm;
1330 if (t0 >= 0) t0 = t - sqrt(t0) * (su*su*su/3.0 + su);
1331 else t0 = t;
1332 }
1333
1334 if (m0 < 0.0) m0 = m0 + p2;
1335 if (ran < 0.0) ran = ran + p2;
1336 if (aper < 0.0) aper = aper + p2;
1337 m0 = rddg * m0;
1338 ran = rddg * ran;
1339 aper = rddg * aper;
1340 inc = rddg * inc;
1341
1342 if (ecc > 1.0) m0 = 0.0;
1343 if (pflg)
1344 {
1345 a = p;
1346 m0 = -1.0;
1347 }
1348 }
1349
1350 /*-------------------- QuickSun --------------------------------------*/
1351
QuickSun(double t)1352 Vec3 QuickSun (double t) // low precision position of the Sun at time t
1353 {
1354 /* Low precision position of the Sun at time t given in
1355 Julian centuries since J2000.
1356 Valid only between 1950 and 2050.
1357
1358 Returns the position vector (in A.U.) in ecliptic coordinates
1359
1360 */
1361 double n, g, l;
1362 Vec3 rs;
1363
1364 n = 36525.0 * t;
1365 l = 280.460 + 0.9856474*n; // mean longitude
1366 g = 357.528 + 0.9856003*n; // mean anomaly
1367 l = 2.0*M_PI * frac(l/360.0);
1368 g = 2.0*M_PI * frac(g/360.0);
1369 l = l + (1.915*sin(g) + 0.020*sin(2.0*g)) * 1.74532925199e-2; // ecl.long.
1370 g = 1.00014 - 0.01671*cos(g) - 0.00014*cos(2.0*g); // radius
1371 rs[0] = g * cos(l);
1372 rs[1] = g * sin(l);
1373 rs[2] = 0;
1374
1375 return rs;
1376 }
1377
1378 /*-------------------- Class Sun200 --------------------------------------*/
1379 /*
1380 Ecliptic coordinates (in A.U.) and velocity (in A.U./day)
1381 of the sun for Equinox of Date given in Julian centuries since J2000.
1382 ======================================================================
1383 */
Sun200()1384 Sun200::Sun200 ()
1385 { }
1386
position(double t)1387 Vec3 Sun200::position (double t) // position of the Sun at time t
1388 {
1389 Vec3 rs, vs;
1390
1391 state (t, rs, vs);
1392
1393 return rs;
1394 }
1395
state(double t,Vec3 & rs,Vec3 & vs)1396 void Sun200::state (double t, Vec3& rs, Vec3& vs)
1397 {
1398 /* State vector rs (position) and vs (velocity) of the Sun in
1399 ecliptic of date coordinates at time t (in Julian Centuries
1400 since J2000).
1401 */
1402 const double p2 = 2.0 * M_PI;
1403
1404 double l, b, r;
1405 int i;
1406
1407 tt = t;
1408
1409 dl = 0.0; dr = 0.0; db = 0.0;
1410 m2 = p2 * frac(0.1387306 + 162.5485917*t);
1411 m3 = p2 * frac(0.9931266 + 99.9973604*t);
1412 m4 = p2 * frac(0.0543250 + 53.1666028*t);
1413 m5 = p2 * frac(0.0551750 + 8.4293972*t);
1414 m6 = p2 * frac(0.8816500 + 3.3938722*t);
1415 d = p2 * frac(0.8274 + 1236.8531*t);
1416 a= p2 * frac(0.3749 + 1325.5524*t);
1417 uu = p2 * frac(0.2591 + 1342.2278*t);
1418
1419 c3[1] = 1.0; s3[1] = 0.0;
1420 c3[2] = cos(m3); s3[2] = sin(m3);
1421 c3[0] = c3[2]; s3[0] = -s3[2];
1422
1423 for (i=3; i<9; i++)
1424 addthe(c3[i-1],s3[i-1],c3[2],s3[2],c3[i],s3[i]);
1425 pertven(); pertmar(); pertjup(); pertsat(); pertmoo();
1426
1427 dl = dl + 6.4 * sin(p2*(0.6983 + 0.0561*t))
1428 + 1.87 * sin(p2*(0.5764 + 0.4174*t))
1429 + 0.27 * sin(p2*(0.4189 + 0.3306*t))
1430 + 0.20 * sin(p2*(0.3581 + 2.4814*t));
1431 l = p2 * frac(0.7859453 + m3/p2 +((6191.2+1.1*t)*t+dl)/1296.0E3);
1432 r = 1.0001398 - 0.0000007*t + dr*1E-6;
1433 b = db * 4.8481368111E-6;
1434
1435 cl= cos(l); sl=sin(l); cb=cos(b); sb=sin(b);
1436 rs[0]=r*cl*cb; rs[1]=r*sl*cb; rs[2]=r*sb;
1437
1438 /* velocity calculation
1439 gms=2.9591202986e-4 AU^3/d^2
1440 e=0.0167086
1441 sqrt(gms/a)=1.7202085e-2 AU/d
1442 sqrt(gms/a)*sqrt(1-e^2)=1.71996836e-2 AU/d
1443 nu=M + 2e*sin(M) = M + 0.0334172 * sin(M) */
1444
1445 uu = m3 + 0.0334172*sin(m3); // eccentric Anomaly E
1446 d = cos(uu);
1447 uu = sin(uu);
1448 a = 1.0 - 0.0167086*d;
1449 vs[0] = -1.7202085e-2*uu/a; // velocity in orbit plane
1450 vs[1] = 1.71996836e-2*d/a;
1451 uu = atan2 (0.9998604*uu, (d-0.0167086)); // true anomaly
1452 d = cos(uu);
1453 uu = sin(uu);
1454 dr = d*vs[0]+uu*vs[1];
1455 dl = (d*vs[1]-uu*vs[0]) / r;
1456
1457 vs[0] = dr*cl*cb - dl*r*sl*cb;
1458 vs[1] = dr*sl*cb + dl*r*cl*cb;
1459 vs[2] = dr*sb;
1460 }
1461
addthe(double c1,double s1,double c2,double s2,double & cc,double & ss)1462 void Sun200::addthe (double c1, double s1, double c2, double s2,
1463 double& cc, double& ss)
1464 {
1465 cc=c1*c2-s1*s2;
1466 ss=s1*c2+c1*s2;
1467 }
1468
term(int i1,int i,int it,double dlc,double dls,double drc,double drs,double dbc,double dbs)1469 void Sun200::term (int i1, int i, int it, double dlc, double dls, double drc,
1470 double drs, double dbc, double dbs)
1471 {
1472 if (it == 0) addthe (c3[i1+1],s3[i1+1],c[i+8],s[i+8],u,v);
1473 else
1474 {
1475 u=u*tt;
1476 v=v*tt;
1477 }
1478 dl = dl + dlc*u + dls*v;
1479 dr = dr + drc*u + drs*v;
1480 db = db + dbc*u + dbs*v;
1481 }
1482
pertven()1483 void Sun200::pertven() // Kepler terms and perturbations by Venus
1484 {
1485 int i;
1486
1487 c[8]=1.0; s[8]=0.0; c[7]=cos(m2); s[7]=-sin(m2);
1488 for (i=7; i>2; i--)
1489 addthe(c[i],s[i],c[7],s[7],c[i-1],s[i-1]);
1490
1491 term (1, 0, 0, -0.22, 6892.76, -16707.37, -0.54, 0.0, 0.0);
1492 term (1, 0, 1, -0.06, -17.35, 42.04, -0.15, 0.0, 0.0);
1493 term (1, 0, 2, -0.01, -0.05, 0.13, -0.02, 0.0, 0.0);
1494 term (2, 0, 0, 0.0, 71.98, -139.57, 0.0, 0.0, 0.0);
1495 term (2, 0, 1, 0.0, -0.36, 0.7, 0.0, 0.0, 0.0);
1496 term (3, 0, 0, 0.0, 1.04, -1.75, 0.0, 0.0, 0.0);
1497 term (0, -1, 0, 0.03, -0.07, -0.16, -0.07, 0.02, -0.02);
1498 term (1, -1, 0, 2.35, -4.23, -4.75, -2.64, 0.0, 0.0);
1499 term (1, -2, 0, -0.1, 0.06, 0.12, 0.2, 0.02, 0.0);
1500 term (2, -1, 0, -0.06, -0.03, 0.2, -0.01, 0.01, -0.09);
1501 term (2, -2, 0, -4.7, 2.9, 8.28, 13.42, 0.01, -0.01);
1502 term (3, -2, 0, 1.8, -1.74, -1.44, -1.57, 0.04, -0.06);
1503 term (3, -3, 0, -0.67, 0.03, 0.11, 2.43, 0.01, 0.0);
1504 term (4, -2, 0, 0.03, -0.03, 0.1, 0.09, 0.01, -0.01);
1505 term (4, -3, 0, 1.51, -0.4, -0.88, -3.36, 0.18, -0.1);
1506 term (4, -4, 0, -0.19, -0.09, -0.38, 0.77, 0.0, 0.0);
1507 term (5, -3, 0, 0.76, -0.68, 0.3, 0.37, 0.01, 0.0);
1508 term (5, -4, 0, -0.14, -0.04, -0.11, 0.43, -0.03, 0.0);
1509 term (5, -5, 0, -0.05, -0.07, -0.31, 0.21, 0.0, 0.0);
1510 term (6, -4, 0, 0.15, -0.04, -0.06, -0.21, 0.01, 0.0);
1511 term (6, -5, 0, -0.03, -0.03, -0.09, 0.09, -0.01, 0.0);
1512 term (6, -6, 0, 0.0, -0.04, -0.18, 0.02, 0.0, 0.0);
1513 term (7, -5, 0, -0.12, -0.03, -0.08, 0.31, -0.02, -0.01);
1514 }
1515
pertmar()1516 void Sun200::pertmar() // Kepler terms and perturbations by Mars
1517 {
1518 int i;
1519
1520 c[7] = cos(m4); s[7] = -sin(m4);
1521 for (i=7; i>0; i--)
1522 addthe(c[i],s[i],c[7],s[7],c[i-1],s[i-1]);
1523 term (1, -1, 0, -0.22, 0.17, -0.21, -0.27, 0.0, 0.0);
1524 term (1, -2, 0, -1.66, 0.62, 0.16, 0.28, 0.0, 0.0);
1525 term (2, -2, 0, 1.96, 0.57, -1.32, 4.55, 0.0, 0.01);
1526 term (2, -3, 0, 0.4, 0.15, -0.17, 0.46, 0.0, 0.0);
1527 term (2, -4, 0, 0.53, 0.26, 0.09, -0.22, 0.0, 0.0);
1528 term (3, -3, 0, 0.05, 0.12, -0.35, 0.15, 0.0, 0.0);
1529 term (3, -4, 0, -0.13, -0.48, 1.06, -0.29, 0.01, 0.0);
1530 term (3, -5, 0, -0.04, -0.2, 0.2, -0.04, 0.0, 0.0);
1531 term (4, -4, 0, 0.0, -0.03, 0.1, 0.04, 0.0, 0.0);
1532 term (4, -5, 0, 0.05, -0.07, 0.2, 0.14, 0.0, 00);
1533 term (4, -6, 0, -0.1, 0.11, -0.23, -0.22, 0.0, 0.0);
1534 term (5, -7, 0, -0.05, 0.0, 0.01, -0.14, 0.0, 0.0);
1535 term (5, -8, 0, 0.05, 0.01, -0.02, 0.1, 0.0, 0.0);
1536 }
1537
pertjup()1538 void Sun200::pertjup() // Kepler terms and perturbations by Jupiter
1539 {
1540 int i;
1541
1542 c[7] = cos(m5); s[7] = -sin(m5);
1543 for (i=7; i>4; i--)
1544 addthe(c[i],s[i],c[7],s[7],c[i-1],s[i-1]);
1545 term (1, -1, 0, 0.01, 0.07, 0.18, -0.02, 0.0, -0.02);
1546 term (0, -1, 0, -0.31, 2.58, 0.52, 0.34, 0.02, 0.0);
1547 term (1, -1, 0, -7.21, -0.06, 0.13, -16.27, 0.0, -0.02);
1548 term (1, -2, 0, -0.54, -1.52, 3.09, -1.12, 0.01, -0.17);
1549 term (1, -3, 0, -0.03, -0.21, 0.38, -0.06, 0.0, -0.02);
1550 term (2, -1, 0, -0.16, 0.05, -0.18, -0.31, 0.01, 0.0);
1551 term (2, -2, 0, 0.14, -2.73, 9.23, 0.48, 0.0, 0.0);
1552 term (2, -3, 0, 0.07, -0.55, 1.83, 0.25, 0.01, 0.0);
1553 term (2, -4, 0, 0.02, -0.08, 0.25, 0.06, 0.0, 0.0);
1554 term (3, -2, 0, 0.01, -0.07, 0.16, 0.04, 0.0, 0.0);
1555 term (3, -3, 0, -0.16, -0.03, 0.08, -0.64, 0.0, 0.0);
1556 term (3, -4, 0, -0.04, -0.01, 0.03, -0.17, 0.0, 0.0);
1557 }
1558
pertsat()1559 void Sun200::pertsat() // Kepler terms and perturbations by Saturn
1560 {
1561 c[7] = cos(m6); s[7] = -sin(m6);
1562 addthe(c[7],s[7],c[7],s[7],c[6],s[6]);
1563 term (0, -1, 0, 0.0, 0.32, 0.01, 0.0, 0.0, 0.0);
1564 term (1, -1, 0, -0.08, -0.41, 0.97, -0.18, 0.0, -0.01);
1565 term (1, -2, 0, 0.04, 0.1, -0.23, 0.1, 0.0, 0.0);
1566 term (2, -2, 0, 0.04, 0.1, -0.35, 0.13, 0.0, 0.0);
1567 }
1568
pertmoo()1569 void Sun200::pertmoo() // corrections for Earth-Moon center of gravity
1570 {
1571 dl = dl + 6.45*sin(d) - 0.42*sin(d-a) + 0.18*sin(d+a) + 0.17*sin(d-m3) - 0.06*sin(d+m3);
1572 dr = dr + 30.76*cos(d) - 3.06*cos(d-a) + 0.85*cos(d+a) - 0.58*cos(d+m3) + 0.57*cos(d-m3);
1573 db = db + 0.576*sin(uu);
1574 }
1575
1576 /*--------------------- Class moon200 ----------------------------------*/
1577
1578 /*
1579 Position vector (in Earth radii) of the Moon referred to Ecliptic
1580 for Equinox of Date.
1581 t is the time in Julian centuries since J2000.
1582 */
1583
Moon200()1584 Moon200::Moon200 ()
1585 { }
1586
position(double t)1587 Vec3 Moon200::position (double t) // position of the Moon at time t
1588 {
1589 Vec3 rm;
1590
1591 const double arc = 206264.81; // 3600*180/pi = ''/r
1592 const double pi2 = M_PI * 2.0;
1593
1594 double lambda, beta, r, fac;
1595
1596 minit(t);
1597 solar1(); solar2(); solar3(); solarn(n); planetary(t);
1598
1599 lambda = pi2 * frac ((l0+dlam/arc) / pi2);
1600 if (lambda < 0) lambda = lambda + pi2;
1601 s = f + ds / arc;
1602
1603 fac = 1.000002708 + 139.978 * dgam;
1604 beta = (fac*(18518.511+1.189+gam1c)*sin(s)-6.24*sin(3*s)+n);
1605 beta = beta * 4.8481368111e-6;
1606 sinpi = sinpi * 0.999953253;
1607 r = arc / sinpi;
1608
1609 rm.assign (r, lambda, beta);
1610 rm = polcar(rm);
1611
1612 return rm;
1613 }
1614
1615
addthe(double c1,double s1,double c2,double s2,double & c,double & s)1616 void Moon200::addthe (double c1, double s1, double c2, double s2,
1617 double& c, double& s)
1618 {
1619 c=c1*c2-s1*s2;
1620 s=s1*c2+c1*s2;
1621 }
1622
sinus(double phi)1623 double Moon200::sinus (double phi)
1624 {
1625 /* sin(phi) in units of 1r = 360ø */
1626 return sin (2.0*M_PI * frac(phi));
1627 }
1628
long_periodic(double t)1629 void Moon200::long_periodic (double t)
1630 {
1631 /* long periodic changes of mean elements */
1632
1633 double s1, s2, s3, s4, s5, s6, s7;
1634
1635 s1 = sinus (0.19833 + 0.05611*t);
1636 s2 = sinus (0.27869 + 0.04508*t);
1637 s3 = sinus (0.16827 - 0.36903*t);
1638 s4 = sinus (0.34734 - 5.37261*t);
1639 s5 = sinus (0.10498 - 5.37899*t);
1640 s6 = sinus (0.42681 - 0.41855*t);
1641 s7 = sinus (0.14943 - 5.37511*t);
1642 dl0 = 0.84*s1 + 0.31*s2 + 14.27*s3 + 7.26*s4 + 0.28*s5 + 0.24*s6;
1643 dl = 2.94*s1 + 0.31*s2 + 14.27*s3 + 9.34*s4 + 1.12*s5 + 0.83*s6;
1644 dls = -6.4*s1 - 1.89*s6;
1645 df = 0.21*s1+0.31*s2+14.27*s3-88.7*s4-15.3*s5+0.24*s6-1.86*s7;
1646 dd = dl0 - dls;
1647 dgam = -3332E-9 * sinus (0.59734 - 5.37261*t)
1648 -539E-9 * sinus (0.35498 - 5.37899*t)
1649 -64E-9 * sinus (0.39943 - 5.37511*t);
1650 }
1651
minit(double t)1652 void Moon200::minit(double t)
1653 {
1654 /* calculate mean elements l (moon), F (node distance)
1655 l' (sun) and D (moon's elongation) */
1656
1657 const double arc = 206264.81; // 3600*180/pi = ''/r
1658 const double pi2 = M_PI * 2.0;
1659
1660 int i, j, max;
1661 double t2, arg, fac;
1662
1663 max = 0; // just to keep the compiler happy
1664 arg = 0;
1665 fac = 0;
1666
1667 t2 = t * t;
1668 dlam = 0; ds = 0; gam1c = 0; sinpi = 3422.7;
1669 long_periodic (t);
1670 l0 = pi2*frac(0.60643382+1336.85522467*t-0.00000313*t2) + dl0/arc;
1671 l = pi2*frac(0.37489701+1325.55240982*t+0.00002565*t2) + dl/arc;
1672 ls = pi2*frac(0.99312619+99.99735956*t-0.00000044*t2) + dls/arc;
1673 f = pi2*frac(0.25909118+1342.22782980*t-0.00000892*t2) + df/arc;
1674 d = pi2*frac(0.82736186+1236.85308708*t-0.00000397*t2) + dd/arc;
1675 for (i=0; i<4; i++)
1676 {
1677 switch (i)
1678 {
1679 case 0: arg=l; max=4; fac=1.000002208; break;
1680 case 1: arg=ls; max=3; fac=0.997504612-0.002495388*t; break;
1681 case 2: arg=f; max=4; fac=1.000002708+139.978*dgam; break;
1682 case 3: arg=d; max=6; fac=1.0; break;
1683 }
1684 co[6][i] = 1.0; co[7][i] = cos (arg) * fac;
1685 si[6][i] = 0.0; si[7][i] = sin (arg) * fac;
1686 for (j = 2; j <= max; j++)
1687 addthe(co[j+5][i],si[j+5][i],co[7][i],si[7][i],co[j+6][i],si[j+6][i]);
1688 for (j = 1; j <= max; j++)
1689 {
1690 co[6-j][i]=co[j+6][i];
1691 si[6-j][i]=-si[j+6][i];
1692 }
1693 }
1694 }
1695
term(int p,int q,int r,int s,double & x,double & y)1696 void Moon200::term (int p, int q, int r, int s, double& x, double& y)
1697 {
1698 // calculate x=cos(p*arg1+q*arg2...); y=sin(p*arg1+q*arg2...)
1699 int i[4];
1700 int k;
1701
1702 i[0]=p; i[1]=q; i[2]=r; i[3]=s; x=1.0; y=0.0;
1703 for (k=0; k<4; k++)
1704 if (i[k]!=0) addthe(x,y,co[i[k]+6][k],si[i[k]+6][k],x,y);
1705 }
1706
addsol(double coeffl,double coeffs,double coeffg,double coeffp,int p,int q,int r,int s)1707 void Moon200::addsol(double coeffl, double coeffs, double coeffg,
1708 double coeffp, int p, int q, int r, int s)
1709 {
1710 double x, y;
1711 term (p,q,r,s,x,y);
1712 dlam = dlam + coeffl*y; ds = ds + coeffs*y;
1713 gam1c = gam1c + coeffg*x; sinpi = sinpi + coeffp*x;
1714 }
1715
solar1()1716 void Moon200::solar1()
1717 {
1718 addsol ( 13.902, 14.06, -0.001, 0.2607, 0, 0, 0, 4);
1719 addsol ( 0.403, -4.01, 0.394, 0.0023, 0, 0, 0, 3);
1720 addsol ( 2369.912, 2373.36, 0.601, 28.2333, 0, 0, 0, 2);
1721 addsol ( -125.154, -112.79, -0.725, -0.9781, 0, 0, 0, 1);
1722 addsol ( 1.979, 6.98, -0.445, 0.0433, 1, 0, 0, 4);
1723 addsol (191.953, 192.72, 0.029, 3.0861, 1, 0, 0, 2);
1724 addsol (-8.466, -13.51, 0.455, -0.1093, 1, 0, 0, 1);
1725 addsol (22639.5, 22609.07, 0.079, 186.5398, 1, 0, 0, 0);
1726 addsol (18.609, 3.59, -0.094, 0.0118, 1, 0, 0, -1);
1727 addsol (-4586.465, -4578.13, -0.077, 34.3117, 1, 0, 0, -2);
1728 addsol (3.215, 5.44, 0.192, -0.0386, 1, 0, 0, -3);
1729 addsol (-38.428, -38.64, 0.001, 0.6008, 1, 0, 0, -4);
1730 addsol (-0.393, -1.43, -0.092, 0.0086, 1, 0, 0, -6);
1731 addsol (-0.289, -1.59, 0.123, -0.0053, 0, 1, 0, 4);
1732 addsol (-24.420, -25.1, 0.04, -0.3, 0, 1, 0, 2);
1733 addsol (18.023, 17.93, 0.007, 0.1494, 0, 1, 0, 1);
1734 addsol (-668.146, -126.98, -1.302, -0.3997, 0, 1, 0, 0);
1735 addsol (0.56, 0.32, -0.001, -0.0037, 0, 1, 0, -1);
1736 addsol (-165.145, -165.06, 0.054, 1.9178, 0, 1, 0, -2);
1737 addsol (-1.877, -6.46, -0.416, 0.0339, 0, 1, 0, -4);
1738 addsol (0.213, 1.02, -0.074, 0.0054, 2, 0, 0, 4);
1739 addsol (14.387, 14.78, -0.017, 0.2833, 2, 0, 0, 2);
1740 addsol (-0.586, -1.2, 0.054, -0.01, 2, 0, 0, 1);
1741 addsol (769.016, 767.96, 0.107, 10.1657, 2, 0, 0, 0);
1742 addsol (1.75, 2.01, -0.018, 0.0155, 2, 0, 0, -1);
1743 addsol (-211.656, -152.53, 5.679, -0.3039, 2, 0, 0, -2);
1744 addsol (1.225, 0.91, -0.03, -0.0088, 2, 0, 0, -3);
1745 addsol (-30.773, -34.07, -0.308, 0.3722, 2, 0, 0, -4);
1746 addsol (-0.57, -1.4, -0.074, 0.0109, 2, 0, 0, -6);
1747 addsol (-2.921, -11.75, 0.787, -0.0484, 1, 1, 0, 2);
1748 addsol (1.267, 1.52, -0.022, 0.0164, 1, 1, 0, 1);
1749 addsol (-109.673, -115.18, 0.461, -0.949, 1, 1, 0, 0);
1750 addsol (-205.962, -182.36, 2.056, 1.4437, 1, 1, 0, -2);
1751 addsol (0.233, 0.36, 0.012, -0.0025, 1, 1, 0, -3);
1752 addsol (-4.391, -9.66, -0.471, 0.0673, 1, 1, 0, -4);
1753 }
1754
solar2()1755 void Moon200::solar2()
1756 {
1757 addsol (0.283, 1.53, -0.111, 0.006, 1, -1, 0, 4);
1758 addsol (14.577, 31.7, -1.54, 0.2302, 1, -1, 0, 2);
1759 addsol (147.687, 138.76, 0.679, 1.1528, 1, -1, 0, 0);
1760 addsol (-1.089, 0.55, 0.021, 0.0, 1, -1, 0, -1);
1761 addsol (28.475, 23.59, -0.443, -0.2257, 1, -1, 0, -2);
1762 addsol (-0.276, -0.38, -0.006, -0.0036, 1, -1, 0, -3);
1763 addsol (0.636, 2.27, 0.146, -0.0102, 1, -1, 0, -4);
1764 addsol (-0.189, -1.68, 0.131, -0.0028, 0, 2, 0, 2);
1765 addsol (-7.486, -0.66, -0.037, -0.0086, 0, 2, 0, 0);
1766 addsol (-8.096, -16.35, -0.74, 0.0918, 0, 2, 0, -2);
1767 addsol (-5.741, -0.04, 0.0, -0.0009, 0, 0, 2, 2);
1768 addsol (0.255, 0.0, 0.0, 0.0, 0, 0, 2, 1);
1769 addsol (-411.608, -0.2, 0.0, -0.0124, 0, 0, 2, 0);
1770 addsol (0.584, 0.84, 0.0, 0.0071, 0, 0, 2, -1);
1771 addsol (-55.173, -52.14, 0.0, -0.1052, 0, 0, 2, -2);
1772 addsol (0.254, 0.25, 0.0, -0.0017, 0, 0, 2, -3);
1773 addsol (0.025, -1.67, 0.0, 0.0031, 0, 0, 2, -4);
1774 addsol (1.06, 2.96, -0.166, 0.0243, 3, 0, 0, 2);
1775 addsol (36.124, 50.64, -1.3, 0.6215, 3, 0, 0, 0);
1776 addsol (-13.193, -16.4, 0.258, -0.1187, 3, 0, 0, -2);
1777 addsol (-1.187, -0.74, 0.042, 0.0074, 3, 0, 0, -4);
1778 addsol (-0.293, -0.31, -0.002, 0.0046, 3, 0, 0, -6);
1779 addsol (-0.29, -1.45, 0.116, -0.0051, 2, 1, 0, 2);
1780 addsol (-7.649, -10.56, 0.259, -0.1038, 2, 1, 0, 0);
1781 addsol (-8.627, -7.59, 0.078, -0.0192, 2, 1, 0, -2);
1782 addsol (-2.74, -2.54, 0.022, 0.0324, 2, 1, 0, -4);
1783 addsol (1.181, 3.32, -0.212, 0.0213, 2, -1, 0, 2);
1784 addsol (9.703, 11.67, -0.151, 0.1268, 2, -1, 0, 0);
1785 addsol (-0.352, -0.37, 0.001, -0.0028, 2, -1, 0, -1);
1786 addsol (-2.494, -1.17, -0.003, -0.0017, 2, -1, 0, -2);
1787 addsol (0.36, 0.2, -0.012, -0.0043, 2, -1, 0, -4);
1788 addsol (-1.167, -1.25, 0.008, -0.0106, 1, 2, 0, 0);
1789 addsol (-7.412, -6.12, 0.117, 0.0484, 1, 2, 0, -2);
1790 addsol (-0.311, -0.65, -0.032, 0.0044, 1, 2, 0, -4);
1791 addsol (0.757, 1.82, -0.105, 0.0112, 1, -2, 0, 2);
1792 addsol (2.58, 2.32, 0.027, 0.0196, 1, -2, 0, 0);
1793 addsol (2.533, 2.4, -0.014, -0.0212, 1, -2, 0, -2);
1794 addsol (-0.344, -0.57, -0.025, 0.0036, 0, 3, 0, -2);
1795 addsol (-0.992, -0.02, 0.0, 0.0, 1, 0, 2, 2);
1796 addsol (-45.099, -0.02, 0.0, -0.0010, 1, 0, 2, 0);
1797 addsol (-0.179, -9.52, 0.0, -0.0833, 1, 0, 2, -2);
1798 addsol (-0.301, -0.33, 0.0, 0.0014, 1, 0, 2, -4);
1799 addsol (-6.382, -3.37, 0.0, -0.0481, 1, 0, -2, 2);
1800 addsol (39.528, 85.13, 0.0, -0.7136, 1, 0, -2, 0);
1801 addsol (9.366, 0.71, 0.0, -0.0112, 1, 0, -2, -2);
1802 addsol (0.202, 0.02, 0.0, 0.0, 1, 0, -2, -4);
1803 }
1804
solar3()1805 void Moon200::solar3()
1806 {
1807 addsol (0.415, 0.1, 0.0, 0.0013, 0, 1, 2, 0);
1808 addsol (-2.152, -2.26, 0.0, -0.0066, 0, 1, 2, -2);
1809 addsol (-1.44, -1.3, 0.0, 0.0014, 0, 1, -2, 2);
1810 addsol (0.384, -0.04, 0.0, 0.0, 0, 1, -2, -2);
1811 addsol (1.938, 3.6, -0.145, 0.0401, 4, 0, 0, 0);
1812 addsol (-0.952, -1.58, 0.052, -0.0130, 4, 0, 0, -2);
1813 addsol (-0.551, -0.94, 0.032, -0.0097, 3, 1, 0, 0);
1814 addsol (-0.482, -0.57, 0.005, -0.0045, 3, 1, 0, -2);
1815 addsol (0.681, 0.96, -0.026, 0.0115, 3, -1, 0, 0);
1816 addsol (-0.297, -0.27, 0.002, -0.0009, 2, 2, 0, -2);
1817 addsol (0.254, 0.21, -0.003, 0.0, 2, -2, 0, -2);
1818 addsol (-0.25, -0.22, 0.004, 0.0014, 1, 3, 0, -2);
1819 addsol (-3.996, 0.0, 0.0, 0.0004, 2, 0, 2, 0);
1820 addsol (0.557, -0.75, 0.0, -0.009, 2, 0, 2, -2);
1821 addsol (-0.459, -0.38, 0.0, -0.0053, 2, 0, -2, 2);
1822 addsol (-1.298, 0.74, 0.0, 0.0004, 2, 0, -2, 0);
1823 addsol (0.538, 1.14, 0.0, -0.0141, 2, 0, -2, -2);
1824 addsol (0.263, 0.02, 0.0, 0.0, 1, 1, 2, 0);
1825 addsol (0.426, 0.07, 0.0, -0.0006, 1, 1, -2, -2);
1826 addsol (-0.304, 0.03, 0.0, 0.0003, 1, -1, 2, 0);
1827 addsol (-0.372, -0.19, 0.0, -0.0027, 1, -1, -2, 2);
1828 addsol (0.418, 0.0, 0.0, 0.0, 0, 0, 4, 0);
1829 addsol (-0.330, -0.04, 0.0, 0.0, 3, 0, 2, 0);
1830 }
1831
addn(double coeffn,int p,int q,int r,int s,double & n,double & x,double & y)1832 void Moon200::addn (double coeffn, int p, int q, int r, int s,
1833 double& n, double&x, double& y)
1834 {
1835 term (p,q,r,s,x,y);
1836 n=n+coeffn*y;
1837 }
1838
solarn(double & n)1839 void Moon200::solarn (double& n)
1840 {
1841 // perturbation N of ecliptic latitude
1842 double x, y;
1843
1844 n = 0.0;
1845 addn (-526.069, 0, 0, 1, -2, n, x, y);
1846 addn (-3.352, 0, 0, 1, -4, n, x, y);
1847 addn (44.297, 1, 0, 1, -2, n, x, y);
1848 addn (-6.0, 1, 0, 1, -4, n, x, y);
1849 addn (20.599, -1, 0, 1, 0, n, x, y);
1850 addn (-30.598, -1, 0, 1, -2, n, x, y);
1851 addn (-24.649, -2, 0, 1, 0, n, x, y);
1852 addn (-2.0, -2, 0, 1, -2, n, x, y);
1853 addn (-22.571, 0, 1, 1, -2, n, x, y);
1854 addn (10.985, 0, -1, 1, -2, n, x, y);
1855 }
1856
planetary(double t)1857 void Moon200::planetary (double t)
1858 {
1859 // perturbations of the ecliptic longitude by Venus and Jupiter
1860
1861 dlam = dlam
1862 + 0.82*sinus(0.7736 -62.5512*t)+0.31*sinus(0.0466-125.1025*t)
1863 + 0.35*sinus(0.5785-25.1042*t)+0.66*sinus(0.4591+1335.8075*t)
1864 + 0.64*sinus(0.313-91.568*t)+1.14*sinus(0.148+1331.2898*t)
1865 + 0.21*sinus(0.5918+1056.5859*t)+0.44*sinus(0.5784+1322.8595*t)
1866 + 0.24*sinus(0.2275-5.7374*t)+0.28*sinus(0.2965+2.6929*t)
1867 + 0.33*sinus(0.3132+6.3368*t);
1868 }
1869
1870 /*-------------------- Class Eclipse --------------------------------------*/
1871 /*
1872 Calculalations for Solar and Lunar Eclipses
1873 */
1874
Eclipse()1875 Eclipse::Eclipse()
1876 {
1877 // some assignments just in case
1878 rs.assign(1,0,0);
1879 rm.assign(1,0,0);
1880 eshadow.assign (1,0,0);
1881 rint.assign (1,0,0);
1882 d_umbra = 0;
1883 ep2 = 0;
1884 }
1885
solar(double jd,double tdut,double & phi,double & lamda)1886 int Eclipse::solar (double jd, double tdut, double& phi, double& lamda)
1887 {
1888 /* Calculates various items about a solar eclipse.
1889 INPUT:
1890 jd = Modified Julian Date (UT)
1891 tdut = TDT - UT in sec
1892
1893 OUTPUT:
1894 phi, lamda: Geographic latitude and longitude (in radians)
1895 of center of shadow if central eclipse
1896
1897 RETURN:
1898 0 : No eclipse
1899 1 : Partial eclipse
1900 2 : Non-central annular eclipse
1901 3 : Non-central total eclipse
1902 4 : Annular eclipse
1903 5 : Total eclipse
1904 */
1905
1906 const double flat = 0.996633; // flatting of the Earth
1907 const double ds = 218.245445; // diameter of Sun in Earth radii
1908 const double dm = 0.544986; // diameter of Moon in Earth radii
1909 double s0, s, dlt, r2, r0;
1910 int phase;
1911 Vec3 ve;
1912
1913 // get the apparent equatorial coordinates of the sun and the moon
1914 equ_sun_moon(jd, tdut);
1915
1916 rs[2]/=flat; // adjust for flatting of the Earth
1917 rm[2]/=flat;
1918 rint.assign ();
1919 lamda = 0;
1920 phi = 0;
1921
1922 // intersect shadow axis with Earth
1923 eshadow = rm - rs;
1924 eshadow = vnorm(eshadow); // direction vector of shadow
1925
1926 s0 = - dot(rm, eshadow); // distance Moon - fundamental plane
1927 r2 = dot (rm,rm);
1928 dlt = s0*s0 + 1.0 - r2;
1929 r0 = 1.0 - dlt;
1930 if (r0 > 0) r0 = sqrt (r0);
1931 else r0 = 0; // distance center of Earth - shadow axis
1932
1933 r2 = abs(rs - rm);
1934 d_umbra = (ds - dm) * s0 / r2 - dm;// diameter of umbra at fundamental plane
1935 d_penumbra = (ds + dm) * s0 / r2 + dm;
1936
1937 // get phase of eclipse
1938 if (r0 < 1.0)
1939 {
1940 if (dlt > 0) dlt = sqrt(dlt);
1941 else dlt = 0;
1942 s = s0 - dlt; // distance Moon - fundamental plane
1943 d_umbra = (ds - dm) * s / r2 - dm; // diameter of umbra at surface
1944 rint = rm + s * eshadow;
1945 rint[2] *= flat; // vector to intersection
1946 ve = carpol(rint);
1947 lamda = ve[1] - lsidtim(jd,0,ep2)*0.261799387799; // geographic coordinates
1948 if (lamda > M_PI) lamda -= 2.0*M_PI;
1949 if (lamda < (-M_PI)) lamda += 2.0*M_PI;
1950 phi = sqrt(rint[0]*rint[0] + rint[1]*rint[1])*0.993305615;
1951 phi = atan2(rint[2],phi);
1952
1953 if (d_umbra > 0) phase = 4; // central annular eclipse
1954 else phase = 5; // central total eclipse
1955 }
1956 else
1957 {
1958 if (r0 < (1.0 + 0.5 * fabs(d_umbra)))
1959 {
1960 if (d_umbra > 0) phase = 2; // non-central annular eclipse
1961 else phase = 3; // non-central total eclipse
1962 }
1963 else
1964 {
1965 if (r0 < (1.0 + 0.5*d_penumbra)) phase = 1; // partial eclipse
1966 else phase = 0; // no eclipse
1967 }
1968 }
1969
1970 rs[2]*=flat; // restore from flatting of the Earth
1971 rm[2]*=flat;
1972
1973 return phase;
1974 }
1975
maxpos(double jd,double tdut,double & phi,double & lamda)1976 void Eclipse::maxpos (double jd, double tdut, double& phi, double& lamda)
1977 {
1978 /* Calculates the geographic position of the place of maximum eclipse
1979 at the time jd for a non-central eclipse. (NOTE that in case of
1980 a central eclipse the maximum position will be calculated by
1981 function solar. Do not use this function in that case as the
1982 result would be incorrect! No check is being done in this routine
1983 whether an eclipse is visible at all. Use maxpos only in case of a
1984 confirmed partial or non-central eclipse.
1985
1986 INPUT:
1987 jd = Modified Julian Date (UT)
1988 tdut = TDT - UT in sec
1989
1990 OUTPUT:
1991 phi, lamda: Geographic latitude and longitude (in radians)
1992 of maximum eclipse at the time
1993 */
1994
1995 const double flat = 0.996633; // flatting of the Earth
1996 double s0;
1997 Vec3 ve;
1998
1999 // get the apparent equatorial coordinates of the sun and the moon
2000 equ_sun_moon(jd, tdut);
2001
2002 rs[2]/=flat; // adjust for flatting of the Earth
2003 rm[2]/=flat;
2004 rint.assign ();
2005 lamda = 0;
2006 phi = 0;
2007
2008 // intersect shadow axis with Earth
2009 eshadow = rm - rs;
2010 eshadow = vnorm(eshadow); // direction vector of shadow
2011
2012 s0 = - dot(rm, eshadow); // distance Moon - fundamental plane
2013 rint = rm + s0 * eshadow;
2014 rint = vnorm(rint); // normalize to 1 Earth radius
2015 rint[2] *= flat; // vector to position of maximum eclipse
2016 ve = carpol(rint);
2017
2018 lamda = ve[1] - lsidtim(jd,0,ep2)*0.261799387799; // geographic coordinates
2019 if (lamda > M_PI) lamda -= 2.0*M_PI;
2020 if (lamda < (-M_PI)) lamda += 2.0*M_PI;
2021 phi = sqrt(rint[0]*rint[0] + rint[1]*rint[1])*0.993305615;
2022 phi = atan2(rint[2],phi);
2023
2024 rs[2]*=flat; // restore from flatting of the Earth
2025 rm[2]*=flat;
2026 }
2027
penumd(double jd,double tdut,Vec3 & vrm,Vec3 & ves,double & dpn,double & pang)2028 void Eclipse::penumd (double jd, double tdut, Vec3& vrm, Vec3& ves,
2029 double& dpn, double& pang)
2030 {
2031 /* Calculates various items needed for finding out the northern and
2032 southern border of the penumbra during a solar eclipse.
2033
2034 INPUT:
2035 jd = Modified Julian Date (UT)
2036 tdut = TDT - UT in sec
2037
2038 OUTPUT:
2039 vrm = position vector of Moon adjusted for flattening
2040 ves = unit vector pointing into the direction of the center of the shadow
2041 dpn = diameter of the penumbra at the fundamental plane
2042 pang = angle of penumbral half-cone in radians
2043 */
2044
2045 const double flat = 0.996633; // flatting of the Earth
2046 const double ds = 218.245445; // diameter of Sun in Earth radii
2047 const double dm = 0.544986; // diameter of Moon in Earth radii
2048 double s0, r2;
2049
2050 // get the apparent equatorial coordinates of the sun and the moon
2051 equ_sun_moon(jd, tdut);
2052 rs[2]/=flat; // adjust for flatting of the Earth
2053 rm[2]/=flat;
2054
2055 // intersect shadow axis with Earth
2056 eshadow = rm - rs;
2057 pang = abs(eshadow);
2058 eshadow = vnorm(eshadow); // direction vector of shadow
2059 ves = eshadow;
2060 vrm = rm;
2061
2062 s0 = - dot(rm, eshadow); // distance Moon - fundamental plane
2063
2064 r2 = abs(rs - rm);
2065 dpn = (ds + dm) * s0 / r2 + dm;
2066
2067 // penumbral angle
2068 pang = asin((dm + ds) / (2.0 * pang));
2069
2070 rs[2]*=flat; // restore from flatting of the Earth
2071 rm[2]*=flat;
2072 }
2073
umbra(double jd,double tdut,Vec3 & vrm,Vec3 & ves,double & dpn,double & pang)2074 void Eclipse::umbra (double jd, double tdut, Vec3& vrm, Vec3& ves,
2075 double& dpn, double& pang)
2076 {
2077 /* Calculates various items needed for finding out the northern and
2078 southern border of the umbra during a solar eclipse.
2079
2080 INPUT:
2081 jd = Modified Julian Date (UT)
2082 tdut = TDT - UT in sec
2083
2084 OUTPUT:
2085 vrm = position vector of Moon adjusted for flattening
2086 ves = unit vector pointing into the direction of the center of the shadow
2087 dpn = diameter of the umbra at the fundamental plane
2088 pang = angle of umbral half-cone in radians
2089 */
2090
2091 const double flat = 0.996633; // flatting of the Earth
2092 const double ds = 218.245445; // diameter of Sun in Earth radii
2093 const double dm = 0.544986; // diameter of Moon in Earth radii
2094 double s0, r2;
2095
2096 // get the apparent equatorial coordinates of the sun and the moon
2097 equ_sun_moon(jd, tdut);
2098 rs[2]/=flat; // adjust for flatting of the Earth
2099 rm[2]/=flat;
2100
2101 // intersect shadow axis with Earth
2102 eshadow = rm - rs;
2103 pang = abs(eshadow);
2104 eshadow = vnorm(eshadow); // direction vector of shadow
2105 ves = eshadow;
2106 vrm = rm;
2107
2108 s0 = - dot(rm, eshadow); // distance Moon - fundamental plane
2109
2110 r2 = abs(rs - rm);
2111 dpn = (ds - dm) * s0 / r2 - dm;
2112
2113 // umbral angle
2114 pang = asin((ds - dm) / (2.0 * pang));
2115
2116 rs[2]*=flat; // restore from flatting of the Earth
2117 rm[2]*=flat;
2118 }
2119
equ_sun_moon(double jd,double tdut)2120 void Eclipse::equ_sun_moon(double jd, double tdut)
2121 {
2122 /* Get the equatorial coordinates of the Sun and the Moon and
2123 store them in rs and rm.
2124 Also store the time t in Julian Centuries and the correction
2125 ep2 for the Apparent Sidereal Time.
2126 jd = Modified Julian Date (UT)
2127 tdut = TDT - UT in sec
2128 */
2129 double ae = 23454.77992; // 149597870.0/6378.14 = 1AE -> Earth Radii
2130 Mat3 mx;
2131
2132 t = julcent (jd) + tdut / 3.15576e9; // =(86400.0 * 36525.0);
2133 rs = sun.position(t);
2134 rm = moon.position(t);
2135 rs = eclequ(t,rs);
2136 rm = eclequ(t,rm);
2137
2138 // correct Moon coordinates for center of figure
2139 mx = zrot(-2.4240684e-6); // +0.5" in longitude
2140 rm = mxvct(mx,rm);
2141 mx = yrot(-1.2120342e-6); // -0.25" in latitude
2142 rm = mxvct(mx,rm);
2143 mx = nutmat(t,ep2); // nutation
2144 rs = mxvct(mx,rs); // apparent coordinates
2145 rs = aberrat(t,rs);
2146 rs *= ae;
2147 rm = mxvct(mx,rm);
2148 }
2149
duration(double jd,double tdut,double & width)2150 double Eclipse::duration (double jd, double tdut, double& width)
2151 {
2152 /* Get the duration of a central eclipse at the center point
2153 for MJD jd and TDT-UT tdut in seconds.
2154 Also return the width of the umbra in km.
2155 A call to solar with the respective jd must have been done
2156 prior to calling this function !!!
2157 */
2158 const double omega = 4.3755e-3; // radial velocity of Earth in rad/min.
2159
2160 double dur, lm, pa, umbold;
2161 Vec3 rold, eold, rsold, rmold;
2162 Mat3 mx;
2163
2164 // save old values for jd
2165 rold = rint;
2166 eold = eshadow;
2167 umbold = d_umbra;
2168 rsold = rs;
2169 rmold = rm;
2170
2171 dur = 0.1; // 0.1 min
2172 if (solar(jd+dur/1440.0,tdut, lm, pa) > 3)
2173 {
2174 mx = zrot(dur*omega);
2175 rint = mxvct(mx,rint);
2176 rint -= rold;
2177 pa = dot (rint, eold);
2178 lm = dot (rint, rint) - pa*pa;
2179 if (lm > 0) lm = sqrt(lm);
2180 else lm = 0;
2181 if (lm > 0) dur = fabs(umbold) / lm * dur * 60.0;
2182 else dur = 0;
2183 }
2184
2185 else dur = -1;
2186
2187 // restore old values for jd
2188 d_umbra = umbold;
2189 eshadow = eold;
2190 eold = rold*rint; // direction perpendicular of apparent shadow movement
2191 rint = rold;
2192 rs = rsold;
2193 rm = rmold;
2194
2195 // get width of umbra at center location
2196 rold = vnorm (rint);
2197 pa = dot(rold,eshadow);
2198 if (pa > 1.0) pa = 1.0;
2199 if (pa < -1.0) pa = -1.0;
2200 if (fabs(pa) < 1.0e-15) lm = fabs(d_umbra); // just in case
2201 else lm = fabs(d_umbra / pa); // projection along shadow vector
2202
2203 rsold = vnorm(rint*eshadow);
2204 rmold = vnorm(eold);
2205 pa = dot(rsold,rmold); // cos of projection shadow - perperndicular
2206 if (pa > 1.0) pa = 1.0;
2207 if (pa < -1.0) pa = -1.0;
2208 umbold = fabs(sin(acos(pa)));
2209
2210 // this is for very slant shadow angles respective to movement
2211 lm = lm * umbold;
2212 umbold = fabs(pa * d_umbra);
2213 if (umbold > lm) width = umbold;
2214 else width = lm;
2215
2216 // this is the normal case
2217 rold = vnorm (eold);
2218 pa = dot(rold,eshadow);
2219 if (pa > 1.0) pa = 1.0;
2220 if (pa < -1.0) pa = -1.0;
2221 pa = fabs(sin(acos(pa)));
2222 if (pa < 0.00001) pa = 0.00001;
2223 lm = d_umbra / pa; // allow for projection
2224 lm = fabs(lm);
2225
2226 if (lm > width) width = lm;
2227 width = width * 6378.14;
2228
2229 return dur;
2230 }
2231
GetRSun()2232 Vec3 Eclipse::GetRSun () // get Earth - Sun vector in Earth radii
2233 {
2234 return rs;
2235 }
2236
GetRMoon()2237 Vec3 Eclipse::GetRMoon () // get Earth - Moon vector in Earth radii
2238 {
2239 return rm;
2240 }
2241
GetEp2()2242 double Eclipse::GetEp2 () // get the ep2 value
2243 {
2244 return ep2;
2245 }
2246
lunar(double jd,double tdut)2247 int Eclipse::lunar (double jd, double tdut)
2248 {
2249 /* check whether lunar eclipse is in progress at
2250 Modified Julian Date jd (UT).
2251 tdut = TDT - UT in sec
2252
2253 RETURN:
2254 0 : No eclipse
2255 1 : Partial Penumbral Eclipse
2256 2 : Total Penumbral Eclipse
2257 3 : Partial Umbral Eclipse
2258 4 : Total Umbral Eclipse
2259 */
2260
2261 const double dm = 0.544986; // diameter of Moon in Earth radii
2262 const double ds = 218.245445; // diameter of Sun in Earth radii
2263
2264 double umbra, penumbra;
2265 double r2, s0, sep;
2266 int phase;
2267 Vec3 v1, v2;
2268
2269 // get position of Sun and Moon
2270 equ_sun_moon(jd, tdut);
2271
2272 // get radius of umbra and penumbra
2273 r2 = abs(rs);
2274 s0 = abs (rm);
2275 umbra = 1.02*fabs((ds - 2.0) * s0 / r2 - 2.0)*0.5; // radius of umbra
2276 penumbra = 1.02*fabs((ds + 2.0) * s0 / r2 + 2.0)*0.5;//radius of penumbra
2277 /* (the factor 1.02 allows for enlargment of shadow due to
2278 Earth's atmosphere) */
2279
2280 // get angular separation of center of shadow and Moon
2281 r2 = abs(rm);
2282 sep = dot(rs,rm)/(abs(rs)*r2);
2283 if (fabs(sep) > 1.0) sep = 1.0;
2284 sep = acos(sep); // in radians
2285 sep = fabs(tan(sep)*r2); // distance of Moon and shadow in Earth radii
2286
2287 // Now check the kind of eclipse
2288 if (sep < (umbra - dm/2.0)) phase = 4;
2289 else if (sep < (umbra + dm/2.0)) phase = 3;
2290 else if (sep < (penumbra - dm/2.0)) phase = 2;
2291 else if (sep < (penumbra + dm/2.0)) phase = 1;
2292 else phase = 0;
2293
2294 return phase;
2295 }
2296
2297
2298