1 /* SWISSEPH
2
3 SWISSEPH modules that may be useful for other applications
4 e.g. chopt.c, venus.c, swetest.c
5
6 Authors: Dieter Koch and Alois Treindl, Astrodienst Zurich
7
8 coordinate transformations
9 obliquity of ecliptic
10 nutation
11 precession
12 delta t
13 sidereal time
14 setting or getting of tidal acceleration of moon
15 chebyshew interpolation
16 ephemeris file name generation
17 cyclic redundancy checksum CRC
18 modulo and normalization functions
19
20 **************************************************************/
21 /* Copyright (C) 1997 - 2021 Astrodienst AG, Switzerland. All rights reserved.
22
23 License conditions
24 ------------------
25
26 This file is part of Swiss Ephemeris.
27
28 Swiss Ephemeris is distributed with NO WARRANTY OF ANY KIND. No author
29 or distributor accepts any responsibility for the consequences of using it,
30 or for whether it serves any particular purpose or works at all, unless he
31 or she says so in writing.
32
33 Swiss Ephemeris is made available by its authors under a dual licensing
34 system. The software developer, who uses any part of Swiss Ephemeris
35 in his or her software, must choose between one of the two license models,
36 which are
37 a) GNU Affero General Public License (AGPL)
38 b) Swiss Ephemeris Professional License
39
40 The choice must be made before the software developer distributes software
41 containing parts of Swiss Ephemeris to others, and before any public
42 service using the developed software is activated.
43
44 If the developer choses the AGPL software license, he or she must fulfill
45 the conditions of that license, which includes the obligation to place his
46 or her whole software project under the AGPL or a compatible license.
47 See https://www.gnu.org/licenses/agpl-3.0.html
48
49 If the developer choses the Swiss Ephemeris Professional license,
50 he must follow the instructions as found in http://www.astro.com/swisseph/
51 and purchase the Swiss Ephemeris Professional Edition from Astrodienst
52 and sign the corresponding license contract.
53
54 The License grants you the right to use, copy, modify and redistribute
55 Swiss Ephemeris, but only under certain conditions described in the License.
56 Among other things, the License requires that the copyright notices and
57 this notice be preserved on all copies.
58
59 Authors of the Swiss Ephemeris: Dieter Koch and Alois Treindl
60
61 The authors of Swiss Ephemeris have no control or influence over any of
62 the derived works, i.e. over software or services created by other
63 programmers which use Swiss Ephemeris functions.
64
65 The names of the authors or of the copyright holder (Astrodienst) must not
66 be used for promoting any software, product or service which uses or contains
67 the Swiss Ephemeris. This copyright notice is the ONLY place where the
68 names of the authors can legally appear, except in cases where they have
69 given special permission in writing.
70
71 The trademarks 'Swiss Ephemeris' and 'Swiss Ephemeris inside' may be used
72 for promoting such software, products or services.
73 */
74
75 #include <string.h>
76 #include <ctype.h>
77 #include "swephexp.h"
78 #include "sweph.h"
79 #include "swephlib.h"
80 #if MSDOS
81 # include <process.h>
82 # define strdup _strdup
83 #endif
84
85 #ifdef TRACE
86 void swi_open_trace(char *serr);
87 TLS FILE *swi_fp_trace_c = NULL;
88 TLS FILE *swi_fp_trace_out = NULL;
89 TLS int32 swi_trace_count = 0;
90 #endif
91
92 static void init_crc32(void);
93 static int init_dt(void);
94 static double adjust_for_tidacc(double ans, double Y, double tid_acc, double tid_acc0, AS_BOOL adjust_after_1955);
95 static double deltat_espenak_meeus_1620(double tjd, double tid_acc);
96 static double deltat_stephenson_etc_2016(double tjd, double tid_acc);
97 static double deltat_longterm_morrison_stephenson(double tjd);
98 static double deltat_stephenson_morrison_2004_1600(double tjd, double tid_acc);
99 static double deltat_stephenson_morrison_1997_1600(double tjd, double tid_acc);
100 static double deltat_aa(double tjd, double tid_acc);
101
102 #define SEFLG_EPHMASK (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH)
103
104 /* Reduce x modulo 360 degrees
105 */
swe_degnorm(double x)106 double CALL_CONV swe_degnorm(double x)
107 {
108 double y;
109 y = fmod(x, 360.0);
110 if (fabs(y) < 1e-13) y = 0; /* Alois fix 11-dec-1999 */
111 if( y < 0.0 ) y += 360.0;
112 return(y);
113 }
114
115 /* Reduce x modulo TWOPI degrees
116 */
swe_radnorm(double x)117 double CALL_CONV swe_radnorm(double x)
118 {
119 double y;
120 y = fmod(x, TWOPI);
121 if (fabs(y) < 1e-13) y = 0; /* Alois fix 11-dec-1999 */
122 if( y < 0.0 ) y += TWOPI;
123 return(y);
124 }
125
swe_deg_midp(double x1,double x0)126 double CALL_CONV swe_deg_midp(double x1, double x0)
127 {
128 double d, y;
129 d = swe_difdeg2n(x1, x0); /* arc from x0 to x1 */
130 y = swe_degnorm(x0 + d / 2);
131 return(y);
132 }
133
swe_rad_midp(double x1,double x0)134 double CALL_CONV swe_rad_midp(double x1, double x0)
135 {
136 return DEGTORAD * swe_deg_midp(x1 * RADTODEG, x0 * RADTODEG);
137 }
138
139 /* Reduce x modulo 2*PI
140 */
swi_mod2PI(double x)141 double swi_mod2PI(double x)
142 {
143 double y;
144 y = fmod(x, TWOPI);
145 if( y < 0.0 ) y += TWOPI;
146 return(y);
147 }
148
149
swi_angnorm(double x)150 double swi_angnorm(double x)
151 {
152 if (x < 0.0 )
153 return x + TWOPI;
154 else if (x >= TWOPI)
155 return x - TWOPI;
156 else
157 return x;
158 }
159
swi_cross_prod(double * a,double * b,double * x)160 void swi_cross_prod(double *a, double *b, double *x)
161 {
162 x[0] = a[1]*b[2] - a[2]*b[1];
163 x[1] = a[2]*b[0] - a[0]*b[2];
164 x[2] = a[0]*b[1] - a[1]*b[0];
165 }
166
167 /* Evaluates a given chebyshev series coef[0..ncf-1]
168 * with ncf terms at x in [-1,1]. Communications of the ACM, algorithm 446,
169 * April 1973 (vol. 16 no.4) by Dr. Roger Broucke.
170 */
swi_echeb(double x,double * coef,int ncf)171 double swi_echeb(double x, double *coef, int ncf)
172 {
173 int j;
174 double x2, br, brp2, brpp;
175 x2 = x * 2.;
176 br = 0.;
177 brp2 = 0.; /* dummy assign to silence gcc warning */
178 brpp = 0.;
179 for (j = ncf - 1; j >= 0; j--) {
180 brp2 = brpp;
181 brpp = br;
182 br = x2 * brpp - brp2 + coef[j];
183 }
184 return (br - brp2) * .5;
185 }
186
187 /*
188 * evaluates derivative of chebyshev series, see echeb
189 */
swi_edcheb(double x,double * coef,int ncf)190 double swi_edcheb(double x, double *coef, int ncf)
191 {
192 double bjpl, xjpl;
193 int j;
194 double x2, bf, bj, dj, xj, bjp2, xjp2;
195 x2 = x * 2.;
196 bf = 0.; /* dummy assign to silence gcc warning */
197 bj = 0.; /* dummy assign to silence gcc warning */
198 xjp2 = 0.;
199 xjpl = 0.;
200 bjp2 = 0.;
201 bjpl = 0.;
202 for (j = ncf - 1; j >= 1; j--) {
203 dj = (double) (j + j);
204 xj = coef[j] * dj + xjp2;
205 bj = x2 * bjpl - bjp2 + xj;
206 bf = bjp2;
207 bjp2 = bjpl;
208 bjpl = bj;
209 xjp2 = xjpl;
210 xjpl = xj;
211 }
212 return (bj - bf) * .5;
213 }
214
215 /*
216 * conversion between ecliptical and equatorial polar coordinates.
217 * for users of SWISSEPH, not used by our routines.
218 * for ecl. to equ. eps must be negative.
219 * for equ. to ecl. eps must be positive.
220 * xpo, xpn are arrays of 3 doubles containing position.
221 * attention: input must be in degrees!
222 */
swe_cotrans(double * xpo,double * xpn,double eps)223 void CALL_CONV swe_cotrans(double *xpo, double *xpn, double eps)
224 {
225 int i;
226 double x[6], e = eps * DEGTORAD;
227 for(i = 0; i <= 1; i++)
228 x[i] = xpo[i];
229 x[0] *= DEGTORAD;
230 x[1] *= DEGTORAD;
231 x[2] = 1;
232 for(i = 3; i <= 5; i++)
233 x[i] = 0;
234 swi_polcart(x, x);
235 swi_coortrf(x, x, e);
236 swi_cartpol(x, x);
237 xpn[0] = x[0] * RADTODEG;
238 xpn[1] = x[1] * RADTODEG;
239 xpn[2] = xpo[2];
240 }
241
242 /*
243 * conversion between ecliptical and equatorial polar coordinates
244 * with speed.
245 * for users of SWISSEPH, not used by our routines.
246 * for ecl. to equ. eps must be negative.
247 * for equ. to ecl. eps must be positive.
248 * xpo, xpn are arrays of 6 doubles containing position and speed.
249 * attention: input must be in degrees!
250 */
swe_cotrans_sp(double * xpo,double * xpn,double eps)251 void CALL_CONV swe_cotrans_sp(double *xpo, double *xpn, double eps)
252 {
253 int i;
254 double x[6], e = eps * DEGTORAD;
255 for (i = 0; i <= 5; i++)
256 x[i] = xpo[i];
257 x[0] *= DEGTORAD;
258 x[1] *= DEGTORAD;
259 x[2] = 1; /* avoids problems with polcart(), if x[2] = 0 */
260 x[3] *= DEGTORAD;
261 x[4] *= DEGTORAD;
262 swi_polcart_sp(x, x);
263 swi_coortrf(x, x, e);
264 swi_coortrf(x+3, x+3, e);
265 swi_cartpol_sp(x, xpn);
266 xpn[0] *= RADTODEG;
267 xpn[1] *= RADTODEG;
268 xpn[2] = xpo[2];
269 xpn[3] *= RADTODEG;
270 xpn[4] *= RADTODEG;
271 xpn[5] = xpo[5];
272 }
273
274 /*
275 * conversion between ecliptical and equatorial cartesian coordinates
276 * for ecl. to equ. eps must be negative
277 * for equ. to ecl. eps must be positive
278 */
swi_coortrf(double * xpo,double * xpn,double eps)279 void swi_coortrf(double *xpo, double *xpn, double eps)
280 {
281 double sineps, coseps;
282 double x[3];
283 sineps = sin(eps);
284 coseps = cos(eps);
285 x[0] = xpo[0];
286 x[1] = xpo[1] * coseps + xpo[2] * sineps;
287 x[2] = -xpo[1] * sineps + xpo[2] * coseps;
288 xpn[0] = x[0];
289 xpn[1] = x[1];
290 xpn[2] = x[2];
291 }
292
293 /*
294 * conversion between ecliptical and equatorial cartesian coordinates
295 * sineps sin(eps)
296 * coseps cos(eps)
297 * for ecl. to equ. sineps must be -sin(eps)
298 */
swi_coortrf2(double * xpo,double * xpn,double sineps,double coseps)299 void swi_coortrf2(double *xpo, double *xpn, double sineps, double coseps)
300 {
301 double x[3];
302 x[0] = xpo[0];
303 x[1] = xpo[1] * coseps + xpo[2] * sineps;
304 x[2] = -xpo[1] * sineps + xpo[2] * coseps;
305 xpn[0] = x[0];
306 xpn[1] = x[1];
307 xpn[2] = x[2];
308 }
309
310 /* conversion of cartesian (x[3]) to polar coordinates (l[3]).
311 * x = l is allowed.
312 * if |x| = 0, then lon, lat and rad := 0.
313 */
swi_cartpol(double * x,double * l)314 void swi_cartpol(double *x, double *l)
315 {
316 double rxy;
317 double ll[3];
318 if (x[0] == 0 && x[1] == 0 && x[2] == 0) {
319 l[0] = l[1] = l[2] = 0;
320 return;
321 }
322 rxy = x[0]*x[0] + x[1]*x[1];
323 ll[2] = sqrt(rxy + x[2]*x[2]);
324 rxy = sqrt(rxy);
325 ll[0] = atan2(x[1], x[0]);
326 if (ll[0] < 0.0) ll[0] += TWOPI;
327 if (rxy == 0) {
328 if (x[2] >= 0)
329 ll[1] = PI / 2;
330 else
331 ll[1] = -(PI / 2);
332 } else {
333 ll[1] = atan(x[2] / rxy);
334 }
335 l[0] = ll[0];
336 l[1] = ll[1];
337 l[2] = ll[2];
338 }
339
340 /* conversion from polar (l[3]) to cartesian coordinates (x[3]).
341 * x = l is allowed.
342 */
swi_polcart(double * l,double * x)343 void swi_polcart(double *l, double *x)
344 {
345 double xx[3];
346 double cosl1;
347 cosl1 = cos(l[1]);
348 xx[0] = l[2] * cosl1 * cos(l[0]);
349 xx[1] = l[2] * cosl1 * sin(l[0]);
350 xx[2] = l[2] * sin(l[1]);
351 x[0] = xx[0];
352 x[1] = xx[1];
353 x[2] = xx[2];
354 }
355
356 /* conversion of position and speed.
357 * from cartesian (x[6]) to polar coordinates (l[6]).
358 * x = l is allowed.
359 * if position is 0, function returns direction of
360 * motion.
361 */
swi_cartpol_sp(double * x,double * l)362 void swi_cartpol_sp(double *x, double *l)
363 {
364 int i;
365 double xx[6], ll[6];
366 double rxy, coslon, sinlon, coslat, sinlat;
367 /* zero position */
368 if (x[0] == 0 && x[1] == 0 && x[2] == 0) {
369 ll[0] = ll[1] = ll[3] = ll[4] = 0;
370 ll[5] = sqrt(square_sum((x+3)));
371 swi_cartpol(x+3, ll);
372 ll[2] = 0;
373 for (i = 0; i <= 5; i++)
374 l[i] = ll[i];
375 return;
376 }
377 /* zero speed */
378 if (x[3] == 0 && x[4] == 0 && x[5] == 0) {
379 l[3] = l[4] = l[5] = 0;
380 swi_cartpol(x, l);
381 return;
382 }
383 /* position */
384 rxy = x[0]*x[0] + x[1]*x[1];
385 ll[2] = sqrt(rxy + x[2]*x[2]);
386 rxy = sqrt(rxy);
387 ll[0] = atan2(x[1], x[0]);
388 if (ll[0] < 0.0) ll[0] += TWOPI;
389 ll[1] = atan(x[2] / rxy);
390 /* speed:
391 * 1. rotate coordinate system by longitude of position about z-axis,
392 * so that new x-axis = position radius projected onto x-y-plane.
393 * in the new coordinate system
394 * vy'/r = dlong/dt, where r = sqrt(x^2 +y^2).
395 * 2. rotate coordinate system by latitude about new y-axis.
396 * vz"/r = dlat/dt, where r = position radius.
397 * vx" = dr/dt
398 */
399 coslon = x[0] / rxy; /* cos(l[0]); */
400 sinlon = x[1] / rxy; /* sin(l[0]); */
401 coslat = rxy / ll[2]; /* cos(l[1]); */
402 sinlat = x[2] / ll[2]; /* sin(ll[1]); */
403 xx[3] = x[3] * coslon + x[4] * sinlon;
404 xx[4] = -x[3] * sinlon + x[4] * coslon;
405 l[3] = xx[4] / rxy; /* speed in longitude */
406 xx[4] = -sinlat * xx[3] + coslat * x[5];
407 xx[5] = coslat * xx[3] + sinlat * x[5];
408 l[4] = xx[4] / ll[2]; /* speed in latitude */
409 l[5] = xx[5]; /* speed in radius */
410 l[0] = ll[0]; /* return position */
411 l[1] = ll[1];
412 l[2] = ll[2];
413 }
414
415 /* conversion of position and speed
416 * from polar (l[6]) to cartesian coordinates (x[6])
417 * x = l is allowed
418 * explanation s. swi_cartpol_sp()
419 */
swi_polcart_sp(double * l,double * x)420 void swi_polcart_sp(double *l, double *x)
421 {
422 double sinlon, coslon, sinlat, coslat;
423 double xx[6], rxy, rxyz;
424 /* zero speed */
425 if (l[3] == 0 && l[4] == 0 && l[5] == 0) {
426 x[3] = x[4] = x[5] = 0;
427 swi_polcart(l, x);
428 return;
429 }
430 /* position */
431 coslon = cos(l[0]);
432 sinlon = sin(l[0]);
433 coslat = cos(l[1]);
434 sinlat = sin(l[1]);
435 xx[0] = l[2] * coslat * coslon;
436 xx[1] = l[2] * coslat * sinlon;
437 xx[2] = l[2] * sinlat;
438 /* speed; explanation s. swi_cartpol_sp(), same method the other way round*/
439 rxyz = l[2];
440 rxy = sqrt(xx[0] * xx[0] + xx[1] * xx[1]);
441 xx[5] = l[5];
442 xx[4] = l[4] * rxyz;
443 x[5] = sinlat * xx[5] + coslat * xx[4]; /* speed z */
444 xx[3] = coslat * xx[5] - sinlat * xx[4];
445 xx[4] = l[3] * rxy;
446 x[3] = coslon * xx[3] - sinlon * xx[4]; /* speed x */
447 x[4] = sinlon * xx[3] + coslon * xx[4]; /* speed y */
448 x[0] = xx[0]; /* return position */
449 x[1] = xx[1];
450 x[2] = xx[2];
451 }
452
swi_dot_prod_unit(double * x,double * y)453 double swi_dot_prod_unit(double *x, double *y)
454 {
455 double dop = x[0]*y[0]+x[1]*y[1]+x[2]*y[2];
456 double e1 = sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2]);
457 double e2 = sqrt(y[0]*y[0]+y[1]*y[1]+y[2]*y[2]);
458 dop /= e1;
459 dop /= e2;
460 if (dop > 1)
461 dop = 1;
462 if (dop < -1)
463 dop = -1;
464 return dop;
465 }
466
467 /* functions for precession and ecliptic obliquity according to Vondrák et alii, 2011 */
468 #define AS2R (DEGTORAD / 3600.0)
469 #define D2PI TWOPI
470 #define EPS0 (84381.406 * AS2R)
471 #define NPOL_PEPS 4
472 #define NPER_PEPS 10
473 #define NPOL_PECL 4
474 #define NPER_PECL 8
475 #define NPOL_PEQU 4
476 #define NPER_PEQU 14
477
478 /* for pre_peps(): */
479 /* polynomials */
480 static const double pepol[NPOL_PEPS][2] = {
481 {+8134.017132, +84028.206305},
482 {+5043.0520035, +0.3624445},
483 {-0.00710733, -0.00004039},
484 {+0.000000271, -0.000000110}
485 };
486
487 /* periodics */
488 static const double peper[5][NPER_PEPS] = {
489 {+409.90, +396.15, +537.22, +402.90, +417.15, +288.92, +4043.00, +306.00, +277.00, +203.00},
490 {-6908.287473, -3198.706291, +1453.674527, -857.748557, +1173.231614, -156.981465, +371.836550, -216.619040, +193.691479, +11.891524},
491 {+753.872780, -247.805823, +379.471484, -53.880558, -90.109153, -353.600190, -63.115353, -28.248187, +17.703387, +38.911307},
492 {-2845.175469, +449.844989, -1255.915323, +886.736783, +418.887514, +997.912441, -240.979710, +76.541307, -36.788069, -170.964086},
493 {-1704.720302, -862.308358, +447.832178, -889.571909, +190.402846, -56.564991, -296.222622, -75.859952, +67.473503, +3.014055}
494 };
495
496 /* for pre_pecl(): */
497 /* polynomials */
498 static const double pqpol[NPOL_PECL][2] = {
499 {+5851.607687, -1600.886300},
500 {-0.1189000, +1.1689818},
501 {-0.00028913, -0.00000020},
502 {+0.000000101, -0.000000437}
503 };
504
505 /* periodics */
506 static const double pqper[5][NPER_PECL] = {
507 {708.15, 2309, 1620, 492.2, 1183, 622, 882, 547},
508 {-5486.751211, -17.127623, -617.517403, 413.44294, 78.614193, -180.732815, -87.676083, 46.140315},
509 // original publication A&A 534, A22 (2011):
510 //{-684.66156, 2446.28388, 399.671049, -356.652376, -186.387003, -316.80007, 198.296071, 101.135679},
511 // typo fixed according to A&A 541, C1 (2012)
512 {-684.66156, 2446.28388, 399.671049, -356.652376, -186.387003, -316.80007, 198.296701, 101.135679},
513 {667.66673, -2354.886252, -428.152441, 376.202861, 184.778874, 335.321713, -185.138669, -120.97283},
514 {-5523.863691, -549.74745, -310.998056, 421.535876, -36.776172, -145.278396, -34.74445, 22.885731}
515 };
516
517 /* for pre_pequ(): */
518 /* polynomials */
519 static const double xypol[NPOL_PEQU][2] = {
520 {+5453.282155, -73750.930350},
521 {+0.4252841, -0.7675452},
522 {-0.00037173, -0.00018725},
523 {-0.000000152, +0.000000231}
524 };
525
526 /* periodics */
527 static const double xyper[5][NPER_PEQU] = {
528 {256.75, 708.15, 274.2, 241.45, 2309, 492.2, 396.1, 288.9, 231.1, 1610, 620, 157.87, 220.3, 1200},
529 {-819.940624, -8444.676815, 2600.009459, 2755.17563, -167.659835, 871.855056, 44.769698, -512.313065, -819.415595, -538.071099, -189.793622, -402.922932, 179.516345, -9.814756},
530 {75004.344875, 624.033993, 1251.136893, -1102.212834, -2660.66498, 699.291817, 153.16722, -950.865637, 499.754645, -145.18821, 558.116553, -23.923029, -165.405086, 9.344131},
531 {81491.287984, 787.163481, 1251.296102, -1257.950837, -2966.79973, 639.744522, 131.600209, -445.040117, 584.522874, -89.756563, 524.42963, -13.549067, -210.157124, -44.919798},
532 {1558.515853, 7774.939698, -2219.534038, -2523.969396, 247.850422, -846.485643, -1393.124055, 368.526116, 749.045012, 444.704518, 235.934465, 374.049623, -171.33018, -22.899655}
533 };
534
swi_ldp_peps(double tjd,double * dpre,double * deps)535 void swi_ldp_peps(double tjd, double *dpre, double *deps)
536 {
537 int i;
538 int npol = NPOL_PEPS;
539 int nper = NPER_PEPS;
540 double t, p, q, w, a, s, c;
541 t = (tjd - J2000) / 36525.0;
542 p = 0;
543 q = 0;
544 /* periodic terms */
545 for (i = 0; i < nper; i++) {
546 w = D2PI * t;
547 a = w / peper[0][i];
548 s = sin(a);
549 c = cos(a);
550 p += c * peper[1][i] + s * peper[3][i];
551 q += c * peper[2][i] + s * peper[4][i];
552 }
553 /* polynomial terms */
554 w = 1;
555 for (i = 0; i < npol; i++) {
556 p += pepol[i][0] * w;
557 q += pepol[i][1] * w;
558 w *= t;
559 }
560 /* both to radians */
561 p *= AS2R;
562 q *= AS2R;
563 /* return */
564 if (dpre != NULL)
565 *dpre = p;
566 if (deps != NULL)
567 *deps = q;
568 //fprintf(stderr, "%.17f\n", *deps / DEGTORAD);
569 }
570
571 /*
572 * Long term high precision precession,
573 * according to Vondrak/Capitaine/Wallace, "New precession expressions, valid
574 * for long time intervals", in A&A 534, A22(2011).
575 */
576 /* precession of the ecliptic */
pre_pecl(double tjd,double * vec)577 static void pre_pecl(double tjd, double *vec)
578 {
579 int i;
580 int npol = NPOL_PECL;
581 int nper = NPER_PECL;
582 double t, p, q, w, a, s, c, z;
583 t = (tjd - J2000) / 36525.0;
584 p = 0;
585 q = 0;
586 /* periodic terms */
587 for (i = 0; i < nper; i++) {
588 w = D2PI * t;
589 a = w / pqper[0][i];
590 s = sin(a);
591 c = cos(a);
592 p += c * pqper[1][i] + s * pqper[3][i];
593 q += c * pqper[2][i] + s * pqper[4][i];
594 }
595 /* polynomial terms */
596 w = 1;
597 for (i = 0; i < npol; i++) {
598 p += pqpol[i][0] * w;
599 q += pqpol[i][1] * w;
600 w *= t;
601 }
602 /* both to radians */
603 p *= AS2R;
604 q *= AS2R;
605 /* ecliptic pole vector */
606 z = 1 - p * p - q * q;
607 if (z < 0)
608 z = 0;
609 else
610 z = sqrt(z);
611 s = sin(EPS0);
612 c = cos(EPS0);
613 vec[0] = p;
614 vec[1] = - q * c - z * s;
615 vec[2] = - q * s + z * c;
616 }
617
618 /* precession of the equator */
pre_pequ(double tjd,double * veq)619 static void pre_pequ(double tjd, double *veq)
620 {
621 int i;
622 int npol = NPOL_PEQU;
623 int nper = NPER_PEQU;
624 double t, x, y, w, a, s, c;
625 t = (tjd - J2000) / 36525.0;
626 x = 0;
627 y = 0;
628 for (i = 0; i < nper; i++) {
629 w = D2PI * t;
630 a = w / xyper[0][i];
631 s = sin(a);
632 c = cos(a);
633 x += c * xyper[1][i] + s * xyper[3][i];
634 y += c * xyper[2][i] + s * xyper[4][i];
635 }
636 /* polynomial terms */
637 w = 1;
638 for (i = 0; i < npol; i++) {
639 x += xypol[i][0] * w;
640 y += xypol[i][1] * w;
641 w *= t;
642 }
643 x *= AS2R;
644 y *= AS2R;
645 /* equator pole vector */
646 veq[0] = x;
647 veq[1] = y;
648 w = x * x + y * y;
649 if (w < 1)
650 veq[2] = sqrt(1 - w);
651 else
652 veq[2] = 0;
653 }
654
655 #if 0
656 static void swi_cross_prod(double *a, double *b, double *x)
657 {
658 x[0] = a[1] * b[2] - a[2] * b[1];
659 x[1] = a[2] * b[0] - a[0] * b[2];
660 x[2] = a[0] * b[1] - a[1] * b[0];
661 }
662 #endif
663
664 /* precession matrix */
pre_pmat(double tjd,double * rp)665 static void pre_pmat(double tjd, double *rp)
666 {
667 double peqr[3], pecl[3], v[3], w, eqx[3];
668 //tjd = 1219339.078000;
669 /*equator pole */
670 pre_pequ(tjd, peqr);
671 /* ecliptic pole */
672 pre_pecl(tjd, pecl);
673 // fprintf(stderr, "%.17f %.17f %.17f\n", peqr[0], peqr[1], peqr[2]);
674 // fprintf(stderr, "%.17f %.17f %.17f\n", pecl[0], pecl[1], pecl[2]);
675 /* equinox */
676 swi_cross_prod(peqr, pecl, v);
677 w = sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
678 eqx[0] = v[0] / w;
679 eqx[1] = v[1] / w;
680 eqx[2] = v[2] / w;
681 swi_cross_prod(peqr, eqx, v);
682 rp[0] = eqx[0];
683 rp[1] = eqx[1];
684 rp[2] = eqx[2];
685 rp[3] = v[0];
686 rp[4] = v[1];
687 rp[5] = v[2];
688 rp[6] = peqr[0];
689 rp[7] = peqr[1];
690 rp[8] = peqr[2];
691 // int i;
692 // for (i = 0; i < 3; i++) {
693 // fprintf(stderr, "(%.17f %.17f %.17f)\n", rp[i*3], rp[i*3+1],rp[i*3+2]);
694 // } /**/
695 }
696
697 /* precession according to Owen 1990:
698 * Owen, William M., Jr., (JPL) "A Theory of the Earth's Precession
699 * Relative to the Invariable Plane of the Solar System", Ph.D.
700 * Dissertation, University of Florida, 1990.
701 * Implemented for time range -18000 to 14000.
702 */
703 /*
704 * p. 177: central time Tc = -160, covering time span -200 <= T <= -120
705 * i.e. -14000 +- 40 centuries
706 * p. 178: central time Tc = -80, covering time span -120 <= T <= -40
707 * i.e. -6000 +- 40 centuries
708 * p. 179: central time Tc = 0, covering time span -40 <= T <= +40
709 * i.e. 2000 +- 40 centuries
710 * p. 180: central time Tc = 80, covering time span 40 <= T <= 120
711 * i.e. 10000 +- 40 centuries
712 * p. 181: central time Tc = 160, covering time span 120 <= T <= 200
713 * i.e. 10000 +- 40 centuries
714 */
715 static const double owen_eps0_coef[5][10] = {
716 {23.699391439256386, 5.2330816033981775e-1, -5.6259493384864815e-2, -8.2033318431602032e-3, 6.6774163554156385e-4, 2.4931584012812606e-5, -3.1313623302407878e-6, 2.0343814827951515e-7, 2.9182026615852936e-8, -4.1118760893281951e-9,},
717 {24.124759551704588, -1.2094875596566286e-1, -8.3914869653015218e-2, 3.5357075322387405e-3, 6.4557467824807032e-4, -2.5092064378707704e-5, -1.7631607274450848e-6, 1.3363622791424094e-7, 1.5577817511054047e-8, -2.4613907093017122e-9,},
718 {23.439103144206208, -4.9386077073143590e-1, -2.3965445283267805e-4, 8.6637485629656489e-3, -5.2828151901367600e-5, -4.3951004595359217e-5, -1.1058785949914705e-6, 6.2431490022621172e-8, 3.4725376218710764e-8, 1.3658853127005757e-9,},
719 {22.724671295125046, -1.6041813558650337e-1, 7.0646783888132504e-2, 1.4967806745062837e-3, -6.6857270989190734e-4, 5.7578378071604775e-6, 3.3738508454638728e-6, -2.2917813537654764e-7, -2.1019907929218137e-8, 4.3139832091694682e-9,},
720 {22.914636050333696, 3.2123508304962416e-1, 3.6633220173792710e-2, -5.9228324767696043e-3, -1.882379107379328e-4, 3.2274552870236244e-5, 4.9052463646336507e-7, -5.9064298731578425e-8, -2.0485712675098837e-8, -6.2163304813908160e-10,},
721 };
722
723 static const double owen_psia_coef[5][10] = {
724 {-218.57864954903122, 51.752257487741612, 1.3304715765661958e-1, 9.2048123521890745e-2, -6.0877528127241278e-3, -7.0013893644531700e-5, -4.9217728385458495e-5, -1.8578234189053723e-6, 7.4396426162029877e-7, -5.9157528981843864e-9,},
725 {-111.94350527506128, 55.175558131675861, 4.7366115762797613e-1, -4.7701750975398538e-2, -9.2445765329325809e-3, 7.0962838707454917e-4, 1.5140455277814658e-4, -7.7813159018954928e-7, -2.4729402281953378e-6, -1.0898887008726418e-7,},
726 {-2.041452011529441e-1, 55.969995858494106, -1.9295093699770936e-1, -5.6819574830421158e-3, 1.1073687302518981e-2, -9.0868489896815619e-5, -1.1999773777895820e-4, 9.9748697306154409e-6, 5.7911493603430550e-7, -2.3647526839778175e-7,},
727 {111.61366860604471, 56.404525305162447, 4.4403302410703782e-1, 7.1490030578883907e-2, -4.9184559079790816e-3, -1.3912698949042046e-3, -6.8490613661884005e-5, 1.2394328562905297e-6, 1.7719847841480384e-6, 2.4889095220628068e-7,},
728 {228.40683531269390, 60.056143904919826, 2.9583200718478960e-2, -1.5710838319490748e-1, -7.0017356811600801e-3, 3.3009615142224537e-3, 2.0318123852537664e-4, -6.5840216067828310e-5, -5.9077673352976155e-6, 1.3983942185303064e-6,},
729 };
730
731 static const double owen_oma_coef[5][10] = {
732 {25.541291140949806, 2.377889511272162e-1, -3.7337334723142133e-1, 2.4579295485161534e-2, 4.3840999514263623e-3, -3.1126873333599556e-4, -9.8443045771748915e-6, -7.9403103080496923e-7, 1.0840116743893556e-9, 9.2865105216887919e-9,},
733 {24.429357654237926, -9.5205745947740161e-1, 8.6738296270534816e-2, 3.0061543426062955e-2, -4.1532480523019988e-3, -3.7920928393860939e-4, 3.5117012399609737e-5, 4.6811877283079217e-6, -8.1836046585546861e-8, -6.1803706664211173e-8,},
734 {23.450465062489337, -9.7259278279739817e-2, 1.1082286925130981e-2, -3.1469883339372219e-2, -1.0041906996819648e-4, 5.6455168475133958e-4, -8.4403910211030209e-6, -3.8269157371098435e-6, 3.1422585261198437e-7, 9.3481729116773404e-9,},
735 {22.581778052947806, -8.7069701538602037e-1, -9.8140710050197307e-2, 2.6025931340678079e-2, 4.8165322168786755e-3, -1.906558772193363e-4, -4.6838759635421777e-5, -1.6608525315998471e-6, -3.2347811293516124e-8, 2.8104728109642000e-9,},
736 {21.518861835737142, 2.0494789509441385e-1, 3.5193604846503161e-1, 1.5305977982348925e-2, -7.5015367726336455e-3, -4.0322553186065610e-4, 1.0655320434844041e-4, 7.1792339586935752e-6, -1.603874697543020e-6, -1.613563462813512e-7,},
737 };
738
739 static const double owen_chia_coef[5][10] = {
740 {8.2378850337329404e-1, -3.7443109739678667, 4.0143936898854026e-1, 8.1822830214590811e-2, -8.5978790792656293e-3, -2.8350488448426132e-5, -4.2474671728156727e-5, -1.6214840884656678e-6, 7.8560442001953050e-7, -1.032016641696707e-8,},
741 {-2.1726062070318606, 7.8470515033132925e-1, 4.4044931004195718e-1, -8.0671247169971653e-2, -8.9672662444325007e-3, 9.2248978383109719e-4, 1.5143472266372874e-4, -1.6387009056475679e-6, -2.4405558979328144e-6, -1.0148113464009015e-7,},
742 {-4.8518673570735556e-1, 1.0016737299946743e-1, -4.7074888613099918e-1, -5.8604054305076092e-3, 1.4300208240553435e-2, -6.7127991650300028e-5, -1.3703764889645475e-4, 9.0505213684444634e-6, 6.0368690647808607e-7, -2.2135404747652171e-7,},
743 {-2.0950740076326087, -9.4447359463206877e-1, 4.0940512860493755e-1, 1.0261699700263508e-1, -5.3133241571955160e-3, -1.6634631550720911e-3, -5.9477519536647907e-5, 2.9651387319208926e-6, 1.6434499452070584e-6, 2.3720647656961084e-7,},
744 {6.3315163285678715e-1, 3.5241082918420464, 2.1223076605364606e-1, -1.5648122502767368e-1, -9.1964075390801980e-3, 3.3896161239812411e-3, 2.1485178626085787e-4, -6.6261759864793735e-5, -5.9257969712852667e-6, 1.3918759086160525e-6,},
745 };
746
get_owen_t0_icof(double tjd,double * t0,int * icof)747 static void get_owen_t0_icof(double tjd, double *t0, int *icof)
748 {
749 int i, j = 0;
750 double t0s[5] = {-3392455.5, -470455.5, 2451544.5, 5373544.5, 8295544.5, };
751 *t0 = t0s[0];
752 for (i = 1; i < 5; i++) {
753 if (tjd < (t0s[i-1] + t0s[i]) / 2) {
754 ;
755 } else {
756 *t0 = t0s[i];
757 j++;
758 }
759 }
760 *icof = j;
761 }
762
763 /* precession matrix Owen 1990 */
owen_pre_matrix(double tjd,double * rp,int iflag)764 static void owen_pre_matrix(double tjd, double *rp, int iflag)
765 {
766 int i, icof = 0;
767 double eps0 = 0, chia = 0, psia = 0, oma = 0;
768 double coseps0, sineps0, coschia, sinchia, cospsia, sinpsia, cosoma, sinoma;
769 double k[10], tau[10];
770 double t0;
771 get_owen_t0_icof(tjd, &t0, &icof);
772 // fprintf(stderr, "%d, %.17f\n", icof, t0);
773 tau[0] = 0;
774 tau[1] = (tjd - t0) / 36525.0 / 40.0;
775 for (i = 2; i <= 9; i++) {
776 tau[i] = tau[1] * tau[i-1];
777 }
778 k[0] = 1;
779 k[1] = tau[1];
780 k[2] = 2 * tau[2] - 1;
781 k[3] = 4 * tau[3] - 3 * tau[1];
782 k[4] = 8 * tau[4] - 8 * tau[2] + 1;
783 k[5] = 16 * tau[5] - 20 * tau[3] + 5 * tau[1];
784 k[6] = 32 * tau[6] - 48 * tau[4] + 18 * tau[2] - 1;
785 k[7] = 64 * tau[7] - 112 * tau[5] + 56 * tau[3] - 7 * tau[1];
786 k[8] = 128 * tau[8] - 256 * tau[6] + 160 * tau[4] - 32 * tau[2] + 1;
787 k[9] = 256 * tau[9] - 576 * tau[7] + 432 * tau[5] - 120 * tau[3] + 9 * tau[1];
788 for (i = 0; i < 10; i++) {
789 //eps += (k[i] * owen_eps0_coef[icof][i]);
790 psia += (k[i] * owen_psia_coef[icof][i]);
791 oma += (k[i] * owen_oma_coef[icof][i]);
792 chia += (k[i] * owen_chia_coef[icof][i]);
793 }
794 if (iflag & (SEFLG_JPLHOR | SEFLG_JPLHOR_APPROX)) {
795 /*
796 * In comparison with JPL Horizons we have an almost constant offset
797 * almost constant offset in ecl. longitude of about -0.000019 deg.
798 * We fix this as follows: */
799 psia += -0.000018560;
800 }
801 eps0 = 84381.448 / 3600.0;
802 //fprintf(stderr, "e=%.17f, ps=%.17f, om=%.17f, ch=%.17f\n", eps0, psia, oma, chia);
803 eps0 *= DEGTORAD;
804 psia *= DEGTORAD;
805 chia *= DEGTORAD;
806 oma *= DEGTORAD;
807 coseps0 = cos(eps0);
808 sineps0 = sin(eps0);
809 coschia = cos(chia);
810 sinchia = sin(chia);
811 cospsia = cos(psia);
812 sinpsia = sin(psia);
813 cosoma = cos(oma);
814 sinoma = sin(oma);
815 rp[0] = coschia * cospsia + sinchia * cosoma * sinpsia;
816 rp[1] = (-coschia * sinpsia + sinchia * cosoma * cospsia) * coseps0 + sinchia * sinoma * sineps0;
817 rp[2] = (-coschia * sinpsia + sinchia * cosoma * cospsia) * sineps0 - sinchia * sinoma * coseps0;
818 rp[3] = -sinchia * cospsia + coschia * cosoma * sinpsia;
819 rp[4] = (sinchia * sinpsia + coschia * cosoma * cospsia) * coseps0 + coschia * sinoma * sineps0;
820 rp[5] = (sinchia * sinpsia + coschia * cosoma * cospsia) * sineps0 - coschia * sinoma * coseps0;
821 rp[6] = sinoma * sinpsia;
822 rp[7] = sinoma * cospsia * coseps0 - cosoma * sineps0;
823 rp[8] = sinoma * cospsia * sineps0 + cosoma * coseps0;
824 /*for (i = 0; i < 3; i++) {
825 fprintf(stderr, "(%.17f %.17f %.17f)\n", rp[i*3], rp[i*3+1],rp[i*3+2]);
826 } */
827 }
epsiln_owen_1986(double tjd,double * eps)828 static void epsiln_owen_1986(double tjd, double *eps)
829 {
830 int i, icof = 0;
831 double k[10], tau[10];
832 double t0;
833 get_owen_t0_icof(tjd, &t0, &icof);
834 *eps = 0;
835 tau[0] = 0;
836 tau[1] = (tjd - t0) / 36525.0 / 40.0;
837 for (i = 2; i <= 9; i++) {
838 tau[i] = tau[1] * tau[i-1];
839 }
840 k[0] = 1;
841 k[1] = tau[1];
842 k[2] = 2 * tau[2] - 1;
843 k[3] = 4 * tau[3] - 3 * tau[1];
844 k[4] = 8 * tau[4] - 8 * tau[2] + 1;
845 k[5] = 16 * tau[5] - 20 * tau[3] + 5 * tau[1];
846 k[6] = 32 * tau[6] - 48 * tau[4] + 18 * tau[2] - 1;
847 k[7] = 64 * tau[7] - 112 * tau[5] + 56 * tau[3] - 7 * tau[1];
848 k[8] = 128 * tau[8] - 256 * tau[6] + 160 * tau[4] - 32 * tau[2] + 1;
849 k[9] = 256 * tau[9] - 576 * tau[7] + 432 * tau[5] - 120 * tau[3] + 9 * tau[1];
850 for (i = 0; i < 10; i++) {
851 *eps += (k[i] * owen_eps0_coef[icof][i]);
852 }
853 //fprintf(stderr, "eps=%.17f\n", *eps);
854 }
855
856 /* Obliquity of the ecliptic at Julian date J
857 *
858 * IAU Coefficients are from:
859 * J. H. Lieske, T. Lederle, W. Fricke, and B. Morando,
860 * "Expressions for the Precession Quantities Based upon the IAU
861 * (1976) System of Astronomical Constants," Astronomy and Astrophysics
862 * 58, 1-16 (1977).
863 *
864 * Before or after 200 years from J2000, the formula used is from:
865 * J. Laskar, "Secular terms of classical planetary theories
866 * using the results of general theory," Astronomy and Astrophysics
867 * 157, 59070 (1986).
868 *
869 * Bretagnon, P. et al.: 2003, "Expressions for Precession Consistent with
870 * the IAU 2000A Model". A&A 400,785
871 *B03 84381.4088 -46.836051*t -1667*10-7*t2 +199911*10-8*t3 -523*10-9*t4 -248*10-10*t5 -3*10-11*t6
872 *C03 84381.406 -46.836769*t -1831*10-7*t2 +20034*10-7*t3 -576*10-9*t4 -434*10-10*t5
873 *
874 * See precess and page B18 of the Astronomical Almanac.
875 */
876 #define OFFSET_EPS_JPLHORIZONS (35.95)
877 #define DCOR_EPS_JPL_TJD0 2437846.5
878 #define NDCOR_EPS_JPL 51
879 double dcor_eps_jpl[] = {
880 36.726, 36.627, 36.595, 36.578, 36.640, 36.659, 36.731, 36.765,
881 36.662, 36.555, 36.335, 36.321, 36.354, 36.227, 36.289, 36.348, 36.257, 36.163,
882 35.979, 35.896, 35.842, 35.825, 35.912, 35.950, 36.093, 36.191, 36.009, 35.943,
883 35.875, 35.771, 35.788, 35.753, 35.822, 35.866, 35.771, 35.732, 35.543, 35.498,
884 35.449, 35.409, 35.497, 35.556, 35.672, 35.760, 35.596, 35.565, 35.510, 35.394,
885 35.385, 35.375, 35.415,
886 };
swi_epsiln(double J,int32 iflag)887 double swi_epsiln(double J, int32 iflag)
888 {
889 double T, eps;
890 double tofs, dofs, t0, t1;
891 int prec_model = swed.astro_models[SE_MODEL_PREC_LONGTERM];
892 int prec_model_short = swed.astro_models[SE_MODEL_PREC_SHORTTERM];
893 int jplhora_model = swed.astro_models[SE_MODEL_JPLHORA_MODE];
894 AS_BOOL is_jplhor = FALSE;
895 if (prec_model == 0) prec_model = SEMOD_PREC_DEFAULT;
896 if (prec_model_short == 0) prec_model_short = SEMOD_PREC_DEFAULT_SHORT;
897 if (jplhora_model == 0) jplhora_model = SEMOD_JPLHORA_DEFAULT;
898 if (iflag & SEFLG_JPLHOR)
899 is_jplhor = TRUE;
900 if ((iflag & SEFLG_JPLHOR_APPROX)
901 && jplhora_model == SEMOD_JPLHORA_3
902 && J <= HORIZONS_TJD0_DPSI_DEPS_IAU1980)
903 is_jplhor = TRUE;
904 T = (J - 2451545.0)/36525.0;
905 if (is_jplhor) {
906 if (J > 2378131.5 && J < 2525323.5) { // between 1.1.1799 and 1.1.2202
907 eps = (((1.813e-3*T-5.9e-4)*T-46.8150)*T+84381.448)*DEGTORAD/3600;
908 } else {
909 epsiln_owen_1986(J, &eps);
910 eps *= DEGTORAD;
911 }
912 } else if ((iflag & SEFLG_JPLHOR_APPROX) && jplhora_model == SEMOD_JPLHORA_2) {
913 eps = (((1.813e-3*T-5.9e-4)*T-46.8150)*T+84381.448)*DEGTORAD/3600;
914 } else if (prec_model_short == SEMOD_PREC_IAU_1976 && fabs(T) <= PREC_IAU_1976_CTIES ) {
915 eps = (((1.813e-3*T-5.9e-4)*T-46.8150)*T+84381.448)*DEGTORAD/3600;
916 } else if (prec_model == SEMOD_PREC_IAU_1976) {
917 eps = (((1.813e-3*T-5.9e-4)*T-46.8150)*T+84381.448)*DEGTORAD/3600;
918 } else if (prec_model_short == SEMOD_PREC_IAU_2000 && fabs(T) <= PREC_IAU_2000_CTIES ) {
919 eps = (((1.813e-3*T-5.9e-4)*T-46.84024)*T+84381.406)*DEGTORAD/3600;
920 } else if (prec_model == SEMOD_PREC_IAU_2000) {
921 eps = (((1.813e-3*T-5.9e-4)*T-46.84024)*T+84381.406)*DEGTORAD/3600;
922 } else if (prec_model_short == SEMOD_PREC_IAU_2006 && fabs(T) <= PREC_IAU_2006_CTIES) {
923 eps = (((((-4.34e-8 * T -5.76e-7) * T +2.0034e-3) * T -1.831e-4) * T -46.836769) * T + 84381.406) * DEGTORAD / 3600.0;
924 } else if (prec_model == SEMOD_PREC_NEWCOMB) {
925 double Tn = (J - 2396758.0)/36525.0;
926 eps = (0.0017 * Tn * Tn * Tn - 0.0085 * Tn * Tn - 46.837 * Tn + 84451.68) * DEGTORAD / 3600.0;
927 } else if (prec_model == SEMOD_PREC_IAU_2006) {
928 eps = (((((-4.34e-8 * T -5.76e-7) * T +2.0034e-3) * T -1.831e-4) * T -46.836769) * T + 84381.406) * DEGTORAD / 3600.0;
929 } else if (prec_model == SEMOD_PREC_BRETAGNON_2003) {
930 eps = ((((((-3e-11 * T - 2.48e-8) * T -5.23e-7) * T +1.99911e-3) * T -1.667e-4) * T -46.836051) * T + 84381.40880) * DEGTORAD / 3600.0;/* */
931 } else if (prec_model == SEMOD_PREC_SIMON_1994) {
932 eps = (((((2.5e-8 * T -5.1e-7) * T +1.9989e-3) * T -1.52e-4) * T -46.80927) * T + 84381.412) * DEGTORAD / 3600.0;/* */
933 } else if (prec_model == SEMOD_PREC_WILLIAMS_1994) {
934 eps = ((((-1.0e-6 * T +2.0e-3) * T -1.74e-4) * T -46.833960) * T + 84381.409) * DEGTORAD / 3600.0;/* */
935 } else if (prec_model == SEMOD_PREC_LASKAR_1986 || prec_model == SEMOD_PREC_WILL_EPS_LASK) {
936 T /= 10.0;
937 eps = ((((((((( 2.45e-10*T + 5.79e-9)*T + 2.787e-7)*T
938 + 7.12e-7)*T - 3.905e-5)*T - 2.4967e-3)*T
939 - 5.138e-3)*T + 1.99925)*T - 0.0155)*T - 468.093)*T
940 + 84381.448;
941 eps *= DEGTORAD/3600.0;
942 } else if (prec_model == SEMOD_PREC_OWEN_1990) {
943 epsiln_owen_1986(J, &eps);
944 eps *= DEGTORAD;
945 //fprintf(stderr, "epso=%.17f\n", eps);
946 } else { /* SEMOD_PREC_VONDRAK_2011 */
947 swi_ldp_peps(J, NULL, &eps);
948 if ((iflag & SEFLG_JPLHOR_APPROX) && jplhora_model != SEMOD_JPLHORA_2) {
949 tofs = (J - DCOR_EPS_JPL_TJD0) / 365.25;
950 dofs = OFFSET_EPS_JPLHORIZONS;
951 if (tofs < 0) {
952 tofs = 0;
953 dofs = dcor_eps_jpl[0];
954 } else if (tofs >= NDCOR_EPS_JPL - 1) {
955 tofs = NDCOR_EPS_JPL;
956 dofs = dcor_eps_jpl[NDCOR_EPS_JPL - 1];
957 } else {
958 t0 = (int) tofs;
959 t1 = t0 + 1;
960 dofs = dcor_eps_jpl[(int)t0];
961 dofs = (tofs - t0) * (dcor_eps_jpl[(int)t0] - dcor_eps_jpl[(int)t1]) + dcor_eps_jpl[(int)t0];
962 }
963 dofs /= (1000.0 * 3600.0);
964 eps += dofs * DEGTORAD;
965 }
966 //fprintf(stderr, "epsv=%.17f\n", eps);
967 }
968 return(eps);
969 }
970
971 /* Precession of the equinox and ecliptic
972 * from epoch Julian date J to or from J2000.0
973 *
974 * Original program by Steve Moshier.
975 * Changes in program structure and implementation of IAU 2003 (P03) and
976 * Vondrak 2011 by Dieter Koch.
977 *
978 * SEMOD_PREC_VONDRAK_2011
979 * J. Vondrák, N. Capitaine, and P. Wallace, "New precession expressions,
980 * valid for long time intervals", A&A 534, A22 (2011)
981 *
982 * SEMOD_PREC_IAU_2006
983 * N. Capitaine, P.T. Wallace, and J. Chapront, "Expressions for IAU 2000
984 * precession quantities", 2003, A&A 412, 567-586 (2003).
985 * This is a "short" term model, that can be combined with other models
986 *
987 * SEMOD_PREC_WILLIAMS_1994
988 * James G. Williams, "Contributions to the Earth's obliquity rate,
989 * precession, and nutation," Astron. J. 108, 711-724 (1994).
990 *
991 * SEMOD_PREC_SIMON_1994
992 * J. L. Simon, P. Bretagnon, J. Chapront, M. Chapront-Touze', G. Francou,
993 * and J. Laskar, "Numerical Expressions for precession formulae and
994 * mean elements for the Moon and the planets," Astronomy and Astrophysics
995 * 282, 663-683 (1994).
996 *
997 * SEMOD_PREC_IAU_1976
998 * IAU Coefficients are from:
999 * J. H. Lieske, T. Lederle, W. Fricke, and B. Morando,
1000 * "Expressions for the Precession Quantities Based upon the IAU
1001 * (1976) System of Astronomical Constants," Astronomy and
1002 * Astrophysics 58, 1-16 (1977).
1003 * This is a "short" term model, that can be combined with other models
1004 *
1005 * SEMOD_PREC_LASKAR_1986
1006 * Newer formulas that cover a much longer time span are from:
1007 * J. Laskar, "Secular terms of classical planetary theories
1008 * using the results of general theory," Astronomy and Astrophysics
1009 * 157, 59070 (1986).
1010 *
1011 * See also:
1012 * P. Bretagnon and G. Francou, "Planetary theories in rectangular
1013 * and spherical variables. VSOP87 solutions," Astronomy and
1014 * Astrophysics 202, 309-315 (1988).
1015 *
1016 * Bretagnon and Francou's expansions for the node and inclination
1017 * of the ecliptic were derived from Laskar's data but were truncated
1018 * after the term in T**6. I have recomputed these expansions from
1019 * Laskar's data, retaining powers up to T**10 in the result.
1020 *
1021 */
1022
precess_1(double * R,double J,int direction,int prec_method)1023 static int precess_1(double *R, double J, int direction, int prec_method)
1024 {
1025 double T, Z = 0, z = 0, TH = 0;
1026 int i;
1027 double x[3];
1028 double sinth, costh, sinZ, cosZ, sinz, cosz, A, B;
1029 if( J == J2000 )
1030 return(0);
1031 T = (J - J2000)/36525.0;
1032 if (prec_method == SEMOD_PREC_IAU_1976) {
1033 Z = (( 0.017998*T + 0.30188)*T + 2306.2181)*T*DEGTORAD/3600;
1034 z = (( 0.018203*T + 1.09468)*T + 2306.2181)*T*DEGTORAD/3600;
1035 TH = ((-0.041833*T - 0.42665)*T + 2004.3109)*T*DEGTORAD/3600;
1036 /*
1037 * precession relative to ecliptic of start epoch is:
1038 * pn = (5029.0966 + 2.22226*T-0.000042*T*T) * t + (1.11161-0.000127*T) * t * t - 0.000113*t*t*t;
1039 * with: t = (tstart - tdate) / 36525.0
1040 * T = (tstart - J2000) / 36525.0
1041 */
1042 } else if (prec_method == SEMOD_PREC_IAU_2000) {
1043 /* AA 2006 B28:*/
1044 Z = (((((- 0.0000002*T - 0.0000327)*T + 0.0179663)*T + 0.3019015)*T + 2306.0809506)*T + 2.5976176)*DEGTORAD/3600;
1045 z = (((((- 0.0000003*T - 0.000047)*T + 0.0182237)*T + 1.0947790)*T + 2306.0803226)*T - 2.5976176)*DEGTORAD/3600;
1046 TH = ((((-0.0000001*T - 0.0000601)*T - 0.0418251)*T - 0.4269353)*T + 2004.1917476)*T*DEGTORAD/3600;
1047 } else if (prec_method == SEMOD_PREC_IAU_2006) {
1048 T = (J - J2000)/36525.0;
1049 Z = (((((- 0.0000003173*T - 0.000005971)*T + 0.01801828)*T + 0.2988499)*T + 2306.083227)*T + 2.650545)*DEGTORAD/3600;
1050 z = (((((- 0.0000002904*T - 0.000028596)*T + 0.01826837)*T + 1.0927348)*T + 2306.077181)*T - 2.650545)*DEGTORAD/3600;
1051 TH = ((((-0.00000011274*T - 0.000007089)*T - 0.04182264)*T - 0.4294934)*T + 2004.191903)*T*DEGTORAD/3600;
1052 } else if (prec_method == SEMOD_PREC_BRETAGNON_2003) {
1053 Z = ((((((-0.00000000013*T - 0.0000003040)*T - 0.000005708)*T + 0.01801752)*T + 0.3023262)*T + 2306.080472)*T + 2.72767)*DEGTORAD/3600;
1054 z = ((((((-0.00000000005*T - 0.0000002486)*T - 0.000028276)*T + 0.01826676)*T + 1.0956768)*T + 2306.076070)*T - 2.72767)*DEGTORAD/3600;
1055 TH = ((((((0.000000000009*T + 0.00000000036)*T -0.0000001127)*T - 0.000007291)*T - 0.04182364)*T - 0.4266980)*T + 2004.190936)*T*DEGTORAD/3600;
1056 #if 0
1057 } else if (prec_method == SEMOD_PREC_NEWCOMB) {
1058 double t1 = (J2000 - 2415020.3135) / 36524.2199;
1059 double T = (J - J2000) / 36524.2199;
1060 double T2 = T * T; double T3 = T2 * T;
1061 Z = (2304.250 + 1.396 * t1) * T + 0.302 * T2 + 0.0179 * T3;
1062 z = (2304.250 + 1.396 * t1) * T + 1.093 * T2 + 0.0192 * T3;
1063 TH =(2004.682 - 0.853 * t1) * T - 0.426 * T2 - 0.0416 * T3;
1064 Z *= (DEGTORAD/3600.0);
1065 z *= (DEGTORAD/3600.0);
1066 TH *= (DEGTORAD/3600.0);
1067 #endif
1068 #if 0
1069 // from Newcomb, "Compendium" (1906), pp. 245f., relative to 1850
1070 /* } else if (prec_method == SEMOD_PREC_NEWCOMB) {
1071 double cties = 36524.2198782; // trop. centuries
1072 double T = (J - B1850) / cties;
1073 double T2 = T * T; double T3 = T2 * T;
1074 double Z1 = 2303.56;
1075 Z = 2303.56 * T + 0.3023 * T2 + 0.018 * T3;
1076 z = 2303.55 * T + 1.094 * T2 + 0.018 * T3;
1077 TH = 2005.11 * T - 0.43 * T2 - 0.041 * T3;
1078 Z *= (DEGTORAD/3600.0);
1079 z *= (DEGTORAD/3600.0);
1080 TH *= (DEGTORAD/3600.0);
1081 */
1082 #endif
1083 #if 0
1084 // Newcomb from Expl. supp. 61 pg. 38
1085 // "Andoyar (Woolard and Clemence) expressions":
1086 } else if (prec_method == SEMOD_PREC_NEWCOMB) {
1087 double mills = 365242.198782; // trop. millennia
1088 double t1 = (J2000 - B1850) / mills;
1089 double t2 = (J - B1850) / mills;
1090 double T = t2 - t1;
1091 double T2 = T * T; double T3 = T2 * T;
1092 double Z1 = 23035.545 + 139.720 * t1 + 0.060 * t1 * t1;
1093 Z = Z1 * T + (30.240 - 0.270 * t1) * T2 + 17.995 * T3;
1094 z = Z1 * T + (109.480 - 0.390 * t1) * T2 + 18.325 * T3;
1095 TH = (20051.12 - 85.29 * t1 - 0.37 * t1 * t1) * T + (-42.65 - 0.37 * t1) * T2 - 41.80 * T3;
1096 Z *= (DEGTORAD/3600.0);
1097 z *= (DEGTORAD/3600.0);
1098 TH *= (DEGTORAD/3600.0);
1099 #endif
1100 #if 1
1101 // Newcomb according to Kinoshita 1975, very close to ExplSuppl/Andoyer;
1102 // one additional digit.
1103 } else if (prec_method == SEMOD_PREC_NEWCOMB) {
1104 double mills = 365242.198782; // trop. millennia
1105 double t1 = (J2000 - B1850) / mills;
1106 double t2 = (J - B1850) / mills;
1107 double T = t2 - t1;
1108 double T2 = T * T; double T3 = T2 * T;
1109 double Z1 = 23035.5548 + 139.720 * t1 + 0.069 * t1 * t1;
1110 Z = Z1 * T + (30.242 - 0.269 * t1) * T2 + 17.996 * T3;
1111 z = Z1 * T + (109.478 - 0.387 * t1) * T2 + 18.324 * T3;
1112 TH = (20051.125 - 85.294 * t1 - 0.365 * t1 * t1) * T + (-42.647 - 0.365 * t1) * T2 - 41.802 * T3;
1113 Z *= (DEGTORAD/3600.0);
1114 z *= (DEGTORAD/3600.0);
1115 TH *= (DEGTORAD/3600.0);
1116 #endif
1117 #if 0
1118 // from Lieske, "Expressions for the Precession Quantities..." (1967), p. 20
1119 } else if (prec_method == SEMOD_PREC_NEWCOMB) {
1120 double cties = 36524.2198782; // trop. centuries
1121 double t1 = (J2000 - J1900) / cties;
1122 double t2 = (J - J1900) / cties;
1123 double T = t2 - t1;
1124 double T2 = T * T; double T3 = T2 * T;
1125 double Z1 = 2304.253 + 1.3972 * t1 + 0.000125 * t1 * t1;
1126 Z = Z1 * T + (0.3023 - 0.000211 * t1) * T2 + 0.0180 * T3;
1127 z = Z1 * T + (1.0949 - 0.00046 * t1) * T2 + 0.0183 * T3;
1128 TH = (2004.684 - 0.8532 * t1 - 0.000317 * t1 * t1) * T + (-0.4266 - 0.00032 * t1) * T2 - 0.0418 * T3;
1129 Z *= (DEGTORAD/3600.0);
1130 z *= (DEGTORAD/3600.0);
1131 TH *= (DEGTORAD/3600.0);
1132 #endif
1133 } else {
1134 return 0;
1135 }
1136 sinth = sin(TH);
1137 costh = cos(TH);
1138 sinZ = sin(Z);
1139 cosZ = cos(Z);
1140 sinz = sin(z);
1141 cosz = cos(z);
1142 A = cosZ*costh;
1143 B = sinZ*costh;
1144 if( direction < 0 ) { /* From J2000.0 to J */
1145 x[0] = (A*cosz - sinZ*sinz)*R[0]
1146 - (B*cosz + cosZ*sinz)*R[1]
1147 - sinth*cosz*R[2];
1148 x[1] = (A*sinz + sinZ*cosz)*R[0]
1149 - (B*sinz - cosZ*cosz)*R[1]
1150 - sinth*sinz*R[2];
1151 x[2] = cosZ*sinth*R[0]
1152 - sinZ*sinth*R[1]
1153 + costh*R[2];
1154 } else { /* From J to J2000.0 */
1155 x[0] = (A*cosz - sinZ*sinz)*R[0]
1156 + (A*sinz + sinZ*cosz)*R[1]
1157 + cosZ*sinth*R[2];
1158 x[1] = - (B*cosz + cosZ*sinz)*R[0]
1159 - (B*sinz - cosZ*cosz)*R[1]
1160 - sinZ*sinth*R[2];
1161 x[2] = - sinth*cosz*R[0]
1162 - sinth*sinz*R[1]
1163 + costh*R[2];
1164 }
1165 for( i=0; i<3; i++ )
1166 R[i] = x[i];
1167 return(0);
1168 }
1169
1170 /* In WILLIAMS and SIMON, Laskar's terms of order higher than t^4
1171 have been retained, because Simon et al mention that the solution
1172 is the same except for the lower order terms. */
1173
1174 /* SEMOD_PREC_WILLIAMS_1994 */
1175 static const double pAcof_williams[] = {
1176 -8.66e-10, -4.759e-8, 2.424e-7, 1.3095e-5, 1.7451e-4, -1.8055e-3,
1177 -0.235316, 0.076, 110.5407, 50287.70000 };
1178 static const double nodecof_williams[] = {
1179 6.6402e-16, -2.69151e-15, -1.547021e-12, 7.521313e-12, 1.9e-10,
1180 -3.54e-9, -1.8103e-7, 1.26e-7, 7.436169e-5,
1181 -0.04207794833, 3.052115282424};
1182 static const double inclcof_williams[] = {
1183 1.2147e-16, 7.3759e-17, -8.26287e-14, 2.503410e-13, 2.4650839e-11,
1184 -5.4000441e-11, 1.32115526e-9, -6.012e-7, -1.62442e-5,
1185 0.00227850649, 0.0 };
1186
1187 /* SEMOD_PREC_SIMON_1994 */
1188 /* Precession coefficients from Simon et al: */
1189 static const double pAcof_simon[] = {
1190 -8.66e-10, -4.759e-8, 2.424e-7, 1.3095e-5, 1.7451e-4, -1.8055e-3,
1191 -0.235316, 0.07732, 111.2022, 50288.200 };
1192 static const double nodecof_simon[] = {
1193 6.6402e-16, -2.69151e-15, -1.547021e-12, 7.521313e-12, 1.9e-10,
1194 -3.54e-9, -1.8103e-7, 2.579e-8, 7.4379679e-5,
1195 -0.0420782900, 3.0521126906};
1196 static const double inclcof_simon[] = {
1197 1.2147e-16, 7.3759e-17, -8.26287e-14, 2.503410e-13, 2.4650839e-11,
1198 -5.4000441e-11, 1.32115526e-9, -5.99908e-7, -1.624383e-5,
1199 0.002278492868, 0.0 };
1200
1201 /* SEMOD_PREC_LASKAR_1986 */
1202 /* Precession coefficients taken from Laskar's paper: */
1203 static const double pAcof_laskar[] = {
1204 -8.66e-10, -4.759e-8, 2.424e-7, 1.3095e-5, 1.7451e-4, -1.8055e-3,
1205 -0.235316, 0.07732, 111.1971, 50290.966 };
1206 /* Node and inclination of the earth's orbit computed from
1207 * Laskar's data as done in Bretagnon and Francou's paper.
1208 * Units are radians.
1209 */
1210 static const double nodecof_laskar[] = {
1211 6.6402e-16, -2.69151e-15, -1.547021e-12, 7.521313e-12, 6.3190131e-10,
1212 -3.48388152e-9, -1.813065896e-7, 2.75036225e-8, 7.4394531426e-5,
1213 -0.042078604317, 3.052112654975 };
1214 static const double inclcof_laskar[] = {
1215 1.2147e-16, 7.3759e-17, -8.26287e-14, 2.503410e-13, 2.4650839e-11,
1216 -5.4000441e-11, 1.32115526e-9, -5.998737027e-7, -1.6242797091e-5,
1217 0.002278495537, 0.0 };
1218
precess_2(double * R,double J,int32 iflag,int direction,int prec_method)1219 static int precess_2(double *R, double J, int32 iflag, int direction, int prec_method)
1220 {
1221 int i;
1222 double T, z;
1223 double eps, sineps, coseps;
1224 double x[3];
1225 const double *p;
1226 double A, B, pA, W;
1227 const double *pAcof, *inclcof, *nodecof;
1228 if( J == J2000 )
1229 return(0);
1230 if (prec_method == SEMOD_PREC_LASKAR_1986) {
1231 pAcof = pAcof_laskar;
1232 nodecof = nodecof_laskar;
1233 inclcof = inclcof_laskar;
1234 } else if (prec_method == SEMOD_PREC_SIMON_1994) {
1235 pAcof = pAcof_simon;
1236 nodecof = nodecof_simon;
1237 inclcof = inclcof_simon;
1238 } else if (prec_method == SEMOD_PREC_WILLIAMS_1994) {
1239 pAcof = pAcof_williams;
1240 nodecof = nodecof_williams;
1241 inclcof = inclcof_williams;
1242 } else { /* default, to satisfy compiler */
1243 pAcof = pAcof_laskar;
1244 nodecof = nodecof_laskar;
1245 inclcof = inclcof_laskar;
1246 }
1247 T = (J - J2000)/36525.0;
1248 /* Implementation by elementary rotations using Laskar's expansions.
1249 * First rotate about the x axis from the initial equator
1250 * to the ecliptic. (The input is equatorial.)
1251 */
1252 if( direction == 1 )
1253 eps = swi_epsiln(J, iflag); /* To J2000 */
1254 else
1255 eps = swi_epsiln(J2000, iflag); /* From J2000 */
1256 sineps = sin(eps);
1257 coseps = cos(eps);
1258 x[0] = R[0];
1259 z = coseps*R[1] + sineps*R[2];
1260 x[2] = -sineps*R[1] + coseps*R[2];
1261 x[1] = z;
1262 /* Precession in longitude */
1263 T /= 10.0; /* thousands of years */
1264 p = pAcof;
1265 pA = *p++;
1266 for( i=0; i<9; i++ ) {
1267 pA = pA * T + *p++;
1268 }
1269 pA *= DEGTORAD/3600 * T;
1270 /* Node of the moving ecliptic on the J2000 ecliptic.
1271 */
1272 p = nodecof;
1273 W = *p++;
1274 for( i=0; i<10; i++ )
1275 W = W * T + *p++;
1276 /* Rotate about z axis to the node.
1277 */
1278 if( direction == 1 )
1279 z = W + pA;
1280 else
1281 z = W;
1282 B = cos(z);
1283 A = sin(z);
1284 z = B * x[0] + A * x[1];
1285 x[1] = -A * x[0] + B * x[1];
1286 x[0] = z;
1287 /* Rotate about new x axis by the inclination of the moving
1288 * ecliptic on the J2000 ecliptic.
1289 */
1290 p = inclcof;
1291 z = *p++;
1292 for( i=0; i<10; i++ )
1293 z = z * T + *p++;
1294 if( direction == 1 )
1295 z = -z;
1296 B = cos(z);
1297 A = sin(z);
1298 z = B * x[1] + A * x[2];
1299 x[2] = -A * x[1] + B * x[2];
1300 x[1] = z;
1301 /* Rotate about new z axis back from the node.
1302 */
1303 if( direction == 1 )
1304 z = -W;
1305 else
1306 z = -W - pA;
1307 B = cos(z);
1308 A = sin(z);
1309 z = B * x[0] + A * x[1];
1310 x[1] = -A * x[0] + B * x[1];
1311 x[0] = z;
1312 /* Rotate about x axis to final equator.
1313 */
1314 if( direction == 1 )
1315 eps = swi_epsiln(J2000, iflag);
1316 else
1317 eps = swi_epsiln(J, iflag);
1318 sineps = sin(eps);
1319 coseps = cos(eps);
1320 z = coseps * x[1] - sineps * x[2];
1321 x[2] = sineps * x[1] + coseps * x[2];
1322 x[1] = z;
1323 for( i=0; i<3; i++ )
1324 R[i] = x[i];
1325 return(0);
1326 }
1327
precess_3(double * R,double J,int direction,int iflag,int prec_meth)1328 static int precess_3(double *R, double J, int direction, int iflag, int prec_meth)
1329 {
1330 //double T;
1331 double x[3], pmat[9];
1332 int i, j;
1333 if( J == J2000 )
1334 return(0);
1335 /* Each precession angle is specified by a polynomial in
1336 * T = Julian centuries from J2000.0. See AA page B18.
1337 */
1338 //T = (J - J2000)/36525.0;
1339 if (prec_meth == SEMOD_PREC_OWEN_1990)
1340 owen_pre_matrix(J, pmat, iflag);
1341 else
1342 pre_pmat(J, pmat);
1343 if (direction == -1) {
1344 for (i = 0, j = 0; i <= 2; i++, j = i * 3) {
1345 x[i] = R[0] * pmat[j + 0] +
1346 R[1] * pmat[j + 1] +
1347 R[2] * pmat[j + 2];
1348 }
1349 } else {
1350 for (i = 0, j = 0; i <= 2; i++, j = i * 3) {
1351 x[i] = R[0] * pmat[i + 0] +
1352 R[1] * pmat[i + 3] +
1353 R[2] * pmat[i + 6];
1354 }
1355 }
1356 for (i = 0; i < 3; i++)
1357 R[i] = x[i];
1358 return(0);
1359 }
1360
1361 /* Subroutine arguments:
1362 *
1363 * R = rectangular equatorial coordinate vector to be precessed.
1364 * The result is written back into the input vector.
1365 * J = Julian date
1366 * direction =
1367 * Precess from J to J2000: direction = 1
1368 * Precess from J2000 to J: direction = -1
1369 * Note that if you want to precess from J1 to J2, you would
1370 * first go from J1 to J2000, then call the program again
1371 * to go from J2000 to J2.
1372 */
swi_precess(double * R,double J,int32 iflag,int direction)1373 int swi_precess(double *R, double J, int32 iflag, int direction )
1374 {
1375 double T = (J - J2000)/36525.0;
1376 int prec_model = swed.astro_models[SE_MODEL_PREC_LONGTERM];
1377 int prec_model_short = swed.astro_models[SE_MODEL_PREC_SHORTTERM];
1378 int jplhora_model = swed.astro_models[SE_MODEL_JPLHORA_MODE];
1379 AS_BOOL is_jplhor = FALSE;
1380 if (prec_model == 0) prec_model = SEMOD_PREC_DEFAULT;
1381 if (prec_model_short == 0) prec_model_short = SEMOD_PREC_DEFAULT_SHORT;
1382 if (jplhora_model == 0) jplhora_model = SEMOD_JPLHORA_DEFAULT;
1383 if (iflag & SEFLG_JPLHOR)
1384 is_jplhor = TRUE;
1385 if ((iflag & SEFLG_JPLHOR_APPROX)
1386 && jplhora_model == SEMOD_JPLHORA_3
1387 && J <= HORIZONS_TJD0_DPSI_DEPS_IAU1980)
1388 is_jplhor = TRUE;
1389 /* JPL Horizons uses precession IAU 1976 and nutation IAU 1980 plus
1390 * some correction to nutation, arriving at extremely high precision */
1391 if (is_jplhor) {
1392 if (J > 2378131.5 && J < 2525323.5) { // between 1.1.1799 and 1.1.2202
1393 return precess_1(R, J, direction, SEMOD_PREC_IAU_1976);
1394 } else {
1395 return precess_3(R, J, direction, iflag, SEMOD_PREC_OWEN_1990);
1396 }
1397 /* Use IAU 1976 formula for a few centuries. */
1398 } else if (prec_model_short == SEMOD_PREC_IAU_1976 && fabs(T) <= PREC_IAU_1976_CTIES) {
1399 return precess_1(R, J, direction, SEMOD_PREC_IAU_1976);
1400 } else if (prec_model == SEMOD_PREC_IAU_1976) {
1401 return precess_1(R, J, direction, SEMOD_PREC_IAU_1976);
1402 /* Use IAU 2000 formula for a few centuries. */
1403 } else if (prec_model_short == SEMOD_PREC_IAU_2000 && fabs(T) <= PREC_IAU_2000_CTIES) {
1404 return precess_1(R, J, direction, SEMOD_PREC_IAU_2000);
1405 } else if (prec_model == SEMOD_PREC_IAU_2000) {
1406 return precess_1(R, J, direction, SEMOD_PREC_IAU_2000);
1407 /* Use IAU 2006 formula for a few centuries. */
1408 } else if (prec_model_short == SEMOD_PREC_IAU_2006 && fabs(T) <= PREC_IAU_2006_CTIES) {
1409 return precess_1(R, J, direction, SEMOD_PREC_IAU_2006);
1410 } else if (prec_model == SEMOD_PREC_IAU_2006) {
1411 return precess_1(R, J, direction, SEMOD_PREC_IAU_2006);
1412 } else if (prec_model == SEMOD_PREC_BRETAGNON_2003) {
1413 return precess_1(R, J, direction, SEMOD_PREC_BRETAGNON_2003);
1414 } else if (prec_model == SEMOD_PREC_NEWCOMB) {
1415 return precess_1(R, J, direction, SEMOD_PREC_NEWCOMB);
1416 } else if (prec_model == SEMOD_PREC_LASKAR_1986) {
1417 return precess_2(R, J, iflag, direction, SEMOD_PREC_LASKAR_1986);
1418 } else if (prec_model == SEMOD_PREC_SIMON_1994) {
1419 return precess_2(R, J, iflag, direction, SEMOD_PREC_SIMON_1994);
1420 } else if (prec_model == SEMOD_PREC_WILLIAMS_1994 || prec_model == SEMOD_PREC_WILL_EPS_LASK) {
1421 return precess_2(R, J, iflag, direction, SEMOD_PREC_WILLIAMS_1994);
1422 } else if (prec_model == SEMOD_PREC_OWEN_1990) {
1423 return precess_3(R, J, direction, iflag, SEMOD_PREC_OWEN_1990);
1424 } else { /* SEMOD_PREC_VONDRAK_2011 */
1425 //int i;
1426 //int retval = precess_3(R, J, direction, iflag, SEMOD_PREC_VONDRAK_2011);
1427 // for (i = 0; i < 3; i++) printf("%.9f, ", R[i]);
1428 return precess_3(R, J, direction, iflag, SEMOD_PREC_VONDRAK_2011);
1429 }
1430 }
1431
1432 /* Nutation in longitude and obliquity
1433 * computed at Julian date J.
1434 *
1435 * References:
1436 * "Summary of 1980 IAU Theory of Nutation (Final Report of the
1437 * IAU Working Group on Nutation)", P. K. Seidelmann et al., in
1438 * Transactions of the IAU Vol. XVIII A, Reports on Astronomy,
1439 * P. A. Wayman, ed.; D. Reidel Pub. Co., 1982.
1440 *
1441 * "Nutation and the Earth's Rotation",
1442 * I.A.U. Symposium No. 78, May, 1977, page 256.
1443 * I.A.U., 1980.
1444 *
1445 * Woolard, E.W., "A redevelopment of the theory of nutation",
1446 * The Astronomical Journal, 58, 1-3 (1953).
1447 *
1448 * This program implements all of the 1980 IAU nutation series.
1449 * Results checked at 100 points against the 1986 AA; all agreed.
1450 *
1451 *
1452 * - S. L. Moshier, November 1987
1453 * October, 1992 - typo fixed in nutation matrix
1454 *
1455 * - D. Koch, November 1995: small changes in structure,
1456 * Corrections to IAU 1980 Series added from Expl. Suppl. p. 116
1457 *
1458 * Each term in the expansion has a trigonometric
1459 * argument given by
1460 * W = i*MM + j*MS + k*FF + l*DD + m*OM
1461 * where the variables are defined below.
1462 * The nutation in longitude is a sum of terms of the
1463 * form (a + bT) * sin(W). The terms for nutation in obliquity
1464 * are of the form (c + dT) * cos(W). The coefficients
1465 * are arranged in the tabulation as follows:
1466 *
1467 * Coefficient:
1468 * i j k l m a b c d
1469 * 0, 0, 0, 0, 1, -171996, -1742, 92025, 89,
1470 * The first line of the table, above, is done separately
1471 * since two of the values do not fit into 16 bit integers.
1472 * The values a and c are arc seconds times 10000. b and d
1473 * are arc seconds per Julian century times 100000. i through m
1474 * are integers. See the program for interpretation of MM, MS,
1475 * etc., which are mean orbital elements of the Sun and Moon.
1476 *
1477 * If terms with coefficient less than X are omitted, the peak
1478 * errors will be:
1479 *
1480 * omit error, omit error,
1481 * a < longitude c < obliquity
1482 * .0005" .0100" .0008" .0094"
1483 * .0046 .0492 .0095 .0481
1484 * .0123 .0880 .0224 .0905
1485 * .0386 .1808 .0895 .1129
1486 */
1487 static const short nt[] = {
1488 /* LS and OC are units of 0.0001"
1489 *LS2 and OC2 are units of 0.00001"
1490 *MM,MS,FF,DD,OM, LS, LS2,OC, OC2 */
1491 0, 0, 0, 0, 2, 2062, 2, -895, 5,
1492 -2, 0, 2, 0, 1, 46, 0, -24, 0,
1493 2, 0,-2, 0, 0, 11, 0, 0, 0,
1494 -2, 0, 2, 0, 2, -3, 0, 1, 0,
1495 1,-1, 0,-1, 0, -3, 0, 0, 0,
1496 0,-2, 2,-2, 1, -2, 0, 1, 0,
1497 2, 0,-2, 0, 1, 1, 0, 0, 0,
1498 0, 0, 2,-2, 2,-13187,-16, 5736,-31,
1499 0, 1, 0, 0, 0, 1426,-34, 54, -1,
1500 0, 1, 2,-2, 2, -517, 12, 224, -6,
1501 0,-1, 2,-2, 2, 217, -5, -95, 3,
1502 0, 0, 2,-2, 1, 129, 1, -70, 0,
1503 2, 0, 0,-2, 0, 48, 0, 1, 0,
1504 0, 0, 2,-2, 0, -22, 0, 0, 0,
1505 0, 2, 0, 0, 0, 17, -1, 0, 0,
1506 0, 1, 0, 0, 1, -15, 0, 9, 0,
1507 0, 2, 2,-2, 2, -16, 1, 7, 0,
1508 0,-1, 0, 0, 1, -12, 0, 6, 0,
1509 -2, 0, 0, 2, 1, -6, 0, 3, 0,
1510 0,-1, 2,-2, 1, -5, 0, 3, 0,
1511 2, 0, 0,-2, 1, 4, 0, -2, 0,
1512 0, 1, 2,-2, 1, 4, 0, -2, 0,
1513 1, 0, 0,-1, 0, -4, 0, 0, 0,
1514 2, 1, 0,-2, 0, 1, 0, 0, 0,
1515 0, 0,-2, 2, 1, 1, 0, 0, 0,
1516 0, 1,-2, 2, 0, -1, 0, 0, 0,
1517 0, 1, 0, 0, 2, 1, 0, 0, 0,
1518 -1, 0, 0, 1, 1, 1, 0, 0, 0,
1519 0, 1, 2,-2, 0, -1, 0, 0, 0,
1520 0, 0, 2, 0, 2, -2274, -2, 977, -5,
1521 1, 0, 0, 0, 0, 712, 1, -7, 0,
1522 0, 0, 2, 0, 1, -386, -4, 200, 0,
1523 1, 0, 2, 0, 2, -301, 0, 129, -1,
1524 1, 0, 0,-2, 0, -158, 0, -1, 0,
1525 -1, 0, 2, 0, 2, 123, 0, -53, 0,
1526 0, 0, 0, 2, 0, 63, 0, -2, 0,
1527 1, 0, 0, 0, 1, 63, 1, -33, 0,
1528 -1, 0, 0, 0, 1, -58, -1, 32, 0,
1529 -1, 0, 2, 2, 2, -59, 0, 26, 0,
1530 1, 0, 2, 0, 1, -51, 0, 27, 0,
1531 0, 0, 2, 2, 2, -38, 0, 16, 0,
1532 2, 0, 0, 0, 0, 29, 0, -1, 0,
1533 1, 0, 2,-2, 2, 29, 0, -12, 0,
1534 2, 0, 2, 0, 2, -31, 0, 13, 0,
1535 0, 0, 2, 0, 0, 26, 0, -1, 0,
1536 -1, 0, 2, 0, 1, 21, 0, -10, 0,
1537 -1, 0, 0, 2, 1, 16, 0, -8, 0,
1538 1, 0, 0,-2, 1, -13, 0, 7, 0,
1539 -1, 0, 2, 2, 1, -10, 0, 5, 0,
1540 1, 1, 0,-2, 0, -7, 0, 0, 0,
1541 0, 1, 2, 0, 2, 7, 0, -3, 0,
1542 0,-1, 2, 0, 2, -7, 0, 3, 0,
1543 1, 0, 2, 2, 2, -8, 0, 3, 0,
1544 1, 0, 0, 2, 0, 6, 0, 0, 0,
1545 2, 0, 2,-2, 2, 6, 0, -3, 0,
1546 0, 0, 0, 2, 1, -6, 0, 3, 0,
1547 0, 0, 2, 2, 1, -7, 0, 3, 0,
1548 1, 0, 2,-2, 1, 6, 0, -3, 0,
1549 0, 0, 0,-2, 1, -5, 0, 3, 0,
1550 1,-1, 0, 0, 0, 5, 0, 0, 0,
1551 2, 0, 2, 0, 1, -5, 0, 3, 0,
1552 0, 1, 0,-2, 0, -4, 0, 0, 0,
1553 1, 0,-2, 0, 0, 4, 0, 0, 0,
1554 0, 0, 0, 1, 0, -4, 0, 0, 0,
1555 1, 1, 0, 0, 0, -3, 0, 0, 0,
1556 1, 0, 2, 0, 0, 3, 0, 0, 0,
1557 1,-1, 2, 0, 2, -3, 0, 1, 0,
1558 -1,-1, 2, 2, 2, -3, 0, 1, 0,
1559 -2, 0, 0, 0, 1, -2, 0, 1, 0,
1560 3, 0, 2, 0, 2, -3, 0, 1, 0,
1561 0,-1, 2, 2, 2, -3, 0, 1, 0,
1562 1, 1, 2, 0, 2, 2, 0, -1, 0,
1563 -1, 0, 2,-2, 1, -2, 0, 1, 0,
1564 2, 0, 0, 0, 1, 2, 0, -1, 0,
1565 1, 0, 0, 0, 2, -2, 0, 1, 0,
1566 3, 0, 0, 0, 0, 2, 0, 0, 0,
1567 0, 0, 2, 1, 2, 2, 0, -1, 0,
1568 -1, 0, 0, 0, 2, 1, 0, -1, 0,
1569
1570 1, 0, 0,-4, 0, -1, 0, 0, 0,
1571 -2, 0, 2, 2, 2, 1, 0, -1, 0,
1572 -1, 0, 2, 4, 2, -2, 0, 1, 0,
1573 2, 0, 0,-4, 0, -1, 0, 0, 0,
1574 1, 1, 2,-2, 2, 1, 0, -1, 0,
1575 1, 0, 2, 2, 1, -1, 0, 1, 0,
1576 -2, 0, 2, 4, 2, -1, 0, 1, 0,
1577 -1, 0, 4, 0, 2, 1, 0, 0, 0,
1578 1,-1, 0,-2, 0, 1, 0, 0, 0,
1579 2, 0, 2,-2, 1, 1, 0, -1, 0,
1580 2, 0, 2, 2, 2, -1, 0, 0, 0,
1581 1, 0, 0, 2, 1, -1, 0, 0, 0,
1582 0, 0, 4,-2, 2, 1, 0, 0, 0,
1583 3, 0, 2,-2, 2, 1, 0, 0, 0,
1584 1, 0, 2,-2, 0, -1, 0, 0, 0,
1585 0, 1, 2, 0, 1, 1, 0, 0, 0,
1586 -1,-1, 0, 2, 1, 1, 0, 0, 0,
1587 0, 0,-2, 0, 1, -1, 0, 0, 0,
1588 0, 0, 2,-1, 2, -1, 0, 0, 0,
1589 0, 1, 0, 2, 0, -1, 0, 0, 0,
1590 1, 0,-2,-2, 0, -1, 0, 0, 0,
1591 0,-1, 2, 0, 1, -1, 0, 0, 0,
1592 1, 1, 0,-2, 1, -1, 0, 0, 0,
1593 1, 0,-2, 2, 0, -1, 0, 0, 0,
1594 2, 0, 0, 2, 0, 1, 0, 0, 0,
1595 0, 0, 2, 4, 2, -1, 0, 0, 0,
1596 0, 1, 0, 1, 0, 1, 0, 0, 0,
1597 #if 1
1598 /*#if NUT_CORR_1987 switch is handled in function calc_nutation_iau1980() */
1599 /* corrections to IAU 1980 nutation series by Herring 1987
1600 * in 0.00001" !!!
1601 * LS OC */
1602 101, 0, 0, 0, 1,-725, 0, 213, 0,
1603 101, 1, 0, 0, 0, 523, 0, 208, 0,
1604 101, 0, 2,-2, 2, 102, 0, -41, 0,
1605 101, 0, 2, 0, 2, -81, 0, 32, 0,
1606 /* LC OS !!! */
1607 102, 0, 0, 0, 1, 417, 0, 224, 0,
1608 102, 1, 0, 0, 0, 61, 0, -24, 0,
1609 102, 0, 2,-2, 2,-118, 0, -47, 0,
1610 /*#endif*/
1611 #endif
1612 ENDMARK,
1613 };
1614
calc_nutation_iau1980(double J,double * nutlo)1615 static int calc_nutation_iau1980(double J, double *nutlo)
1616 {
1617 /* arrays to hold sines and cosines of multiple angles */
1618 double ss[5][8];
1619 double cc[5][8];
1620 double arg;
1621 double args[5];
1622 double f, g, T, T2;
1623 double MM, MS, FF, DD, OM;
1624 double cu, su, cv, sv, sw, s;
1625 double C, D;
1626 int i, j, k, k1, m, n;
1627 int ns[5];
1628 const short *p;
1629 int nut_model = swed.astro_models[SE_MODEL_NUT];
1630 if (nut_model == 0) nut_model = SEMOD_NUT_DEFAULT;
1631 /* Julian centuries from 2000 January 1.5,
1632 * barycentric dynamical time
1633 */
1634 T = (J - 2451545.0) / 36525.0;
1635 T2 = T * T;
1636 /* Fundamental arguments in the FK5 reference system.
1637 * The coefficients, originally given to 0.001",
1638 * are converted here to degrees.
1639 */
1640 /* longitude of the mean ascending node of the lunar orbit
1641 * on the ecliptic, measured from the mean equinox of date
1642 */
1643 OM = -6962890.539 * T + 450160.280 + (0.008 * T + 7.455) * T2;
1644 OM = swe_degnorm(OM/3600) * DEGTORAD;
1645 /* mean longitude of the Sun minus the
1646 * mean longitude of the Sun's perigee
1647 */
1648 MS = 129596581.224 * T + 1287099.804 - (0.012 * T + 0.577) * T2;
1649 MS = swe_degnorm(MS/3600) * DEGTORAD;
1650 /* mean longitude of the Moon minus the
1651 * mean longitude of the Moon's perigee
1652 */
1653 MM = 1717915922.633 * T + 485866.733 + (0.064 * T + 31.310) * T2;
1654 MM = swe_degnorm(MM/3600) * DEGTORAD;
1655 /* mean longitude of the Moon minus the
1656 * mean longitude of the Moon's node
1657 */
1658 FF = 1739527263.137 * T + 335778.877 + (0.011 * T - 13.257) * T2;
1659 FF = swe_degnorm(FF/3600) * DEGTORAD;
1660 /* mean elongation of the Moon from the Sun.
1661 */
1662 DD = 1602961601.328 * T + 1072261.307 + (0.019 * T - 6.891) * T2;
1663 DD = swe_degnorm(DD/3600) * DEGTORAD;
1664 args[0] = MM;
1665 ns[0] = 3;
1666 args[1] = MS;
1667 ns[1] = 2;
1668 args[2] = FF;
1669 ns[2] = 4;
1670 args[3] = DD;
1671 ns[3] = 4;
1672 args[4] = OM;
1673 ns[4] = 2;
1674 /* Calculate sin( i*MM ), etc. for needed multiple angles
1675 */
1676 for (k = 0; k <= 4; k++) {
1677 arg = args[k];
1678 n = ns[k];
1679 su = sin(arg);
1680 cu = cos(arg);
1681 ss[k][0] = su; /* sin(L) */
1682 cc[k][0] = cu; /* cos(L) */
1683 sv = 2.0*su*cu;
1684 cv = cu*cu - su*su;
1685 ss[k][1] = sv; /* sin(2L) */
1686 cc[k][1] = cv;
1687 for( i=2; i<n; i++ ) {
1688 s = su*cv + cu*sv;
1689 cv = cu*cv - su*sv;
1690 sv = s;
1691 ss[k][i] = sv; /* sin( i+1 L ) */
1692 cc[k][i] = cv;
1693 }
1694 }
1695 /* first terms, not in table: */
1696 C = (-0.01742*T - 17.1996)*ss[4][0]; /* sin(OM) */
1697 D = ( 0.00089*T + 9.2025)*cc[4][0]; /* cos(OM) */
1698 for(p = &nt[0]; *p != ENDMARK; p += 9) {
1699 if (nut_model != SEMOD_NUT_IAU_CORR_1987 && (p[0] == 101 || p[0] == 102))
1700 continue;
1701 /* argument of sine and cosine */
1702 k1 = 0;
1703 cv = 0.0;
1704 sv = 0.0;
1705 for( m=0; m<5; m++ ) {
1706 j = p[m];
1707 if (j > 100)
1708 j = 0; /* p[0] is a flag */
1709 if( j ) {
1710 k = j;
1711 if( j < 0 )
1712 k = -k;
1713 su = ss[m][k-1]; /* sin(k*angle) */
1714 if( j < 0 )
1715 su = -su;
1716 cu = cc[m][k-1];
1717 if( k1 == 0 ) { /* set first angle */
1718 sv = su;
1719 cv = cu;
1720 k1 = 1;
1721 }
1722 else { /* combine angles */
1723 sw = su*cv + cu*sv;
1724 cv = cu*cv - su*sv;
1725 sv = sw;
1726 }
1727 }
1728 }
1729 /* longitude coefficient, in 0.0001" */
1730 f = p[5] * 0.0001;
1731 if( p[6] != 0 )
1732 f += 0.00001 * T * p[6];
1733 /* obliquity coefficient, in 0.0001" */
1734 g = p[7] * 0.0001;
1735 if( p[8] != 0 )
1736 g += 0.00001 * T * p[8];
1737 if (*p >= 100) { /* coefficients in 0.00001" */
1738 f *= 0.1;
1739 g *= 0.1;
1740 }
1741 /* accumulate the terms */
1742 if (*p != 102) {
1743 C += f * sv;
1744 D += g * cv;
1745 }
1746 else { /* cos for nutl and sin for nuto */
1747 C += f * cv;
1748 D += g * sv;
1749 }
1750 /*
1751 if (i >= 105) {
1752 printf("%4.10f, %4.10f\n",f*sv,g*cv);
1753 }
1754 */
1755 }
1756 /*
1757 printf("%4.10f, %4.10f, %4.10f, %4.10f\n",MS*RADTODEG,FF*RADTODEG,DD*RADTODEG,OM*RADTODEG);
1758 printf( "nutation: in longitude %.9f\", in obliquity %.9f\"\n", C, D );
1759 */
1760 /* Save answers, expressed in radians */
1761 nutlo[0] = DEGTORAD * C / 3600.0;
1762 nutlo[1] = DEGTORAD * D / 3600.0;
1763 return(0);
1764 }
1765
1766 /* Nutation IAU 2000A model
1767 * (MHB2000 luni-solar and planetary nutation, without free core nutation)
1768 *
1769 * Function returns nutation in longitude and obliquity in radians with
1770 * respect to the equinox of date. For the obliquity of the ecliptic
1771 * the calculation of Lieske & al. (1977) must be used.
1772 *
1773 * The precision in recent years is about 0.001 arc seconds.
1774 *
1775 * The calculation includes luni-solar and planetary nutation.
1776 * Free core nutation, which cannot be predicted, is omitted,
1777 * the error being of the order of a few 0.0001 arc seconds.
1778 *
1779 * References:
1780 *
1781 * Capitaine, N., Wallace, P.T., Chapront, J., A & A 432, 366 (2005).
1782 *
1783 * Chapront, J., Chapront-Touze, M. & Francou, G., A & A 387, 700 (2002).
1784 *
1785 * Lieske, J.H., Lederle, T., Fricke, W. & Morando, B., "Expressions
1786 * for the precession quantities based upon the IAU (1976) System of
1787 * Astronomical Constants", A & A 58, 1-16 (1977).
1788 *
1789 * Mathews, P.M., Herring, T.A., Buffet, B.A., "Modeling of nutation
1790 * and precession New nutation series for nonrigid Earth and
1791 * insights into the Earth's interior", J.Geophys.Res., 107, B4,
1792 * 2002.
1793 *
1794 * Simon, J.-L., Bretagnon, P., Chapront, J., Chapront-Touze, M.,
1795 * Francou, G., Laskar, J., A & A 282, 663-683 (1994).
1796 *
1797 * Souchay, J., Loysel, B., Kinoshita, H., Folgueira, M., A & A Supp.
1798 * Ser. 135, 111 (1999).
1799 *
1800 * Wallace, P.T., "Software for Implementing the IAU 2000
1801 * Resolutions", in IERS Workshop 5.1 (2002).
1802 *
1803 * Nutation IAU 2000A series in:
1804 * Kaplan, G.H., United States Naval Observatory Circular No. 179 (Oct. 2005)
1805 * aa.usno.navy.mil/publications/docs/Circular_179.html
1806 *
1807 * MHB2000 code at
1808 * - ftp://maia.usno.navy.mil/conv2000/chapter5/IAU2000A.
1809 * - http://www.iau-sofa.rl.ac.uk/2005_0901/Downloads.html
1810 */
1811
1812 #include "swenut2000a.h"
calc_nutation_iau2000ab(double J,double * nutlo)1813 static int calc_nutation_iau2000ab(double J, double *nutlo)
1814 {
1815 int i, j, k, inls;
1816 double M, SM, F, D, OM;
1817 double AL, ALSU, AF, AD, AOM, APA;
1818 double ALME, ALVE, ALEA, ALMA, ALJU, ALSA, ALUR, ALNE;
1819 double darg, sinarg, cosarg;
1820 double dpsi = 0, deps = 0;
1821 double T = (J - J2000 ) / 36525.0;
1822 int nut_model = swed.astro_models[SE_MODEL_NUT];
1823 if (nut_model == 0) nut_model = SEMOD_NUT_DEFAULT;
1824 /* luni-solar nutation */
1825 /* Fundamental arguments, Simon & al. (1994) */
1826 /* Mean anomaly of the Moon. */
1827 M = swe_degnorm(( 485868.249036 +
1828 T*( 1717915923.2178 +
1829 T*( 31.8792 +
1830 T*( 0.051635 +
1831 T*( - 0.00024470 ))))) / 3600.0) * DEGTORAD;
1832 /* Mean anomaly of the Sun */
1833 SM = swe_degnorm((1287104.79305 +
1834 T*( 129596581.0481 +
1835 T*( - 0.5532 +
1836 T*( 0.000136 +
1837 T*( - 0.00001149 ))))) / 3600.0) * DEGTORAD;
1838 /* Mean argument of the latitude of the Moon. */
1839 F = swe_degnorm(( 335779.526232 +
1840 T*( 1739527262.8478 +
1841 T*( - 12.7512 +
1842 T*( - 0.001037 +
1843 T*( 0.00000417 ))))) / 3600.0) * DEGTORAD;
1844 /* Mean elongation of the Moon from the Sun. */
1845 D = swe_degnorm((1072260.70369 +
1846 T*( 1602961601.2090 +
1847 T*( - 6.3706 +
1848 T*( 0.006593 +
1849 T*( - 0.00003169 ))))) / 3600.0) * DEGTORAD;
1850 /* Mean longitude of the ascending node of the Moon. */
1851 OM = swe_degnorm(( 450160.398036 +
1852 T*( - 6962890.5431 +
1853 T*( 7.4722 +
1854 T*( 0.007702 +
1855 T*( - 0.00005939 ))))) / 3600.0) * DEGTORAD;
1856 /* luni-solar nutation series, in reverse order, starting with small terms */
1857 if (nut_model == SEMOD_NUT_IAU_2000B)
1858 inls = NLS_2000B;
1859 else
1860 inls = NLS;
1861 for (i = inls - 1; i >= 0; i--) {
1862 j = i * 5;
1863 darg = swe_radnorm((double) nls[j + 0] * M +
1864 (double) nls[j + 1] * SM +
1865 (double) nls[j + 2] * F +
1866 (double) nls[j + 3] * D +
1867 (double) nls[j + 4] * OM);
1868 sinarg = sin(darg);
1869 cosarg = cos(darg);
1870 k = i * 6;
1871 dpsi += (cls[k+0] + cls[k+1] * T) * sinarg + cls[k+2] * cosarg;
1872 deps += (cls[k+3] + cls[k+4] * T) * cosarg + cls[k+5] * sinarg;
1873 }
1874 nutlo[0] = dpsi * O1MAS2DEG;
1875 nutlo[1] = deps * O1MAS2DEG;
1876 if (nut_model == SEMOD_NUT_IAU_2000A) {
1877 /* planetary nutation
1878 * note: The MHB2000 code computes the luni-solar and planetary nutation
1879 * in different routines, using slightly different Delaunay
1880 * arguments in the two cases. This behaviour is faithfully
1881 * reproduced here. Use of the Simon et al. expressions for both
1882 * cases leads to negligible changes, well below 0.1 microarcsecond.*/
1883 /* Mean anomaly of the Moon.*/
1884 AL = swe_radnorm(2.35555598 + 8328.6914269554 * T);
1885 /* Mean anomaly of the Sun.*/
1886 ALSU = swe_radnorm(6.24006013 + 628.301955 * T);
1887 /* Mean argument of the latitude of the Moon. */
1888 AF = swe_radnorm(1.627905234 + 8433.466158131 * T);
1889 /* Mean elongation of the Moon from the Sun. */
1890 AD = swe_radnorm(5.198466741 + 7771.3771468121 * T);
1891 /* Mean longitude of the ascending node of the Moon. */
1892 AOM = swe_radnorm(2.18243920 - 33.757045 * T);
1893 /* Planetary longitudes, Mercury through Neptune (Souchay et al. 1999). */
1894 ALME = swe_radnorm(4.402608842 + 2608.7903141574 * T);
1895 ALVE = swe_radnorm(3.176146697 + 1021.3285546211 * T);
1896 ALEA = swe_radnorm(1.753470314 + 628.3075849991 * T);
1897 ALMA = swe_radnorm(6.203480913 + 334.0612426700 * T);
1898 ALJU = swe_radnorm(0.599546497 + 52.9690962641 * T);
1899 ALSA = swe_radnorm(0.874016757 + 21.3299104960 * T);
1900 ALUR = swe_radnorm(5.481293871 + 7.4781598567 * T);
1901 ALNE = swe_radnorm(5.321159000 + 3.8127774000 * T);
1902 /* General accumulated precession in longitude. */
1903 APA = (0.02438175 + 0.00000538691 * T) * T;
1904 /* planetary nutation series (in reverse order).*/
1905 dpsi = 0;
1906 deps = 0;
1907 for (i = NPL - 1; i >= 0; i--) {
1908 j = i * 14;
1909 darg = swe_radnorm((double) npl[j + 0] * AL +
1910 (double) npl[j + 1] * ALSU +
1911 (double) npl[j + 2] * AF +
1912 (double) npl[j + 3] * AD +
1913 (double) npl[j + 4] * AOM +
1914 (double) npl[j + 5] * ALME +
1915 (double) npl[j + 6] * ALVE +
1916 (double) npl[j + 7] * ALEA +
1917 (double) npl[j + 8] * ALMA +
1918 (double) npl[j + 9] * ALJU +
1919 (double) npl[j +10] * ALSA +
1920 (double) npl[j +11] * ALUR +
1921 (double) npl[j +12] * ALNE +
1922 (double) npl[j +13] * APA);
1923 k = i * 4;
1924 sinarg = sin(darg);
1925 cosarg = cos(darg);
1926 dpsi += (double) icpl[k+0] * sinarg + (double) icpl[k+1] * cosarg;
1927 deps += (double) icpl[k+2] * sinarg + (double) icpl[k+3] * cosarg;
1928 }
1929 nutlo[0] += dpsi * O1MAS2DEG;
1930 nutlo[1] += deps * O1MAS2DEG;
1931 #if 1
1932 /* changes required by adoption of P03 precession
1933 * according to Capitaine et al. A & A 412, 366 (2005) = IAU 2006 */
1934 dpsi = -8.1 * sin(OM) - 0.6 * sin(2 * F - 2 * D + 2 * OM);
1935 dpsi += T * (47.8 * sin(OM) + 3.7 * sin(2 * F - 2 * D + 2 * OM) + 0.6 * sin(2 * F + 2 * OM) - 0.6 * sin(2 * OM));
1936 deps = T * (-25.6 * cos(OM) - 1.6 * cos(2 * F - 2 * D + 2 * OM));
1937 nutlo[0] += dpsi / (3600.0 * 1000000.0);
1938 nutlo[1] += deps / (3600.0 * 1000000.0);
1939 #endif
1940 }
1941 nutlo[0] *= DEGTORAD;
1942 nutlo[1] *= DEGTORAD;
1943 return 0;
1944 }
1945
1946 /* an incomplete implementation of nutation Woolard 1953 */
calc_nutation_woolard(double J,double * nutlo)1947 static int calc_nutation_woolard(double J, double *nutlo)
1948 {
1949 double deps, dpsi;
1950 double ls, ld; /* sun's mean longitude, moon's mean longitude */
1951 double ms, md; /* sun's mean anomaly, moon's mean anomaly */
1952 double nm; /* longitude of moon's ascending node */
1953 double t, t2; /* number of Julian centuries of 36525 days since
1954 * Jan 0.5 1900.
1955 */
1956 double tls, tnm, tld; /* twice above */
1957 double a, b; /* temps */
1958 double mjd = J - J1900;
1959 t = mjd/36525.;
1960 t2 = t*t;
1961 a = 100.0021358*t;
1962 b = 360.*(a-(long)a);
1963 ls = 279.697+.000303*t2+b;
1964 a = 1336.855231*t;
1965 b = 360.*(a-(long)a);
1966 ld = 270.434-.001133*t2+b;
1967 a = 99.99736056000026*t;
1968 b = 360.*(a-(long)a);
1969 ms = 358.476-.00015*t2+b;
1970 a = 13255523.59*t;
1971 b = 360.*(a-(long)a);
1972 md = 296.105+.009192*t2+b;
1973 a = 5.372616667*t;
1974 b = 360.*(a-(long)a);
1975 nm = 259.183+.002078*t2-b;
1976 /* convert to radian forms for use with trig functions.
1977 */
1978 tls = 2*ls * DEGTORAD;
1979 nm = nm * DEGTORAD;
1980 tnm = 2*nm;
1981 ms = ms * DEGTORAD;
1982 tld = 2*ld * DEGTORAD;
1983 md = md * DEGTORAD;
1984 /* find delta psi and eps, in arcseconds.
1985 */
1986 dpsi = (-17.2327-.01737*t)*sin(nm)+(-1.2729-.00013*t)*sin(tls)
1987 +.2088*sin(tnm)-.2037*sin(tld)+(.1261-.00031*t)*sin(ms)
1988 +.0675*sin(md)-(.0497-.00012*t)*sin(tls+ms)
1989 -.0342*sin(tld-nm)-.0261*sin(tld+md)+.0214*sin(tls-ms)
1990 -.0149*sin(tls-tld+md)+.0124*sin(tls-nm)+.0114*sin(tld-md);
1991 deps = (9.21+.00091*t)*cos(nm)+(.5522-.00029*t)*cos(tls)
1992 -.0904*cos(tnm)+.0884*cos(tld)+.0216*cos(tls+ms)
1993 +.0183*cos(tld-nm)+.0113*cos(tld+md)-.0093*cos(tls-ms)
1994 -.0066*cos(tls-nm);
1995 /* convert to radians.
1996 */
1997 dpsi = dpsi/3600.0 * DEGTORAD;
1998 deps = deps/3600.0 * DEGTORAD;
1999 nutlo[1] = deps;
2000 nutlo[0] = dpsi;
2001 return OK;
2002 }
2003
bessel(double * v,int n,double t)2004 static double bessel(double *v, int n, double t)
2005 {
2006 int i, iy, k;
2007 double ans, p, B, d[6];
2008 if (t <= 0) {
2009 ans = v[0];
2010 goto done;
2011 }
2012 if (t >= n - 1) {
2013 ans = v[n - 1];
2014 goto done;
2015 }
2016 p = floor(t);
2017 iy = (int) t;
2018 /* Zeroth order estimate is value at start of year */
2019 ans = v[iy];
2020 k = iy + 1;
2021 if (k >= n)
2022 goto done;
2023 /* The fraction of tabulation interval */
2024 p = t - p;
2025 ans += p * (v[k] - v[iy]);
2026 if( (iy - 1 < 0) || (iy + 2 >= n) )
2027 goto done; /* can't do second differences */
2028 /* Make table of first differences */
2029 k = iy - 2;
2030 for (i = 0; i < 5; i++) {
2031 if((k < 0) || (k + 1 >= n))
2032 d[i] = 0;
2033 else
2034 d[i] = v[k+1] - v[k];
2035 k += 1;
2036 }
2037 /* Compute second differences */
2038 for (i = 0; i < 4; i++ )
2039 d[i] = d[i+1] - d[i];
2040 B = 0.25 * p * (p - 1.0);
2041 ans += B * (d[1] + d[2]);
2042 #if DEMO
2043 printf("B %.4lf, ans %.4lf\n", B, ans);
2044 #endif
2045 if (iy + 2 >= n)
2046 goto done;
2047 /* Compute third differences */
2048 for (i = 0; i < 3; i++ )
2049 d[i] = d[i + 1] - d[i];
2050 B = 2.0 * B / 3.0;
2051 ans += (p - 0.5) * B * d[1];
2052 #if DEMO
2053 printf("B %.4lf, ans %.4lf\n", B * (p - 0.5), ans);
2054 #endif
2055 if ((iy - 2 < 0) || (iy + 3 > n))
2056 goto done;
2057 /* Compute fourth differences */
2058 for (i = 0; i < 2; i++)
2059 d[i] = d[i + 1] - d[i];
2060 B = 0.125 * B * (p + 1.0) * (p - 2.0);
2061 ans += B * (d[0] + d[1]);
2062 #if DEMO
2063 printf("B %.4lf, ans %.4lf\n", B, ans);
2064 #endif
2065 done:
2066 return ans;
2067 }
2068
calc_nutation(double J,int32 iflag,double * nutlo)2069 static int calc_nutation(double J, int32 iflag, double *nutlo)
2070 {
2071 int n;
2072 double dpsi, deps, J2;
2073 int nut_model = swed.astro_models[SE_MODEL_NUT];
2074 int jplhora_model = swed.astro_models[SE_MODEL_JPLHORA_MODE];
2075 AS_BOOL is_jplhor = FALSE;
2076 if (nut_model == 0) nut_model = SEMOD_NUT_DEFAULT;
2077 if (jplhora_model == 0) jplhora_model = SEMOD_JPLHORA_DEFAULT;
2078 if (iflag & SEFLG_JPLHOR)
2079 is_jplhor = TRUE;
2080 if ((iflag & SEFLG_JPLHOR_APPROX) &&
2081 jplhora_model == SEMOD_JPLHORA_3
2082 && J <= HORIZONS_TJD0_DPSI_DEPS_IAU1980)
2083 is_jplhor = TRUE;
2084 if (is_jplhor) {
2085 calc_nutation_iau1980(J, nutlo);
2086 if (iflag & SEFLG_JPLHOR) {
2087 n = (int) (swed.eop_tjd_end - swed.eop_tjd_beg + 0.000001);
2088 J2 = J;
2089 if (J < swed.eop_tjd_beg_horizons)
2090 J2 = swed.eop_tjd_beg_horizons;
2091 dpsi = bessel(swed.dpsi, n + 1, J2 - swed.eop_tjd_beg);
2092 deps = bessel(swed.deps, n + 1, J2 - swed.eop_tjd_beg);
2093 nutlo[0] += dpsi / 3600.0 * DEGTORAD;
2094 nutlo[1] += deps / 3600.0 * DEGTORAD;
2095 #if 0
2096 printf("tjd=%f, dpsi=%f, deps=%f\n", J, dpsi * 1000, deps * 1000);
2097 #endif
2098 } else {
2099 nutlo[0] += DPSI_IAU1980_TJD0 / 3600.0 * DEGTORAD;
2100 nutlo[1] += DEPS_IAU1980_TJD0 / 3600.0 * DEGTORAD;
2101 }
2102 } else if (nut_model == SEMOD_NUT_IAU_1980 || nut_model == SEMOD_NUT_IAU_CORR_1987) {
2103 calc_nutation_iau1980(J, nutlo);
2104 } else if (nut_model == SEMOD_NUT_IAU_2000A || nut_model == SEMOD_NUT_IAU_2000B) {
2105 calc_nutation_iau2000ab(J, nutlo);
2106 if ((iflag & SEFLG_JPLHOR_APPROX) && jplhora_model == SEMOD_JPLHORA_2) {
2107 nutlo[0] += -41.7750 / 3600.0 / 1000.0 * DEGTORAD;
2108 nutlo[1] += -6.8192 / 3600.0 / 1000.0 * DEGTORAD;
2109 }
2110 } else if (nut_model == SEMOD_NUT_WOOLARD) {
2111 calc_nutation_woolard(J, nutlo);
2112 }
2113 return OK;
2114 }
2115
quadratic_intp(double ym,double y0,double yp,double x)2116 static double quadratic_intp(double ym, double y0, double yp, double x)
2117 {
2118 double a, b, c, y;
2119 c = y0;
2120 b = (yp - ym) / 2.0;
2121 a = (yp + ym) / 2.0 - c;
2122 y = a * x * x + b * x + c;
2123 return y;
2124 }
2125
swi_nutation(double tjd,int32 iflag,double * nutlo)2126 int swi_nutation(double tjd, int32 iflag, double *nutlo)
2127 {
2128 int retc = OK;
2129 double dnut[2], dx;
2130 if (!swed.do_interpolate_nut) {
2131 retc = calc_nutation(tjd, iflag, nutlo);
2132 // from interpolation, with three data points in 1-day steps;
2133 // maximum error is about 3 mas
2134 } else {
2135 // precalculated data points available
2136 if (tjd < swed.interpol.tjd_nut2 && tjd > swed.interpol.tjd_nut0) {
2137 dx = (tjd - swed.interpol.tjd_nut0) - 1.0;
2138 nutlo[0] = quadratic_intp(swed.interpol.nut_dpsi0, swed.interpol.nut_dpsi1, swed.interpol.nut_dpsi2, dx);
2139 nutlo[1] = quadratic_intp(swed.interpol.nut_deps0, swed.interpol.nut_deps1, swed.interpol.nut_deps2, dx);
2140 } else {
2141 swed.interpol.tjd_nut0 = tjd - 1.0; // one day earlier
2142 swed.interpol.tjd_nut2 = tjd + 1.0; // one day later
2143 retc = calc_nutation(swed.interpol.tjd_nut0, iflag, dnut);
2144 if (retc == ERR) return ERR;
2145 swed.interpol.nut_dpsi0 = dnut[0];
2146 swed.interpol.nut_deps0 = dnut[1];
2147 retc = calc_nutation(swed.interpol.tjd_nut2, iflag, dnut);
2148 if (retc == ERR) return ERR;
2149 swed.interpol.nut_dpsi2 = dnut[0];
2150 swed.interpol.nut_deps2 = dnut[1];
2151 retc = calc_nutation(tjd, iflag, nutlo);
2152 if (retc == ERR) return ERR;
2153 swed.interpol.nut_dpsi1 = nutlo[0];
2154 swed.interpol.nut_deps1 = nutlo[1];
2155 }
2156 }
2157 return retc;
2158 }
2159
2160 #define OFFSET_JPLHORIZONS (-52.3)
2161 #define DCOR_RA_JPL_TJD0 2437846.5
2162 #define NDCOR_RA_JPL 51
2163 double dcor_ra_jpl[] = {
2164 -51.257, -51.103, -51.065, -51.503, -51.224, -50.796, -51.161, -51.181,
2165 -50.932, -51.064, -51.182, -51.386, -51.416, -51.428, -51.586, -51.766, -52.038, -52.370,
2166 -52.553, -52.397, -52.340, -52.676, -52.348, -51.964, -52.444, -52.364, -51.988, -52.212,
2167 -52.370, -52.523, -52.541, -52.496, -52.590, -52.629, -52.788, -53.014, -53.053, -52.902,
2168 -52.850, -53.087, -52.635, -52.185, -52.588, -52.292, -51.796, -51.961, -52.055, -52.134,
2169 -52.165, -52.141, -52.255,
2170 };
2171
swi_approx_jplhor(double * x,double tjd,int32 iflag,AS_BOOL backward)2172 static void swi_approx_jplhor(double *x, double tjd, int32 iflag, AS_BOOL backward)
2173 {
2174 double t0, t1;
2175 double t = (tjd - DCOR_RA_JPL_TJD0) / 365.25;
2176 double dofs = OFFSET_JPLHORIZONS;
2177 int jplhora_model = swed.astro_models[SE_MODEL_JPLHORA_MODE];
2178 if (jplhora_model == 0) jplhora_model = SEMOD_JPLHORA_DEFAULT;
2179 if (!(iflag & SEFLG_JPLHOR_APPROX))
2180 return;
2181 if (jplhora_model == SEMOD_JPLHORA_2)
2182 return;
2183 if (t < 0) {
2184 t = 0;
2185 dofs = dcor_ra_jpl[0];
2186 } else if (t >= NDCOR_RA_JPL - 1) {
2187 t = NDCOR_RA_JPL;
2188 dofs = dcor_ra_jpl[NDCOR_RA_JPL - 1];
2189 } else {
2190 t0 = (int) t;
2191 t1 = t0 + 1;
2192 dofs = dcor_ra_jpl[(int)t0];
2193 dofs = (t - t0) * (dcor_ra_jpl[(int)t0] - dcor_ra_jpl[(int)t1]) + dcor_ra_jpl[(int)t0];
2194 }
2195 dofs /= (1000.0 * 3600.0);
2196 swi_cartpol(x, x);
2197 if (backward)
2198 x[0] -= dofs * DEGTORAD;
2199 else
2200 x[0] += dofs * DEGTORAD;
2201 swi_polcart(x, x);
2202 }
2203
2204 /* GCRS to J2000 */
swi_bias(double * x,double tjd,int32 iflag,AS_BOOL backward)2205 void swi_bias(double *x, double tjd, int32 iflag, AS_BOOL backward)
2206 {
2207 #if 0
2208 double DAS2R = 1.0 / 3600.0 * DEGTORAD;
2209 double dpsi_bias = -0.041775 * DAS2R;
2210 double deps_bias = -0.0068192 * DAS2R;
2211 double dra0 = -0.0146 * DAS2R;
2212 double deps2000 = 84381.448 * DAS2R;
2213 #endif
2214 double xx[6], rb[3][3];
2215 int i;
2216 int bias_model = swed.astro_models[SE_MODEL_BIAS];
2217 int jplhora_model = swed.astro_models[SE_MODEL_JPLHORA_MODE];
2218 if (bias_model == 0) bias_model = SEMOD_BIAS_DEFAULT;
2219 if (jplhora_model == 0) jplhora_model = SEMOD_JPLHORA_DEFAULT;
2220 if (bias_model == SEMOD_BIAS_NONE)
2221 return;
2222 if (iflag & SEFLG_JPLHOR_APPROX) {
2223 if (jplhora_model == SEMOD_JPLHORA_2)
2224 return;
2225 if (jplhora_model == SEMOD_JPLHORA_3 && tjd < DPSI_DEPS_IAU1980_TJD0_HORIZONS)
2226 return;
2227 }
2228 /* #if FRAME_BIAS_IAU2006 * frame bias 2006 */
2229 if (bias_model == SEMOD_BIAS_IAU2006) {
2230 rb[0][0] = +0.99999999999999412;
2231 rb[1][0] = -0.00000007078368961;
2232 rb[2][0] = +0.00000008056213978;
2233 rb[0][1] = +0.00000007078368695;
2234 rb[1][1] = +0.99999999999999700;
2235 rb[2][1] = +0.00000003306428553;
2236 rb[0][2] = -0.00000008056214212;
2237 rb[1][2] = -0.00000003306427981;
2238 rb[2][2] = +0.99999999999999634;
2239 /* #else * frame bias 2000, makes no difference in result */
2240 } else {
2241 rb[0][0] = +0.9999999999999942;
2242 rb[1][0] = -0.0000000707827974;
2243 rb[2][0] = +0.0000000805621715;
2244 rb[0][1] = +0.0000000707827948;
2245 rb[1][1] = +0.9999999999999969;
2246 rb[2][1] = +0.0000000330604145;
2247 rb[0][2] = -0.0000000805621738;
2248 rb[1][2] = -0.0000000330604088;
2249 rb[2][2] = +0.9999999999999962;
2250 }
2251 /*#endif*/
2252 #if 0
2253 rb[0][0] = +0.9999999999999968;
2254 rb[1][0] = +0.0000000000000000;
2255 rb[2][0] = +0.0000000805621715;
2256 rb[0][1] = -0.0000000000000027;
2257 rb[1][1] = +0.9999999999999994;
2258 rb[2][1] = +0.0000000330604145;
2259 rb[0][2] = -0.0000000805621715;
2260 rb[1][2] = -0.0000000330604145;
2261 rb[2][2] = +0.9999999999999962;
2262 #endif
2263 if (backward) {
2264 swi_approx_jplhor(x, tjd, iflag, TRUE);
2265 for (i = 0; i <= 2; i++) {
2266 xx[i] = x[0] * rb[i][0] +
2267 x[1] * rb[i][1] +
2268 x[2] * rb[i][2];
2269 if (iflag & SEFLG_SPEED)
2270 xx[i+3] = x[3] * rb[i][0] +
2271 x[4] * rb[i][1] +
2272 x[5] * rb[i][2];
2273 }
2274 } else {
2275 for (i = 0; i <= 2; i++) {
2276 xx[i] = x[0] * rb[0][i] +
2277 x[1] * rb[1][i] +
2278 x[2] * rb[2][i];
2279 if (iflag & SEFLG_SPEED)
2280 xx[i+3] = x[3] * rb[0][i] +
2281 x[4] * rb[1][i] +
2282 x[5] * rb[2][i];
2283 }
2284 swi_approx_jplhor(xx, tjd, iflag, FALSE);
2285 }
2286 for (i = 0; i <= 2; i++) x[i] = xx[i];
2287 if (iflag & SEFLG_SPEED)
2288 for (i = 3; i <= 5; i++) x[i] = xx[i];
2289 }
2290
2291 /* GCRS to FK5 */
swi_icrs2fk5(double * x,int32 iflag,AS_BOOL backward)2292 void swi_icrs2fk5(double *x, int32 iflag, AS_BOOL backward)
2293 {
2294 #if 0
2295 double DAS2R = 1.0 / 3600.0 * DEGTORAD;
2296 double dra0 = -0.0229 * DAS2R;
2297 double dxi0 = 0.0091 * DAS2R;
2298 double det0 = -0.0199 * DAS2R;
2299 #endif
2300 double xx[6], rb[3][3];
2301 int i;
2302 rb[0][0] = +0.9999999999999928;
2303 rb[0][1] = +0.0000001110223287;
2304 rb[0][2] = +0.0000000441180557;
2305 rb[1][0] = -0.0000001110223330;
2306 rb[1][1] = +0.9999999999999891;
2307 rb[1][2] = +0.0000000964779176;
2308 rb[2][0] = -0.0000000441180450;
2309 rb[2][1] = -0.0000000964779225;
2310 rb[2][2] = +0.9999999999999943;
2311 if (backward) {
2312 for (i = 0; i <= 2; i++) {
2313 xx[i] = x[0] * rb[i][0] +
2314 x[1] * rb[i][1] +
2315 x[2] * rb[i][2];
2316 if (iflag & SEFLG_SPEED)
2317 xx[i+3] = x[3] * rb[i][0] +
2318 x[4] * rb[i][1] +
2319 x[5] * rb[i][2];
2320 }
2321 } else {
2322 for (i = 0; i <= 2; i++) {
2323 xx[i] = x[0] * rb[0][i] +
2324 x[1] * rb[1][i] +
2325 x[2] * rb[2][i];
2326 if (iflag & SEFLG_SPEED)
2327 xx[i+3] = x[3] * rb[0][i] +
2328 x[4] * rb[1][i] +
2329 x[5] * rb[2][i];
2330 }
2331 }
2332 for (i = 0; i <= 5; i++) x[i] = xx[i];
2333 }
2334
2335 /* DeltaT = Ephemeris Time - Universal Time, in days.
2336 *
2337 * Before 1955 we use the data developed by
2338 * Stephenson, Morrison, and Hohenkerk (2016),
2339 *
2340 * 1955 - today + a couple of years:
2341 * ---------------------------------
2342 * The tabulated values of deltaT from the Astronomical
2343 * Alamanc (AA 1997 etc. pp. K8-K9) are used. Some
2344 * more recent values have been taken from IERS
2345 * (http://maia.usno.navy.mil/ser7/deltat.data).
2346 * Bessel's interpolation formula is implemented to obtain fourth
2347 * order interpolated values at intermediate times.
2348 * The values are adjusted depending on the ephemeris used
2349 * and its inherent value of secular tidal acceleration ndot.
2350 *
2351 * future:
2352 * ---------------------------------
2353 * For the time after the last tabulated value, we use the formula
2354 * of Stephenson (1997; p. 507), with a modification that avoids a jump
2355 * at the end of the tabulated period. A linear term is added that
2356 * makes a slow transition from the table to the formula over a period
2357 * of 100 years. (Need not be updated, when table will be enlarged.)
2358 *
2359 * References:
2360 *
2361 * Stephenson, F. R., and L. V. Morrison, "Long-term changes
2362 * in the rotation of the Earth: 700 B.C. to A.D. 1980,"
2363 * Philosophical Transactions of the Royal Society of London
2364 * Series A 313, 47-70 (1984)
2365 *
2366 * Borkowski, K. M., "ELP2000-85 and the Dynamical Time
2367 * - Universal Time relation," Astronomy and Astrophysics
2368 * 205, L8-L10 (1988)
2369 * Borkowski's formula is derived from partly doubtful eclipses
2370 * going back to 2137 BC and uses lunar position based on tidal
2371 * coefficient of -23.9 arcsec/cy^2.
2372 *
2373 * Chapront-Touze, Michelle, and Jean Chapront, _Lunar Tables
2374 * and Programs from 4000 B.C. to A.D. 8000_, Willmann-Bell 1991
2375 * Their table agrees with the one here, but the entries are
2376 * rounded to the nearest whole second.
2377 *
2378 * Stephenson, F. R., and M. A. Houlden, _Atlas of Historical
2379 * Eclipse Maps_, Cambridge U. Press (1986)
2380 *
2381 * Stephenson, F.R. & Morrison, L.V., "Long-Term Fluctuations in
2382 * the Earth's Rotation: 700 BC to AD 1990", Philosophical
2383 * Transactions of the Royal Society of London,
2384 * Ser. A, 351 (1995), 165-202.
2385 *
2386 * Stephenson, F. Richard, _Historical Eclipses and Earth's
2387 * Rotation_, Cambridge U. Press (1997)
2388 *
2389 * Morrison, L. V., and F.R. Stephenson, "Historical Values of the Earth's
2390 * Clock Error DT and the Calculation of Eclipses", JHA xxxv (2004),
2391 * pp.327-336
2392 *
2393 * Stephenson, F.R., Morrison, L.V., and Hohenkerk, C.Y., "Measurement of the
2394 * Earth's Rotation: 720 BC to AD 2015", Royal Society Proceedings A
2395 * 7 Dec 2016,
2396 * http://rspa.royalsocietypublishing.org/lookup/doi/10.1098/rspa.2016.0404
2397 *
2398 * Table from AA for 1620 through today
2399 * Note, Stephenson and Morrison's table starts at the year 1630.
2400 * The Chapronts' table does not agree with the Almanac prior to 1630.
2401 * The actual accuracy decreases rapidly prior to 1780.
2402 *
2403 * Jean Meeus, Astronomical Algorithms, 2nd edition, 1998.
2404 *
2405 * For a comprehensive collection of publications and formulae, see:
2406 * http://www.phys.uu.nl/~vgent/deltat/deltat_modern.htm
2407 * http://www.phys.uu.nl/~vgent/deltat/deltat_old.htm
2408 *
2409 * For future values of delta t, the following data from the
2410 * Earth Orientation Department of the US Naval Observatory can be used:
2411 * (TAI-UTC) from: ftp://maia.usno.navy.mil/ser7/tai-utc.dat
2412 * (UT1-UTC) from: ftp://maia.usno.navy.mil/ser7/finals.all (cols. 59-68)
2413 * or: ftp://ftp.iers.org/products/eop/rapid/standard/finals.data
2414 * file description in: ftp://maia.usno.navy.mil/ser7/readme.finals
2415 * Delta T = TAI-UT1 + 32.184 sec = (TAI-UTC) - (UT1-UTC) + 32.184 sec
2416 *
2417 * Also, there is the following file:
2418 * http://maia.usno.navy.mil/ser7/deltat.data, but it is about 3 months
2419 * behind (on 3 feb 2009); and predictions:
2420 * http://maia.usno.navy.mil/ser7/deltat.preds
2421 *
2422 * Last update of table dt[]: Dieter Koch, 18 dec 2013.
2423 * ATTENTION: Whenever updating this table, do not forget to adjust
2424 * the macros TABEND and TABSIZ !
2425 */
2426 #define TABSTART 1620
2427 #define TABEND 2028
2428 #define TABSIZ (TABEND-TABSTART+1)
2429 /* we make the table greater for additional values read from external file */
2430 #define TABSIZ_SPACE (TABSIZ+100)
2431 static TLS double dt[TABSIZ_SPACE] = {
2432 /* 1620.0 - 1659.0 */
2433 124.00, 119.00, 115.00, 110.00, 106.00, 102.00, 98.00, 95.00, 91.00, 88.00,
2434 85.00, 82.00, 79.00, 77.00, 74.00, 72.00, 70.00, 67.00, 65.00, 63.00,
2435 62.00, 60.00, 58.00, 57.00, 55.00, 54.00, 53.00, 51.00, 50.00, 49.00,
2436 48.00, 47.00, 46.00, 45.00, 44.00, 43.00, 42.00, 41.00, 40.00, 38.00,
2437 /* 1660.0 - 1699.0 */
2438 37.00, 36.00, 35.00, 34.00, 33.00, 32.00, 31.00, 30.00, 28.00, 27.00,
2439 26.00, 25.00, 24.00, 23.00, 22.00, 21.00, 20.00, 19.00, 18.00, 17.00,
2440 16.00, 15.00, 14.00, 14.00, 13.00, 12.00, 12.00, 11.00, 11.00, 10.00,
2441 10.00, 10.00, 9.00, 9.00, 9.00, 9.00, 9.00, 9.00, 9.00, 9.00,
2442 /* 1700.0 - 1739.0 */
2443 9.00, 9.00, 9.00, 9.00, 9.00, 9.00, 9.00, 9.00, 10.00, 10.00,
2444 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 11.00, 11.00, 11.00,
2445 11.00, 11.00, 11.00, 11.00, 11.00, 11.00, 11.00, 11.00, 11.00, 11.00,
2446 11.00, 11.00, 11.00, 11.00, 12.00, 12.00, 12.00, 12.00, 12.00, 12.00,
2447 /* 1740.0 - 1779.0 */
2448 12.00, 12.00, 12.00, 12.00, 13.00, 13.00, 13.00, 13.00, 13.00, 13.00,
2449 13.00, 14.00, 14.00, 14.00, 14.00, 14.00, 14.00, 14.00, 15.00, 15.00,
2450 15.00, 15.00, 15.00, 15.00, 15.00, 16.00, 16.00, 16.00, 16.00, 16.00,
2451 16.00, 16.00, 16.00, 16.00, 16.00, 17.00, 17.00, 17.00, 17.00, 17.00,
2452 /* 1780.0 - 1799.0 */
2453 17.00, 17.00, 17.00, 17.00, 17.00, 17.00, 17.00, 17.00, 17.00, 17.00,
2454 17.00, 17.00, 16.00, 16.00, 16.00, 16.00, 15.00, 15.00, 14.00, 14.00,
2455 /* 1800.0 - 1819.0 */
2456 13.70, 13.40, 13.10, 12.90, 12.70, 12.60, 12.50, 12.50, 12.50, 12.50,
2457 12.50, 12.50, 12.50, 12.50, 12.50, 12.50, 12.50, 12.40, 12.30, 12.20,
2458 /* 1820.0 - 1859.0 */
2459 12.00, 11.70, 11.40, 11.10, 10.60, 10.20, 9.60, 9.10, 8.60, 8.00,
2460 7.50, 7.00, 6.60, 6.30, 6.00, 5.80, 5.70, 5.60, 5.60, 5.60,
2461 5.70, 5.80, 5.90, 6.10, 6.20, 6.30, 6.50, 6.60, 6.80, 6.90,
2462 7.10, 7.20, 7.30, 7.40, 7.50, 7.60, 7.70, 7.70, 7.80, 7.80,
2463 /* 1860.0 - 1899.0 */
2464 7.88, 7.82, 7.54, 6.97, 6.40, 6.02, 5.41, 4.10, 2.92, 1.82,
2465 1.61, .10, -1.02, -1.28, -2.69, -3.24, -3.64, -4.54, -4.71, -5.11,
2466 -5.40, -5.42, -5.20, -5.46, -5.46, -5.79, -5.63, -5.64, -5.80, -5.66,
2467 -5.87, -6.01, -6.19, -6.64, -6.44, -6.47, -6.09, -5.76, -4.66, -3.74,
2468 /* 1900.0 - 1939.0 */
2469 -2.72, -1.54, -.02, 1.24, 2.64, 3.86, 5.37, 6.14, 7.75, 9.13,
2470 10.46, 11.53, 13.36, 14.65, 16.01, 17.20, 18.24, 19.06, 20.25, 20.95,
2471 21.16, 22.25, 22.41, 23.03, 23.49, 23.62, 23.86, 24.49, 24.34, 24.08,
2472 24.02, 24.00, 23.87, 23.95, 23.86, 23.93, 23.73, 23.92, 23.96, 24.02,
2473 /* 1940.0 - 1949.0 */
2474 24.33, 24.83, 25.30, 25.70, 26.24, 26.77, 27.28, 27.78, 28.25, 28.71,
2475 /* 1950.0 - 1959.0 */
2476 29.15, 29.57, 29.97, 30.36, 30.72, 31.07, 31.35, 31.68, 32.18, 32.68,
2477 /* 1960.0 - 1969.0 */
2478 33.15, 33.59, 34.00, 34.47, 35.03, 35.73, 36.54, 37.43, 38.29, 39.20,
2479 /* 1970.0 - 1979.0 */
2480 /* from 1974 on values (with 4-digit precision) were calculated from IERS data */
2481 40.18, 41.17, 42.23, 43.37, 44.4841, 45.4761, 46.4567, 47.5214, 48.5344, 49.5862,
2482 /* 1980.0 - 1989.0 */
2483 50.5387, 51.3808, 52.1668, 52.9565, 53.7882, 54.3427, 54.8713, 55.3222, 55.8197, 56.3000,
2484 /* 1990.0 - 1999.0 */
2485 56.8553, 57.5653, 58.3092, 59.1218, 59.9845, 60.7854, 61.6287, 62.2951, 62.9659, 63.4673,
2486 /* 2000.0 - 2009.0 */
2487 63.8285, 64.0908, 64.2998, 64.4734, 64.5736, 64.6876, 64.8452, 65.1464, 65.4574, 65.7768,
2488 /* 2010.0 - 2018.0 */
2489 66.0699, 66.3246, 66.6030, 66.9069, 67.2810, 67.6439, 68.1024, 68.5927, 68.9676, 69.2202,
2490 /* 2020.0 - */
2491 69.3612,
2492 /* Extrapolated values:
2493 * 2021 - 2028 */
2494 69.4271, 70.00, 70.50, 71.00, 71.50, 72.00, 72.50, 73.00,
2495 };
2496
2497 #define TAB2_SIZ 27
2498 #define TAB2_START (-1000)
2499 #define TAB2_END 1600
2500 #define TAB2_STEP 100
2501 #define LTERM_EQUATION_YSTART 1820
2502 #define LTERM_EQUATION_COEFF 32
2503 /* Table for -1000 through 1600, from Morrison & Stephenson (2004). */
2504 static const short dt2[TAB2_SIZ] = {
2505 /*-1000 -900 -800 -700 -600 -500 -400 -300 -200 -100*/
2506 25400,23700,22000,21000,19040,17190,15530,14080,12790,11640,
2507 /* 0 100 200 300 400 500 600 700 800 900*/
2508 10580, 9600, 8640, 7680, 6700, 5710, 4740, 3810, 2960, 2200,
2509 /* 1000 1100 1200 1300 1400 1500 1600, */
2510 1570, 1090, 740, 490, 320, 200, 120,
2511 };
2512
2513 /* Table for -500 through 1600, from Stephenson & Morrison (1995).
2514 *
2515 * The first value for -550 has been added from Borkowski
2516 * in order to make this table fit with the Borkowski formula
2517 * for times before -550.
2518 */
2519 #define TAB97_SIZ 43
2520 #define TAB97_START (-500)
2521 #define TAB97_END (1600)
2522 #define TAB97_STEP (50)
2523 static const short dt97[TAB97_SIZ] = {
2524 /* -500 -450 -400 -350 -300 -250 -200 -150 -100 -50*/
2525 16800,16000,15300,14600,14000,13400,12800,12200,11600,11100,
2526 /* 0 50 100 150 200 250 300 350 400 450*/
2527 10600,10100, 9600, 9100, 8600, 8200, 7700, 7200, 6700, 6200,
2528 /* 500 550 600 650 700 750 800 850 900 950*/
2529 5700, 5200, 4700, 4300, 3800, 3400, 3000, 2600, 2200, 1900,
2530 /* 1000 1050 1100 1150 1200 1250 1300 1350 1400 1450*/
2531 1600, 1350, 1100, 900, 750, 600, 470, 380, 300, 230,
2532 /* 1500 1550 1600 */
2533 180, 140, 110,
2534 };
2535
2536 /* returns DeltaT (ET - UT) in days
2537 * double tjd = julian day in UT
2538 * delta t is adjusted to the tidal acceleration that is compatible
2539 * with the ephemeris flag contained in iflag and with the ephemeris
2540 * files made accessible through swe_set_ephe_path() or swe_set_jplfile().
2541 * If iflag = -1, then the default tidal acceleration is ussed (i.e.
2542 * that of DE431).
2543 */
2544 #define DEMO 0
calc_deltat(double tjd,int32 iflag,double * deltat,char * serr)2545 static int32 calc_deltat(double tjd, int32 iflag, double *deltat, char *serr)
2546 {
2547 double ans = 0;
2548 double B, Y, Ygreg, dd;
2549 int iy;
2550 int32 retc;
2551 int deltat_model = swed.astro_models[SE_MODEL_DELTAT];
2552 double tid_acc;
2553 int32 denum, denumret;
2554 int32 epheflag, otherflag;
2555 //fprintf(stderr, "dmod=%f, %.f\n", (double) deltat_model, (double) SEMOD_DELTAT_DEFAULT);
2556 if (deltat_model == 0) deltat_model = SEMOD_DELTAT_DEFAULT;
2557 epheflag = iflag & SEFLG_EPHMASK;
2558 otherflag = iflag & ~SEFLG_EPHMASK;
2559 /* with iflag == -1, we use default tid_acc */
2560 if (iflag == -1) {
2561 retc = swi_get_tid_acc(tjd, 0, 9999, &denumret, &tid_acc, serr); /* for default tid_acc */
2562 /* otherwise we use tid_acc consistent with epheflag */
2563 } else {
2564 denum = swed.jpldenum;
2565 if (epheflag & SEFLG_SWIEPH) denum = swed.fidat[SEI_FILE_MOON].sweph_denum;
2566 if (swi_init_swed_if_start() == 1 && !(epheflag & SEFLG_MOSEPH)) {
2567 if (serr != NULL)
2568 strcpy(serr, "Please call swe_set_ephe_path() or swe_set_jplfile() before calling swe_deltat_ex()");
2569 retc = swi_set_tid_acc(tjd, epheflag, denum, NULL); /* _set_ saves tid_acc in swed */
2570 } else {
2571 retc = swi_set_tid_acc(tjd, epheflag, denum, serr); /* _set_ saves tid_acc in swed */
2572 }
2573 tid_acc = swed.tid_acc;
2574 }
2575 iflag = otherflag | retc;
2576 Y = 2000.0 + (tjd - J2000)/365.25;
2577 Ygreg = 2000.0 + (tjd - J2000)/365.2425;
2578 /* Model for epochs before 1955, currently default in Swiss Ephemeris:
2579 * Stephenson/Morrison/Hohenkerk 2016
2580 * (we switch over to Astronomical Almanac K8-K9 and IERS at 1 Jan. 1955.
2581 * To make the curve continuous we apply some linear term over
2582 * 1000 days before this date.)
2583 * Delta T according to Stephenson/Morrison/Hohenkerk 2016 is based on
2584 * ancient, medieval, and modern observations of eclipses and occultations.
2585 * Values of Deltat T before 1955 depend on this kind of observations.
2586 * For more recent data we want to use the data provided by IERS
2587 * (or Astronomical Almanac K8-K9).
2588 */
2589 if (deltat_model == SEMOD_DELTAT_STEPHENSON_ETC_2016 && tjd < 2435108.5) { // tjd < 2432521.453645833) {
2590 *deltat = deltat_stephenson_etc_2016(tjd, tid_acc);
2591 if (tjd >= 2434108.5) {
2592 *deltat += (1.0 - (2435108.5 - tjd) / 1000.0) * 0.6610218 / 86400.0;
2593 }
2594 return iflag;
2595 }
2596 /* Model used SE 1.77 - 2.05.01, for epochs before 1633:
2597 * Polynomials by Espenak & Meeus 2006,
2598 * derived from Stephenson & Morrison 2004.
2599 * deltat_model == SEMOD_DELTAT_ESPENAK_MEEUS_2006:
2600 * This method is used only for epochs before 1633. (For more recent
2601 * epochs, we use the data provided by Astronomical Almanac K8-K9.)
2602 */
2603 if (deltat_model == SEMOD_DELTAT_ESPENAK_MEEUS_2006 && tjd < 2317746.13090277789) {
2604 *deltat = deltat_espenak_meeus_1620(tjd, tid_acc);
2605 return iflag;
2606 }
2607 /* delta t model used in SE 1.72 - 1.76:
2608 * Stephenson & Morrison 2004;
2609 * before 1620 */
2610 if (deltat_model == SEMOD_DELTAT_STEPHENSON_MORRISON_2004 && Y < TABSTART) {
2611 // before 1600:
2612 if (Y < TAB2_END) {
2613 *deltat = deltat_stephenson_morrison_2004_1600(tjd, tid_acc);
2614 return iflag;
2615 } else {
2616 /* between 1600 and 1620:
2617 * linear interpolation between
2618 * end of table dt2 and start of table dt */
2619 if (Y >= TAB2_END) {
2620 B = TABSTART - TAB2_END;
2621 iy = (TAB2_END - TAB2_START) / TAB2_STEP;
2622 dd = (Y - TAB2_END) / B;
2623 ans = dt2[iy] + dd * (dt[0] - dt2[iy]);
2624 ans = adjust_for_tidacc(ans, Ygreg, tid_acc, SE_TIDAL_26, FALSE);
2625 *deltat = ans / 86400.0;
2626 return iflag;
2627 }
2628 }
2629 }
2630 /* delta t model used in SE 1.64 - 1.71:
2631 * Stephenson 1997;
2632 * before 1620 */
2633 if (deltat_model == SEMOD_DELTAT_STEPHENSON_1997 && Y < TABSTART) {
2634 // before 1600:
2635 if (Y < TAB97_END) {
2636 *deltat = deltat_stephenson_morrison_1997_1600(tjd, tid_acc);
2637 return iflag;
2638 } else {
2639 /* between 1600 and 1620:
2640 * linear interpolation between
2641 * end of table dt2 and start of table dt */
2642 if (Y >= TAB97_END) {
2643 B = TABSTART - TAB97_END;
2644 iy = (TAB97_END - TAB97_START) / TAB97_STEP;
2645 dd = (Y - TAB97_END) / B;
2646 ans = dt97[iy] + dd * (dt[0] - dt97[iy]);
2647 ans = adjust_for_tidacc(ans, Ygreg, tid_acc, SE_TIDAL_26, FALSE);
2648 *deltat = ans / 86400.0;
2649 return iflag;
2650 }
2651 }
2652 }
2653 /* delta t model used before SE 1.64:
2654 * Stephenson/Morrison 1984 with Borkowski 1988;
2655 * before 1620 */
2656 if (deltat_model == SEMOD_DELTAT_STEPHENSON_MORRISON_1984 && Y < TABSTART) {
2657 if( Y >= 948.0 ) {
2658 /* Stephenson and Morrison, stated domain is 948 to 1600:
2659 * 25.5(centuries from 1800)^2 - 1.9159(centuries from 1955)^2 */
2660 B = 0.01 * (Y - 2000.0);
2661 ans = (23.58 * B + 100.3) * B + 101.6;
2662 } else {
2663 /* Borkowski, before 948 and between 1600 and 1620 */
2664 B = 0.01 * (Y - 2000.0) + 3.75;
2665 ans = 35.0 * B * B + 40.;
2666 }
2667 *deltat = ans / 86400.0;
2668 return iflag;
2669 }
2670 /* 1620 - today + a few years (tabend):
2671 * Tabulated values of deltaT from Astronomical Almanac
2672 * (AA 1997etc., pp. K8-K9) and from IERS
2673 * (http://maia.usno.navy.mil/ser7/deltat.data).
2674 */
2675 if (Y >= TABSTART) {
2676 *deltat = deltat_aa(tjd, tid_acc);
2677 return iflag;
2678 }
2679 #ifdef TRACE
2680 swi_open_trace(NULL);
2681 if (swi_trace_count < TRACE_COUNT_MAX) {
2682 if (swi_fp_trace_c != NULL) {
2683 fputs("\n/*SWE_DELTAT*/\n", swi_fp_trace_c);
2684 fprintf(swi_fp_trace_c, " tjd = %.9f;", tjd);
2685 fprintf(swi_fp_trace_c, " iflag = %d;", iflag);
2686 fprintf(swi_fp_trace_c, " t = swe_deltat_ex(tjd, iflag, NULL);\n");
2687 fputs(" printf(\"swe_deltat: %f\\t%f\\t\\n\", ", swi_fp_trace_c);
2688 fputs("tjd, t);\n", swi_fp_trace_c);
2689 fflush(swi_fp_trace_c);
2690 }
2691 if (swi_fp_trace_out != NULL) {
2692 fprintf(swi_fp_trace_out, "swe_deltat: %f\t%f\t\n", tjd, ans);
2693 fflush(swi_fp_trace_out);
2694 }
2695 }
2696 #endif
2697 *deltat = ans / 86400.0;
2698 return iflag;
2699 }
2700
swe_deltat_ex(double tjd,int32 iflag,char * serr)2701 double CALL_CONV swe_deltat_ex(double tjd, int32 iflag, char *serr)
2702 {
2703 double deltat;
2704 if (swed.delta_t_userdef_is_set)
2705 return swed.delta_t_userdef;
2706 if (serr != NULL)
2707 *serr = '\0';
2708 calc_deltat(tjd, iflag, &deltat, serr);
2709 return deltat;
2710 }
2711
swe_deltat(double tjd)2712 double CALL_CONV swe_deltat(double tjd)
2713 {
2714 int32 iflag = swi_guess_ephe_flag();
2715 return swe_deltat_ex(tjd, iflag, NULL); /* with default tidal acceleration/default ephemeris */
2716 }
2717
2718 /* The tabulated values of deltaT, in hundredths of a second,
2719 * were taken from The Astronomical Almanac 1997etc., pp. K8-K9.
2720 * Some more recent values are taken from IERS
2721 * http://maia.usno.navy.mil/ser7/deltat.data .
2722 * Bessel's interpolation formula is implemented to obtain fourth
2723 * order interpolated values at intermediate times.
2724 * The values are adjusted depending on the ephemeris used
2725 * and its inherent value of secular tidal acceleration ndot.
2726 * Note by Dieter Jan. 2017:
2727 * Bessel interpolation assumes equidistant sampling points. However the
2728 * sampling points are not equidistant, because they are for first January of
2729 * every year and years can have either 365 or 366 days. The interpolation uses
2730 * a step width of 365.25 days. As a consequence, in three out of four years
2731 * the interpolation does not reproduce the exact values of the sampling points
2732 * on the days they refer to. */
deltat_aa(double tjd,double tid_acc)2733 static double deltat_aa(double tjd, double tid_acc)
2734 {
2735 double ans = 0, ans2 = 0, ans3;
2736 double p, B, B2, Y, dd;
2737 double d[6];
2738 int i, iy, k;
2739 /* read additional values from swedelta.txt */
2740 int tabsiz = init_dt();
2741 int tabend = TABSTART + tabsiz - 1;
2742 int deltat_model = swed.astro_models[SE_MODEL_DELTAT];
2743 if (deltat_model == 0) deltat_model = SEMOD_DELTAT_DEFAULT;
2744 Y = 2000.0 + (tjd - 2451544.5)/365.25;
2745 if (Y <= tabend) {
2746 /* Index into the table.
2747 */
2748 p = floor(Y);
2749 iy = (int) (p - TABSTART);
2750 /* Zeroth order estimate is value at start of year */
2751 ans = dt[iy];
2752 k = iy + 1;
2753 if( k >= tabsiz )
2754 goto done; /* No data, can't go on. */
2755 /* The fraction of tabulation interval */
2756 p = Y - p;
2757 /* First order interpolated value */
2758 ans += p*(dt[k] - dt[iy]);
2759 if( (iy-1 < 0) || (iy+2 >= tabsiz) )
2760 goto done; /* can't do second differences */
2761 /* Make table of first differences */
2762 k = iy - 2;
2763 for( i=0; i<5; i++ ) {
2764 if( (k < 0) || (k+1 >= tabsiz) )
2765 d[i] = 0;
2766 else
2767 d[i] = dt[k+1] - dt[k];
2768 k += 1;
2769 }
2770 /* Compute second differences */
2771 for( i=0; i<4; i++ )
2772 d[i] = d[i+1] - d[i];
2773 B = 0.25*p*(p-1.0);
2774 ans += B*(d[1] + d[2]);
2775 #if DEMO
2776 printf( "B %.4lf, ans %.4lf\n", B, ans );
2777 #endif
2778 if( iy+2 >= tabsiz )
2779 goto done;
2780 /* Compute third differences */
2781 for( i=0; i<3; i++ )
2782 d[i] = d[i+1] - d[i];
2783 B = 2.0*B/3.0;
2784 ans += (p-0.5)*B*d[1];
2785 #if DEMO
2786 printf( "B %.4lf, ans %.4lf\n", B*(p-0.5), ans );
2787 #endif
2788 if( (iy-2 < 0) || (iy+3 > tabsiz) )
2789 goto done;
2790 /* Compute fourth differences */
2791 for( i=0; i<2; i++ )
2792 d[i] = d[i+1] - d[i];
2793 B = 0.125*B*(p+1.0)*(p-2.0);
2794 ans += B*(d[0] + d[1]);
2795 #if DEMO
2796 printf( "B %.4lf, ans %.4lf\n", B, ans );
2797 #endif
2798 done:
2799 ans = adjust_for_tidacc(ans, Y, tid_acc, SE_TIDAL_26, FALSE);
2800 return ans / 86400.0;
2801 }
2802 /* today - future:
2803 * 3rd degree polynomial based on data given by
2804 * Stephenson/Morrison/Hohenkerk 2016 here:
2805 * http://astro.ukho.gov.uk/nao/lvm/
2806 */
2807 if (deltat_model == SEMOD_DELTAT_STEPHENSON_ETC_2016) {
2808 B = (Y - 2000);
2809 if (Y < 2500) {
2810 ans = B * B * B * 121.0 / 30000000.0 + B * B / 1250.0 + B * 521.0 / 3000.0 + 64.0;
2811 /* for slow transition from tablulated data */
2812 B2 = (tabend - 2000);
2813 ans2 = B2 * B2 * B2 * 121.0 / 30000000.0 + B2 * B2 / 1250.0 + B2 * 521.0 / 3000.0 + 64.0;
2814 /* we use a parable after 2500 */
2815 } else {
2816 B = 0.01 * (Y - 2000);
2817 ans = B * B * 32.5 + 42.5;
2818 }
2819 /*
2820 * Formula Stephenson (1997; p. 507),
2821 * with modification to avoid jump at end of AA table,
2822 * similar to what Meeus 1998 had suggested.
2823 * Slow transition within 100 years.
2824 */
2825 } else {
2826 B = 0.01 * (Y - 1820);
2827 ans = -20 + 31 * B * B;
2828 /* for slow transition from tablulated data */
2829 B2 = 0.01 * (tabend - 1820);
2830 ans2 = -20 + 31 * B2 * B2;
2831 }
2832 /* slow transition from tabulated values to Stephenson formula: */
2833 if (Y <= tabend+100) {
2834 ans3 = dt[tabsiz-1];
2835 dd = (ans2 - ans3);
2836 ans += dd * (Y - (tabend + 100)) * 0.01;
2837 }
2838 return ans / 86400.0;
2839 }
2840
deltat_longterm_morrison_stephenson(double tjd)2841 static double deltat_longterm_morrison_stephenson(double tjd)
2842 {
2843 double Ygreg = 2000.0 + (tjd - J2000)/365.2425;
2844 double u = (Ygreg - 1820) / 100.0;
2845 return (-20 + 32 * u * u);
2846 }
2847
deltat_stephenson_morrison_1997_1600(double tjd,double tid_acc)2848 static double deltat_stephenson_morrison_1997_1600(double tjd, double tid_acc)
2849 {
2850 double ans = 0, ans2, ans3;
2851 double p, B, Y, dd;
2852 int iy;
2853 Y = 2000.0 + (tjd - J2000)/365.25;
2854 /* before -500:
2855 * formula by Stephenson (1997; p. 508) but adjusted to fit the starting
2856 * point of table dt97 (Stephenson 1997). */
2857 if( Y < TAB97_START ) {
2858 B = (Y - 1735) * 0.01;
2859 ans = -20 + 35 * B * B;
2860 ans = adjust_for_tidacc(ans, Y, tid_acc, SE_TIDAL_26, FALSE);
2861 /* transition from formula to table over 100 years */
2862 if (Y >= TAB97_START - 100) {
2863 /* starting value of table dt97: */
2864 ans2 = adjust_for_tidacc(dt97[0], TAB97_START, tid_acc, SE_TIDAL_26, FALSE);
2865 /* value of formula at epoch TAB97_START */
2866 B = (TAB97_START - 1735) * 0.01;
2867 ans3 = -20 + 35 * B * B;
2868 ans3 = adjust_for_tidacc(ans3, Y, tid_acc, SE_TIDAL_26, FALSE);
2869 dd = ans3 - ans2;
2870 B = (Y - (TAB97_START - 100)) * 0.01;
2871 /* fit to starting point of table dt97. */
2872 ans = ans - dd * B;
2873 }
2874 }
2875 /* between -500 and 1600:
2876 * linear interpolation between values of table dt97 (Stephenson 1997) */
2877 if (Y >= TAB97_START && Y < TAB2_END) {
2878 p = floor(Y);
2879 iy = (int) ((p - TAB97_START) / 50.0);
2880 dd = (Y - (TAB97_START + 50 * iy)) / 50.0;
2881 ans = dt97[iy] + (dt97[iy+1] - dt97[iy]) * dd;
2882 /* correction for tidal acceleration used by our ephemeris */
2883 ans = adjust_for_tidacc(ans, Y, tid_acc, SE_TIDAL_26, FALSE);
2884 }
2885 ans /= 86400.0;
2886 return ans;
2887 }
2888
2889 /* Stephenson & Morrison (2004) */
deltat_stephenson_morrison_2004_1600(double tjd,double tid_acc)2890 static double deltat_stephenson_morrison_2004_1600(double tjd, double tid_acc)
2891 {
2892 double ans = 0, ans2, ans3;
2893 double p, B, dd;
2894 double tjd0;
2895 int iy;
2896 double Y = 2000.0 + (tjd - J2000)/365.2425;
2897 /* double Y = 2000.0 + (tjd - J2000)/365.25;*/
2898 /* before -1000:
2899 * formula by Stephenson & Morrison (2004; p. 335) but adjusted to fit the
2900 * starting point of table dt2. */
2901 if( Y < TAB2_START ) { // before -1000
2902 ans = deltat_longterm_morrison_stephenson(tjd);
2903 ans = adjust_for_tidacc(ans, Y, tid_acc, SE_TIDAL_26, FALSE);
2904 /* transition from formula to table over 100 years */
2905 if (Y >= TAB2_START - 100) {
2906 /* starting value of table dt2: */
2907 ans2 = adjust_for_tidacc(dt2[0], TAB2_START, tid_acc, SE_TIDAL_26, FALSE);
2908 /* value of formula at epoch TAB2_START */
2909 /* B = (TAB2_START - LTERM_EQUATION_YSTART) * 0.01;
2910 ans3 = -20 + LTERM_EQUATION_COEFF * B * B;*/
2911 tjd0 = (TAB2_START - 2000) * 365.2425 + J2000;
2912 ans3 = deltat_longterm_morrison_stephenson(tjd0);
2913 ans3 = adjust_for_tidacc(ans3, Y, tid_acc, SE_TIDAL_26, FALSE);
2914 dd = ans3 - ans2;
2915 B = (Y - (TAB2_START - 100)) * 0.01;
2916 /* fit to starting point of table dt2. */
2917 ans = ans - dd * B;
2918 }
2919 }
2920 /* between -1000 and 1600:
2921 * linear interpolation between values of table dt2 (Stephenson & Morrison 2004) */
2922 if (Y >= TAB2_START && Y < TAB2_END) {
2923 double Yjul = 2000 + (tjd - 2451557.5) / 365.25;
2924 p = floor(Yjul);
2925 iy = (int) ((p - TAB2_START) / TAB2_STEP);
2926 dd = (Yjul - (TAB2_START + TAB2_STEP * iy)) / TAB2_STEP;
2927 ans = dt2[iy] + (dt2[iy+1] - dt2[iy]) * dd;
2928 /* correction for tidal acceleration used by our ephemeris */
2929 ans = adjust_for_tidacc(ans, Y, tid_acc, SE_TIDAL_26, FALSE);
2930 }
2931 ans /= 86400.0;
2932 return ans;
2933 }
2934
2935 /*
2936 * These coefficients represent the spline approximation discussed in the
2937 * paper "Measurement of the Earth's Rotation: 720 BC to AD 2015",
2938 * Stephenson, F.R., Morrison, L.V., and Hohenkerk, C.Y., published by
2939 * Royal Society Proceedings A and available from their website at
2940 * http://rspa.royalsocietypublishing.org/lookup/doi/10.1098/rspa.2016.0404.
2941 * Year numbers have been replaced by Julian day numbers by D. Koch.
2942 */
2943 #define NDTCF16 54
2944 double dtcf16[NDTCF16][6] =
2945 {
2946 /*00*/ {1458085.5, 1867156.5, 20550.593,-21268.478, 11863.418, -4541.129}, /* ybeg=-720, yend= 400 */
2947 /*01*/ {1867156.5, 2086302.5, 6604.404, -5981.266, -505.093, 1349.609}, /* ybeg= 400, yend=1000 */
2948 /*02*/ {2086302.5, 2268923.5, 1467.654, -2452.187, 2460.927, -1183.759}, /* ybeg=1000, yend=1500 */
2949 /*03*/ {2268923.5, 2305447.5, 292.635, -216.322, -43.614, 56.681}, /* ybeg=1500, yend=1600 */
2950 /*04*/ {2305447.5, 2323710.5, 89.380, -66.754, 31.607, -10.497}, /* ybeg=1600, yend=1650 */
2951 /*05*/ {2323710.5, 2349276.5, 43.736, -49.043, 0.227, 15.811}, /* ybeg=1650, yend=1720 */
2952 /*06*/ {2349276.5, 2378496.5, 10.730, -1.321, 62.250, -52.946}, /* ybeg=1720, yend=1800 */
2953 /*07*/ {2378496.5, 2382148.5, 18.714, -4.457, -1.509, 2.507}, /* ybeg=1800, yend=1810 */
2954 /*08*/ {2382148.5, 2385800.5, 15.255, 0.046, 6.012, -4.634}, /* ybeg=1810, yend=1820 */
2955 /*09*/ {2385800.5, 2389453.5, 16.679, -1.831, -7.889, 3.799}, /* ybeg=1820, yend=1830 */
2956 /*10*/ {2389453.5, 2393105.5, 10.758, -6.211, 3.509, -0.388}, /* ybeg=1830, yend=1840 */
2957 /*11*/ {2393105.5, 2396758.5, 7.668, -0.357, 2.345, -0.338}, /* ybeg=1840, yend=1850 */
2958 /*12*/ {2396758.5, 2398584.5, 9.317, 1.659, 0.332, -0.932}, /* ybeg=1850, yend=1855 */
2959 /*13*/ {2398584.5, 2400410.5, 10.376, -0.472, -2.463, 1.596}, /* ybeg=1855, yend=1860 */
2960 /*14*/ {2400410.5, 2402237.5, 9.038, -0.610, 2.325, -2.497}, /* ybeg=1860, yend=1865 */
2961 /*15*/ {2402237.5, 2404063.5, 8.256, -3.450, -5.166, 2.729}, /* ybeg=1865, yend=1870 */
2962 /*16*/ {2404063.5, 2405889.5, 2.369, -5.596, 3.020, -0.919}, /* ybeg=1870, yend=1875 */
2963 /*17*/ {2405889.5, 2407715.5, -1.126, -2.312, 0.264, -0.037}, /* ybeg=1875, yend=1880 */
2964 /*18*/ {2407715.5, 2409542.5, -3.211, -1.894, 0.154, 0.562}, /* ybeg=1880, yend=1885 */
2965 /*19*/ {2409542.5, 2411368.5, -4.388, 0.101, 1.841, -1.438}, /* ybeg=1885, yend=1890 */
2966 /*20*/ {2411368.5, 2413194.5, -3.884, -0.531, -2.473, 1.870}, /* ybeg=1890, yend=1895 */
2967 /*21*/ {2413194.5, 2415020.5, -5.017, 0.134, 3.138, -0.232}, /* ybeg=1895, yend=1900 */
2968 /*22*/ {2415020.5, 2416846.5, -1.977, 5.715, 2.443, -1.257}, /* ybeg=1900, yend=1905 */
2969 /*23*/ {2416846.5, 2418672.5, 4.923, 6.828, -1.329, 0.720}, /* ybeg=1905, yend=1910 */
2970 /*24*/ {2418672.5, 2420498.5, 11.142, 6.330, 0.831, -0.825}, /* ybeg=1910, yend=1915 */
2971 /*25*/ {2420498.5, 2422324.5, 17.479, 5.518, -1.643, 0.262}, /* ybeg=1915, yend=1920 */
2972 /*26*/ {2422324.5, 2424151.5, 21.617, 3.020, -0.856, 0.008}, /* ybeg=1920, yend=1925 */
2973 /*27*/ {2424151.5, 2425977.5, 23.789, 1.333, -0.831, 0.127}, /* ybeg=1925, yend=1930 */
2974 /*28*/ {2425977.5, 2427803.5, 24.418, 0.052, -0.449, 0.142}, /* ybeg=1930, yend=1935 */
2975 /*29*/ {2427803.5, 2429629.5, 24.164, -0.419, -0.022, 0.702}, /* ybeg=1935, yend=1940 */
2976 /*30*/ {2429629.5, 2431456.5, 24.426, 1.645, 2.086, -1.106}, /* ybeg=1940, yend=1945 */
2977 /*31*/ {2431456.5, 2433282.5, 27.050, 2.499, -1.232, 0.614}, /* ybeg=1945, yend=1950 */
2978 /*32*/ {2433282.5, 2434378.5, 28.932, 1.127, 0.220, -0.277}, /* ybeg=1950, yend=1953 */
2979 /*33*/ {2434378.5, 2435473.5, 30.002, 0.737, -0.610, 0.631}, /* ybeg=1953, yend=1956 */
2980 /*34*/ {2435473.5, 2436569.5, 30.760, 1.409, 1.282, -0.799}, /* ybeg=1956, yend=1959 */
2981 /*35*/ {2436569.5, 2437665.5, 32.652, 1.577, -1.115, 0.507}, /* ybeg=1959, yend=1962 */
2982 /*36*/ {2437665.5, 2438761.5, 33.621, 0.868, 0.406, 0.199}, /* ybeg=1962, yend=1965 */
2983 /*37*/ {2438761.5, 2439856.5, 35.093, 2.275, 1.002, -0.414}, /* ybeg=1965, yend=1968 */
2984 /*38*/ {2439856.5, 2440952.5, 37.956, 3.035, -0.242, 0.202}, /* ybeg=1968, yend=1971 */
2985 /*39*/ {2440952.5, 2442048.5, 40.951, 3.157, 0.364, -0.229}, /* ybeg=1971, yend=1974 */
2986 /*40*/ {2442048.5, 2443144.5, 44.244, 3.198, -0.323, 0.172}, /* ybeg=1974, yend=1977 */
2987 /*41*/ {2443144.5, 2444239.5, 47.291, 3.069, 0.193, -0.192}, /* ybeg=1977, yend=1980 */
2988 /*42*/ {2444239.5, 2445335.5, 50.361, 2.878, -0.384, 0.081}, /* ybeg=1980, yend=1983 */
2989 /*43*/ {2445335.5, 2446431.5, 52.936, 2.354, -0.140, -0.166}, /* ybeg=1983, yend=1986 */
2990 /*44*/ {2446431.5, 2447527.5, 54.984, 1.577, -0.637, 0.448}, /* ybeg=1986, yend=1989 */
2991 /*45*/ {2447527.5, 2448622.5, 56.373, 1.649, 0.709, -0.277}, /* ybeg=1989, yend=1992 */
2992 /*46*/ {2448622.5, 2449718.5, 58.453, 2.235, -0.122, 0.111}, /* ybeg=1992, yend=1995 */
2993 /*47*/ {2449718.5, 2450814.5, 60.677, 2.324, 0.212, -0.315}, /* ybeg=1995, yend=1998 */
2994 /*48*/ {2450814.5, 2451910.5, 62.899, 1.804, -0.732, 0.112}, /* ybeg=1998, yend=2001 */
2995 /*49*/ {2451910.5, 2453005.5, 64.082, 0.675, -0.396, 0.193}, /* ybeg=2001, yend=2004 */
2996 /*50*/ {2453005.5, 2454101.5, 64.555, 0.463, 0.184, -0.008}, /* ybeg=2004, yend=2007 */
2997 /*51*/ {2454101.5, 2455197.5, 65.194, 0.809, 0.161, -0.101}, /* ybeg=2007, yend=2010 */
2998 /*52*/ {2455197.5, 2456293.5, 66.063, 0.828, -0.142, 0.168}, /* ybeg=2010, yend=2013 */
2999 /*53*/ {2456293.5, 2457388.5, 66.917, 1.046, 0.360, -0.282}, /* ybeg=2013, yend=2016 */
3000 };
deltat_stephenson_etc_2016(double tjd,double tid_acc)3001 static double deltat_stephenson_etc_2016(double tjd, double tid_acc)
3002 {
3003 double t, dt, Ygreg;
3004 int i, irec = -1;
3005 Ygreg = 2000.0 + (tjd - J2000)/365.2425;
3006 // after the year -720 get value from spline curve
3007 for (i = 0; i < NDTCF16; i++) {
3008 if (tjd < dtcf16[i][0]) break;
3009 if (tjd < dtcf16[i][1]) {
3010 irec = i;
3011 break;
3012 }
3013 }
3014 if (irec >= 0) {
3015 t = (tjd - dtcf16[irec][0]) / (dtcf16[irec][1] - dtcf16[irec][0]);
3016 dt = dtcf16[irec][2] + dtcf16[irec][3] * t + dtcf16[irec][4] * t * t + dtcf16[irec][5] * t * t * t;
3017 // for earlier epochs, use long term parabola
3018 } else if (Ygreg < -720) {
3019 t = (Ygreg - 1825) / 100.0;
3020 dt = -320 + 32.5 * t * t;
3021 dt -= 179.7337208; // to make curve continous on 1 Jan -720 (D. Koch)
3022 // future
3023 } else {
3024 t = (Ygreg - 1825) / 100.0;
3025 dt = -320 + 32.5 * t * t;
3026 dt += 269.4790417; // to make curve continous on 1 Jan 2016 (D. Koch)
3027 }
3028 /* The parameter adjust_after_1955 must be TRUE here, because the
3029 * Stephenson 2016 curve is based on occultation data alone,
3030 * not on IERS data.
3031 * Note, however, the current function deltat_stephenson_etc_2016()
3032 * is called only for dates before 1 Jan 1955. */
3033 dt = adjust_for_tidacc(dt, Ygreg, tid_acc, SE_TIDAL_STEPHENSON_2016, TRUE);
3034 dt /= 86400.0;
3035 return dt;
3036 }
3037
deltat_espenak_meeus_1620(double tjd,double tid_acc)3038 static double deltat_espenak_meeus_1620(double tjd, double tid_acc)
3039 {
3040 double ans = 0;
3041 double Ygreg;
3042 double u;
3043 /* double Y = 2000.0 + (tjd - J2000)/365.25;*/
3044 Ygreg = 2000.0 + (tjd - J2000)/365.2425;
3045 if (Ygreg < -500) {
3046 ans = deltat_longterm_morrison_stephenson(tjd);
3047 } else if (Ygreg < 500) {
3048 u = Ygreg / 100.0;
3049 ans = (((((0.0090316521 * u + 0.022174192) * u - 0.1798452) * u - 5.952053) * u+ 33.78311) * u - 1014.41) * u + 10583.6;
3050 } else if (Ygreg < 1600) {
3051 u = (Ygreg - 1000) / 100.0;
3052 ans = (((((0.0083572073 * u - 0.005050998) * u - 0.8503463) * u + 0.319781) * u + 71.23472) * u - 556.01) * u + 1574.2;
3053 } else if (Ygreg < 1700) {
3054 u = Ygreg - 1600;
3055 ans = 120 - 0.9808 * u - 0.01532 * u * u + u * u * u / 7129.0;
3056 } else if (Ygreg < 1800) {
3057 u = Ygreg - 1700;
3058 ans = (((-u / 1174000.0 + 0.00013336) * u - 0.0059285) * u + 0.1603) * u + 8.83;
3059 } else if (Ygreg < 1860) {
3060 u = Ygreg - 1800;
3061 ans = ((((((0.000000000875 * u - 0.0000001699) * u + 0.0000121272) * u - 0.00037436) * u + 0.0041116) * u + 0.0068612) * u - 0.332447) * u + 13.72;
3062 } else if (Ygreg < 1900) {
3063 u = Ygreg - 1860;
3064 ans = ((((u / 233174.0 - 0.0004473624) * u + 0.01680668) * u - 0.251754) * u + 0.5737) * u + 7.62;
3065 } else if (Ygreg < 1920) {
3066 u = Ygreg - 1900;
3067 ans = (((-0.000197 * u + 0.0061966) * u - 0.0598939) * u + 1.494119) * u -2.79;
3068 } else if (Ygreg < 1941) {
3069 u = Ygreg - 1920;
3070 ans = 21.20 + 0.84493 * u - 0.076100 * u * u + 0.0020936 * u * u * u;
3071 } else if (Ygreg < 1961) {
3072 u = Ygreg - 1950;
3073 ans = 29.07 + 0.407 * u - u * u / 233.0 + u * u * u / 2547.0;
3074 } else if (Ygreg < 1986) {
3075 u = Ygreg - 1975;
3076 ans = 45.45 + 1.067 * u - u * u / 260.0 - u * u * u / 718.0;
3077 } else if (Ygreg < 2005) {
3078 u = Ygreg - 2000;
3079 ans = ((((0.00002373599 * u + 0.000651814) * u + 0.0017275) * u - 0.060374) * u + 0.3345) * u + 63.86;
3080 }
3081 ans = adjust_for_tidacc(ans, Ygreg, tid_acc, SE_TIDAL_26, FALSE);
3082 ans /= 86400.0;
3083 return ans;
3084 }
3085
3086 /* Read delta t values from external file.
3087 * record structure: year(whitespace)delta_t in 0.01 sec.
3088 */
init_dt(void)3089 static int init_dt(void)
3090 {
3091 FILE *fp;
3092 int year;
3093 int tab_index;
3094 int tabsiz;
3095 int i;
3096 char s[AS_MAXCH];
3097 char *sp;
3098 if (!swed.init_dt_done) {
3099 swed.init_dt_done = TRUE;
3100 /* no error message if file is missing */
3101 if ((fp = swi_fopen(-1, "swe_deltat.txt", swed.ephepath, NULL)) == NULL
3102 && (fp = swi_fopen(-1, "sedeltat.txt", swed.ephepath, NULL)) == NULL)
3103 return TABSIZ;
3104 while(fgets(s, AS_MAXCH, fp) != NULL) {
3105 sp = s;
3106 while (strchr(" \t", *sp) != NULL && *sp != '\0')
3107 sp++; /* was *sp++ fixed by Alois 2-jul-2003 */
3108 if (*sp == '#' || *sp == '\n')
3109 continue;
3110 year = atoi(s);
3111 tab_index = year - TABSTART;
3112 /* table space is limited. no error msg, if exceeded */
3113 if (tab_index >= TABSIZ_SPACE)
3114 continue;
3115 sp += 4;
3116 while (strchr(" \t", *sp) != NULL && *sp != '\0')
3117 sp++; /* was *sp++ fixed by Alois 2-jul-2003 */
3118 /*dt[tab_index] = (short) (atof(sp) * 100 + 0.5);*/
3119 dt[tab_index] = atof(sp);
3120 }
3121 fclose(fp);
3122 }
3123 /* find table size */
3124 tabsiz = 2001 - TABSTART + 1;
3125 for (i = tabsiz - 1; i < TABSIZ_SPACE; i++) {
3126 if (dt[i] == 0)
3127 break;
3128 else
3129 tabsiz++;
3130 }
3131 tabsiz--;
3132 return tabsiz;
3133 }
3134
3135 /* Astronomical Almanac table is corrected by adding the expression
3136 * -0.000091 (ndot + 26)(year-1955)^2 seconds
3137 * to entries prior to 1955 (AA page K8), where ndot is the secular
3138 * tidal term in the mean motion of the Moon.
3139 *
3140 * Entries after 1955 are referred to atomic time standards and
3141 * are not affected by errors in Lunar or planetary theory.
3142 */
adjust_for_tidacc(double ans,double Y,double tid_acc,double tid_acc0,AS_BOOL adjust_after_1955)3143 static double adjust_for_tidacc(double ans, double Y, double tid_acc, double tid_acc0, AS_BOOL adjust_after_1955)
3144 {
3145 double B;
3146 if( Y < 1955.0 || adjust_after_1955) {
3147 B = (Y - 1955.0);
3148 ans += -0.000091 * (tid_acc - tid_acc0) * B * B;
3149 }
3150 return ans;
3151 }
3152
3153 /* returns tidal acceleration used in swe_deltat() and swe_deltat_ex() */
swe_get_tid_acc()3154 double CALL_CONV swe_get_tid_acc()
3155 {
3156 return swed.tid_acc;
3157 }
3158
3159 /* function sets tidal acceleration of the Moon.
3160 * t_acc can be either
3161 * - the value of the tidal acceleration in arcsec/cty^2
3162 * of the Moon will be set consistent with that ephemeris.
3163 * - SE_TIDAL_AUTOMATIC,
3164 */
swe_set_tid_acc(double t_acc)3165 void CALL_CONV swe_set_tid_acc(double t_acc)
3166 {
3167 if (t_acc == SE_TIDAL_AUTOMATIC) {
3168 swed.tid_acc = SE_TIDAL_DEFAULT;
3169 swed.is_tid_acc_manual = FALSE;
3170 return;
3171 }
3172 swed.tid_acc = t_acc;
3173 swed.is_tid_acc_manual = TRUE;
3174 }
3175
swe_set_delta_t_userdef(double dt)3176 void CALL_CONV swe_set_delta_t_userdef(double dt)
3177 {
3178 if (dt == SE_DELTAT_AUTOMATIC) {
3179 swed.delta_t_userdef_is_set = FALSE;
3180 } else {
3181 swed.delta_t_userdef_is_set = TRUE;
3182 swed.delta_t_userdef = dt;
3183 }
3184 }
3185
swi_guess_ephe_flag()3186 int32 swi_guess_ephe_flag()
3187 {
3188 int32 iflag = SEFLG_SWIEPH;
3189 /* if jpl file is open, assume SEFLG_JPLEPH */
3190 if (swed.jpl_file_is_open) {
3191 iflag = SEFLG_JPLEPH;
3192 } else {
3193 iflag = SEFLG_SWIEPH;
3194 }
3195 return iflag;
3196 }
3197
swi_get_tid_acc(double tjd_ut,int32 iflag,int32 denum,int32 * denumret,double * tid_acc,char * serr)3198 int32 swi_get_tid_acc(double tjd_ut, int32 iflag, int32 denum, int32 *denumret, double *tid_acc, char *serr)
3199 {
3200 iflag &= SEFLG_EPHMASK;
3201 if (swed.is_tid_acc_manual) {
3202 *tid_acc = swed.tid_acc;
3203 return iflag;
3204 }
3205 if (denum == 0) {
3206 if (iflag & SEFLG_MOSEPH) {
3207 *tid_acc = SE_TIDAL_DE404;
3208 *denumret = 404;
3209 return iflag;
3210 }
3211 if (iflag & SEFLG_JPLEPH) {
3212 if (swed.jpl_file_is_open) {
3213 denum = swed.jpldenum;
3214 }
3215 }
3216 /* SEFLG_SWIEPH wanted or SEFLG_JPLEPH failed: */
3217 if (iflag & SEFLG_SWIEPH) {
3218 if (swed.fidat[SEI_FILE_MOON].fptr != NULL) {
3219 denum = swed.fidat[SEI_FILE_MOON].sweph_denum;
3220 }
3221 }
3222 }
3223 switch(denum) {
3224 case 200: *tid_acc = SE_TIDAL_DE200; break;
3225 case 403: *tid_acc = SE_TIDAL_DE403; break;
3226 case 404: *tid_acc = SE_TIDAL_DE404; break;
3227 case 405: *tid_acc = SE_TIDAL_DE405; break;
3228 case 406: *tid_acc = SE_TIDAL_DE406; break;
3229 case 421: *tid_acc = SE_TIDAL_DE421; break;
3230 case 422: *tid_acc = SE_TIDAL_DE422; break;
3231 case 430: *tid_acc = SE_TIDAL_DE430; break;
3232 case 431: *tid_acc = SE_TIDAL_DE431; break;
3233 case 440: *tid_acc = SE_TIDAL_DE441; break;
3234 case 441: *tid_acc = SE_TIDAL_DE441; break;
3235 default: denum = SE_DE_NUMBER; *tid_acc = SE_TIDAL_DEFAULT; break;
3236 }
3237 *denumret = denum;
3238 iflag &= SEFLG_EPHMASK;
3239 return iflag;
3240 }
3241
swi_set_tid_acc(double tjd_ut,int32 iflag,int32 denum,char * serr)3242 int32 swi_set_tid_acc(double tjd_ut, int32 iflag, int32 denum, char *serr)
3243 {
3244 int32 retc = iflag;
3245 int32 denumret;
3246 /* manual tid_acc overrides automatic tid_acc */
3247 if (swed.is_tid_acc_manual)
3248 return retc;
3249 retc = swi_get_tid_acc(tjd_ut, iflag, denum, &denumret, &(swed.tid_acc), serr);
3250 #if TRACE
3251 swi_open_trace(NULL);
3252 if (swi_trace_count < TRACE_COUNT_MAX) {
3253 if (swi_fp_trace_c != NULL) {
3254 fputs("\n/*SWE_SET_TID_ACC*/\n", swi_fp_trace_c);
3255 fprintf(swi_fp_trace_c, " t = %.9f;\n", swed.tid_acc);
3256 fprintf(swi_fp_trace_c, " swe_set_tid_acc(t);\n");
3257 fputs(" printf(\"swe_set_tid_acc: %f\\t\\n\", ", swi_fp_trace_c);
3258 fputs("t);\n", swi_fp_trace_c);
3259 fflush(swi_fp_trace_c);
3260 }
3261 if (swi_fp_trace_out != NULL) {
3262 fprintf(swi_fp_trace_out, "swe_set_tid_acc: %f\t\n", swed.tid_acc);
3263 fflush(swi_fp_trace_out);
3264 }
3265 }
3266 #endif
3267 return retc;
3268 }
3269
3270 /*
3271 * The time range of DE431 requires a new calculation of sidereal time that
3272 * gives sensible results for the remote past and future.
3273 * The algorithm is based on the formula of the mean earth by Simon & alii,
3274 * "Precession formulae and mean elements for the Moon and the Planets",
3275 * A&A 282 (1994), p. 675/678.
3276 * The longitude of the mean earth relative to the mean equinox J2000
3277 * is calculated and then precessed to the equinox of date, using the
3278 * default precession model of the Swiss Ephmeris. Afte that,
3279 * sidereal time is derived.
3280 * The algoritm provides exact agreement for epoch 1 Jan. 2003 with the
3281 * definition of sidereal time as given in the IERS Convention 2010.
3282 */
3283 /*#define SIDT_LTERM TRUE
3284 #if SIDT_LTERM*/
sidtime_long_term(double tjd_ut,double eps,double nut)3285 static double sidtime_long_term(double tjd_ut, double eps, double nut)
3286 {
3287 double tsid = 0, tjd_et;
3288 double dlon, xs[6], xobl[6], dhour, nutlo[2];
3289 double dlt = AUNIT / CLIGHT / 86400.0;
3290 double t, t2, t3;
3291 tjd_et = tjd_ut + swe_deltat_ex(tjd_ut, -1, NULL);
3292 t = (tjd_et - J2000) / 365250.0;
3293 t2 = t * t; t3 = t * t2;
3294 /* mean longitude of earth J2000 */
3295 dlon = 100.46645683 + (1295977422.83429 * t - 2.04411 * t2 - 0.00523 * t3) / 3600.0;
3296 /* light time sun-earth */
3297 dlon = swe_degnorm(dlon - dlt * 360.0 / 365.2425);
3298 xs[0] = dlon * DEGTORAD; xs[1] = 0; xs[2] = 1;
3299 /* to mean equator J2000, cartesian */
3300 xobl[0] = 23.45; xobl[1] = 23.45;
3301 xobl[1] = swi_epsiln(J2000 + swe_deltat_ex(J2000, -1, NULL), 0) * RADTODEG;
3302 swi_polcart(xs, xs);
3303 swi_coortrf(xs, xs, -xobl[1] * DEGTORAD);
3304 /* precess to mean equinox of date */
3305 swi_precess(xs, tjd_et, 0, -1);
3306 /* to mean equinox of date */
3307 xobl[1] = swi_epsiln(tjd_et, 0) * RADTODEG;
3308 swi_nutation(tjd_et, 0, nutlo);
3309 xobl[0] = xobl[1] + nutlo[1] * RADTODEG;
3310 xobl[2] = nutlo[0] * RADTODEG;
3311 swi_coortrf(xs, xs, xobl[1] * DEGTORAD);
3312 swi_cartpol(xs, xs);
3313 xs[0] *= RADTODEG;
3314 dhour = fmod(tjd_ut - 0.5, 1) * 360;
3315 /* mean to true (if nut != 0) */
3316 if (eps == 0)
3317 xs[0] += xobl[2] * cos(xobl[0] * DEGTORAD);
3318 else
3319 xs[0] += nut * cos(eps * DEGTORAD);
3320 /* add hour */
3321 xs[0] = swe_degnorm(xs[0] + dhour);
3322 tsid = xs[0] / 15;
3323 return tsid;
3324 }
3325 /*#endif*/
3326
3327 /* Apparent Sidereal Time at Greenwich with equation of the equinoxes
3328 * ERA-based expression for for Greenwich Sidereal Time (GST) based
3329 * on the IAU 2006 precession and IAU 2000A_R06 nutation
3330 * ftp://maia.usno.navy.mil/conv2010/chapter5/tab5.2e.txt
3331 *
3332 * returns sidereal time in hours.
3333 *
3334 * program returns sidereal hours since sidereal midnight
3335 * tjd julian day UT
3336 * eps obliquity of ecliptic, degrees
3337 * nut nutation, degrees
3338 */
3339 /* C'_{s,j})_i C'_{c,j})_i */
3340 #define SIDTNTERM 33
3341 static const double stcf[SIDTNTERM * 2] = {
3342 2640.96,-0.39,
3343 63.52,-0.02,
3344 11.75,0.01,
3345 11.21,0.01,
3346 -4.55,0.00,
3347 2.02,0.00,
3348 1.98,0.00,
3349 -1.72,0.00,
3350 -1.41,-0.01,
3351 -1.26,-0.01,
3352 -0.63,0.00,
3353 -0.63,0.00,
3354 0.46,0.00,
3355 0.45,0.00,
3356 0.36,0.00,
3357 -0.24,-0.12,
3358 0.32,0.00,
3359 0.28,0.00,
3360 0.27,0.00,
3361 0.26,0.00,
3362 -0.21,0.00,
3363 0.19,0.00,
3364 0.18,0.00,
3365 -0.10,0.05,
3366 0.15,0.00,
3367 -0.14,0.00,
3368 0.14,0.00,
3369 -0.14,0.00,
3370 0.14,0.00,
3371 0.13,0.00,
3372 -0.11,0.00,
3373 0.11,0.00,
3374 0.11,0.00,
3375 };
3376 #define SIDTNARG 14
3377 /* l l' F D Om L_Me L_Ve L_E L_Ma L_J L_Sa L_U L_Ne p_A*/
3378 static const int stfarg[SIDTNTERM * SIDTNARG] = {
3379 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3380 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3381 0, 0, 2, -2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3382 0, 0, 2, -2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3383 0, 0, 2, -2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3384 0, 0, 2, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3385 0, 0, 2, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3386 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3387 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3388 0, 1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3389 1, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3390 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3391 0, 1, 2, -2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3392 0, 1, 2, -2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3393 0, 0, 4, -4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3394 0, 0, 1, -1, 1, 0, -8, 12, 0, 0, 0, 0, 0, 0,
3395 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3396 0, 0, 2, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3397 1, 0, 2, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3398 1, 0, 2, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3399 0, 0, 2, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3400 0, 1, -2, 2, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3401 0, 1, -2, 2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3402 0, 0, 0, 0, 0, 0, 8, -13, 0, 0, 0, 0, 0, -1,
3403 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3404 2, 0, -2, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3405 1, 0, 0, -2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3406 0, 1, 2, -2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3407 1, 0, 0, -2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3408 0, 0, 4, -2, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3409 0, 0, 2, -2, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3410 1, 0, -2, 0, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3411 1, 0, -2, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3412 };
sidtime_non_polynomial_part(double tt)3413 static double sidtime_non_polynomial_part(double tt)
3414 {
3415 int i, j;
3416 double delm[SIDTNARG];
3417 double dadd, darg;
3418 /* L Mean anomaly of the Moon.*/
3419 delm[0] = swe_radnorm(2.35555598 + 8328.6914269554 * tt);
3420 /* LSU Mean anomaly of the Sun.*/
3421 delm[1] = swe_radnorm(6.24006013 + 628.301955 * tt);
3422 /* F Mean argument of the latitude of the Moon. */
3423 delm[2] = swe_radnorm(1.627905234 + 8433.466158131 * tt);
3424 /* D Mean elongation of the Moon from the Sun. */
3425 delm[3] = swe_radnorm(5.198466741 + 7771.3771468121 * tt);
3426 /* OM Mean longitude of the ascending node of the Moon. */
3427 delm[4] = swe_radnorm(2.18243920 - 33.757045 * tt);
3428 /* Planetary longitudes, Mercury through Neptune (Souchay et al. 1999).
3429 * LME, LVE, LEA, LMA, LJU, LSA, LUR, LNE */
3430 delm[5] = swe_radnorm(4.402608842 + 2608.7903141574 * tt);
3431 delm[6] = swe_radnorm(3.176146697 + 1021.3285546211 * tt);
3432 delm[7] = swe_radnorm(1.753470314 + 628.3075849991 * tt);
3433 delm[8] = swe_radnorm(6.203480913 + 334.0612426700 * tt);
3434 delm[9] = swe_radnorm(0.599546497 + 52.9690962641 * tt);
3435 delm[10] = swe_radnorm(0.874016757 + 21.3299104960 * tt);
3436 delm[11] = swe_radnorm(5.481293871 + 7.4781598567 * tt);
3437 delm[12] = swe_radnorm(5.321159000 + 3.8127774000 * tt);
3438 /* PA General accumulated precession in longitude. */
3439 delm[13] = (0.02438175 + 0.00000538691 * tt) * tt;
3440 dadd = -0.87 * sin(delm[4]) * tt;
3441 for (i = 0; i < SIDTNTERM; i++) {
3442 darg = 0;
3443 for (j = 0; j < SIDTNARG; j++) {
3444 darg += stfarg[i * SIDTNARG + j] * delm[j];
3445 }
3446 dadd += stcf[i * 2] * sin(darg) + stcf[i * 2 + 1] * cos(darg);
3447 }
3448 dadd /= (3600.0 * 1000000.0);
3449 return dadd;
3450 }
3451
3452 /*
3453 * SEMOD_SIDT_IAU_2006
3454 * N. Capitaine, P.T. Wallace, and J. Chapront, "Expressions for IAU 2000
3455 * precession quantities", 2003, A&A 412, 567-586 (2003), p. 582.
3456 * This is a "short" term model, that can be combined with other models
3457 */
3458 /*#define SIDT_IERS_CONV_2010 TRUE*/
3459 /* sidtime_long_term() is not used between the following two dates */
3460 #define SIDT_LTERM_T0 2396758.5 /* 1 Jan 1850 */
3461 #define SIDT_LTERM_T1 2469807.5 /* 1 Jan 2050 */
3462 #define SIDT_LTERM_OFS0 (0.000378172 / 15.0)
3463 #define SIDT_LTERM_OFS1 (0.001385646 / 15.0)
swe_sidtime0(double tjd,double eps,double nut)3464 double CALL_CONV swe_sidtime0(double tjd, double eps, double nut)
3465 {
3466 double jd0; /* Julian day at midnight Universal Time */
3467 double secs; /* Time of day, UT seconds since UT midnight */
3468 double eqeq, jd, tu, tt, msday, jdrel;
3469 double gmst, dadd;
3470 int prec_model_short = swed.astro_models[SE_MODEL_PREC_SHORTTERM];
3471 int sidt_model = swed.astro_models[SE_MODEL_SIDT];
3472 if (prec_model_short == 0) prec_model_short = SEMOD_PREC_DEFAULT_SHORT;
3473 if (sidt_model == 0) sidt_model = SEMOD_SIDT_DEFAULT;
3474 swi_init_swed_if_start();
3475 if (sidt_model == SEMOD_SIDT_LONGTERM) {
3476 if (tjd <= SIDT_LTERM_T0 || tjd >= SIDT_LTERM_T1) {
3477 gmst = sidtime_long_term(tjd, eps, nut);
3478 if (tjd <= SIDT_LTERM_T0)
3479 gmst -= SIDT_LTERM_OFS0;
3480 else if (tjd >= SIDT_LTERM_T1)
3481 gmst -= SIDT_LTERM_OFS1;
3482 if (gmst >= 24) gmst -= 24;
3483 if (gmst < 0) gmst += 24;
3484 goto sidtime_done;
3485 }
3486 }
3487 /* Julian day at given UT */
3488 jd = tjd;
3489 jd0 = floor(jd);
3490 secs = tjd - jd0;
3491 if( secs < 0.5 ) {
3492 jd0 -= 0.5;
3493 secs += 0.5;
3494 } else {
3495 jd0 += 0.5;
3496 secs -= 0.5;
3497 }
3498 secs *= 86400.0;
3499 tu = (jd0 - J2000)/36525.0; /* UT1 in centuries after J2000 */
3500 if (sidt_model == SEMOD_SIDT_IERS_CONV_2010 || sidt_model == SEMOD_SIDT_LONGTERM) {
3501 /* ERA-based expression for Greenwich Sidereal Time (GST) based
3502 * on the IAU 2006 precession */
3503 jdrel = tjd - J2000;
3504 tt = (tjd + swe_deltat_ex(tjd, -1, NULL) - J2000) / 36525.0;
3505 gmst = swe_degnorm((0.7790572732640 + 1.00273781191135448 * jdrel) * 360);
3506 gmst += (0.014506 + tt * (4612.156534 + tt * (1.3915817 + tt * (-0.00000044 + tt * (-0.000029956 + tt * -0.0000000368))))) / 3600.0;
3507 dadd = sidtime_non_polynomial_part(tt);
3508 gmst = swe_degnorm(gmst + dadd);
3509 /*printf("gmst iers=%f \n", gmst);*/
3510 gmst = gmst / 15.0 * 3600.0;
3511 /* sidt_model == SEMOD_SIDT_IAU_2006, older standards according to precession model */
3512 } else if (sidt_model == SEMOD_SIDT_IAU_2006) {
3513 tt = (jd0 + swe_deltat_ex(jd0, -1, NULL) - J2000)/36525.0; /* TT in centuries after J2000 */
3514 gmst = (((-0.000000002454*tt - 0.00000199708)*tt - 0.0000002926)*tt + 0.092772110)*tt*tt + 307.4771013*(tt-tu) + 8640184.79447825*tu + 24110.5493771;
3515 /* mean solar days per sidereal day at date tu;
3516 * for the derivative of gmst, we can assume UT1 =~ TT */
3517 msday = 1 + ((((-0.000000012270*tt - 0.00000798832)*tt - 0.0000008778)*tt + 0.185544220)*tt + 8640184.79447825)/(86400.*36525.);
3518 gmst += msday * secs;
3519 /* SEMOD_SIDT_IAU_1976 */
3520 } else { /* IAU 1976 formula */
3521 /* Greenwich Mean Sidereal Time at 0h UT of date */
3522 gmst = (( -6.2e-6*tu + 9.3104e-2)*tu + 8640184.812866)*tu + 24110.54841;
3523 /* mean solar days per sidereal day at date tu, = 1.00273790934 in 1986 */
3524 msday = 1.0 + ((-1.86e-5*tu + 0.186208)*tu + 8640184.812866)/(86400.*36525.);
3525 gmst += msday * secs;
3526 }
3527 /* Local apparent sidereal time at given UT at Greenwich */
3528 eqeq = 240.0 * nut * cos(eps * DEGTORAD);
3529 gmst = gmst + eqeq /* + 240.0*tlong */;
3530 /* Sidereal seconds modulo 1 sidereal day */
3531 gmst = gmst - 86400.0 * floor( gmst/86400.0 );
3532 /* return in hours */
3533 gmst /= 3600;
3534 goto sidtime_done;
3535 sidtime_done:
3536 #ifdef TRACE
3537 swi_open_trace(NULL);
3538 if (swi_trace_count < TRACE_COUNT_MAX) {
3539 if (swi_fp_trace_c != NULL) {
3540 fputs("\n/*SWE_SIDTIME0*/\n", swi_fp_trace_c);
3541 fprintf(swi_fp_trace_c, " tjd = %.9f;", tjd);
3542 fprintf(swi_fp_trace_c, " eps = %.9f;", eps);
3543 fprintf(swi_fp_trace_c, " nut = %.9f;\n", nut);
3544 fprintf(swi_fp_trace_c, " t = swe_sidtime0(tjd, eps, nut);\n");
3545 fputs(" printf(\"swe_sidtime0: %f\\tsidt = %f\\teps = %f\\tnut = %f\\t\\n\", ", swi_fp_trace_c);
3546 fputs("tjd, t, eps, nut);\n", swi_fp_trace_c);
3547 fflush(swi_fp_trace_c);
3548 }
3549 if (swi_fp_trace_out != NULL) {
3550 fprintf(swi_fp_trace_out, "swe_sidtime0: %f\tsidt = %f\teps = %f\tnut = %f\t\n", tjd, gmst, eps, nut);
3551 fflush(swi_fp_trace_out);
3552 }
3553 }
3554 #endif
3555 return gmst;
3556 }
3557
swe_set_interpolate_nut(AS_BOOL do_interpolate)3558 void CALL_CONV swe_set_interpolate_nut(AS_BOOL do_interpolate)
3559 {
3560 if (swed.do_interpolate_nut == do_interpolate)
3561 return;
3562 if (do_interpolate)
3563 swed.do_interpolate_nut = TRUE;
3564 else
3565 swed.do_interpolate_nut = FALSE;
3566 swed.interpol.tjd_nut0 = 0;
3567 swed.interpol.tjd_nut2 = 0;
3568 swed.interpol.nut_dpsi0 = 0;
3569 swed.interpol.nut_dpsi1 = 0;
3570 swed.interpol.nut_dpsi2 = 0;
3571 swed.interpol.nut_deps0 = 0;
3572 swed.interpol.nut_deps1 = 0;
3573 swed.interpol.nut_deps2 = 0;
3574 }
3575
3576 /* sidereal time, without eps and nut as parameters.
3577 * tjd must be UT !!!
3578 * for more informsation, see comment with swe_sidtime0()
3579 */
swe_sidtime(double tjd_ut)3580 double CALL_CONV swe_sidtime(double tjd_ut)
3581 {
3582 int i;
3583 double eps, nutlo[2], tsid;
3584 double tjde;
3585 /* delta t adjusted to default tidal acceleration of the moon */
3586 tjde = tjd_ut + swe_deltat_ex(tjd_ut, -1, NULL);
3587 swi_init_swed_if_start();
3588 eps = swi_epsiln(tjde, 0) * RADTODEG;
3589 swi_nutation(tjde, 0, nutlo);
3590 for (i = 0; i < 2; i++)
3591 nutlo[i] *= RADTODEG;
3592 tsid = swe_sidtime0(tjd_ut, eps + nutlo[1], nutlo[0]);
3593 return tsid;
3594 }
3595
3596 /* SWISSEPH
3597 * generates name of ephemeris file
3598 * file name looks as follows:
3599 * swephpl.m30, where
3600 *
3601 * "sweph" "swiss ephemeris"
3602 * "pl","mo","as" planet, moon, or asteroid
3603 * "m" or "_" BC or AD
3604 *
3605 * "30" start century
3606 * tjd = ephemeris file for which julian day
3607 * ipli = number of planet
3608 * fname = ephemeris file name
3609 */
swi_gen_filename(double tjd,int ipli,char * fname)3610 void swi_gen_filename(double tjd, int ipli, char *fname)
3611 {
3612 int icty;
3613 int ncties = (int) NCTIES;
3614 short gregflag;
3615 int jmon, jday, jyear, sgn;
3616 double jut;
3617 char *sform;
3618 switch(ipli) {
3619 case SEI_MOON:
3620 strcpy(fname, "semo");
3621 break;
3622 case SEI_EMB:
3623 case SEI_MERCURY:
3624 case SEI_VENUS:
3625 case SEI_MARS:
3626 case SEI_JUPITER:
3627 case SEI_SATURN:
3628 case SEI_URANUS:
3629 case SEI_NEPTUNE:
3630 case SEI_PLUTO:
3631 case SEI_SUNBARY:
3632 strcpy(fname, "sepl");
3633 break;
3634 case SEI_CERES:
3635 case SEI_PALLAS:
3636 case SEI_JUNO:
3637 case SEI_VESTA:
3638 case SEI_CHIRON:
3639 case SEI_PHOLUS:
3640 strcpy(fname, "seas");
3641 break;
3642 default: /* asteroid or planetary moon */
3643 if (ipli > SE_PLMOON_OFFSET && ipli < SE_AST_OFFSET) {
3644 sprintf(fname, "sat%ssepm%d.%s", DIR_GLUE, ipli, SE_FILE_SUFFIX);
3645 } else {
3646 sform = "ast%d%sse%05d.%s";
3647 if (ipli - SE_AST_OFFSET > 99999)
3648 sform = "ast%d%ss%06d.%s";
3649 sprintf(fname, sform, (ipli - SE_AST_OFFSET) / 1000, DIR_GLUE, ipli - SE_AST_OFFSET, SE_FILE_SUFFIX);
3650 }
3651 return; /* asteroids or planetary moons: only one file 3000 bc - 3000 ad */
3652 /* break; */
3653 }
3654 /* century of tjd */
3655 /* if tjd > 1600 then gregorian calendar */
3656 if (tjd >= 2305447.5) {
3657 gregflag = TRUE;
3658 swe_revjul(tjd, gregflag, &jyear, &jmon, &jday, &jut);
3659 /* else julian calendar */
3660 } else {
3661 gregflag = FALSE;
3662 swe_revjul(tjd, gregflag, &jyear, &jmon, &jday, &jut);
3663 }
3664 /* start century of file containing tjd */
3665 if (jyear < 0)
3666 sgn = -1;
3667 else
3668 sgn = 1;
3669 icty = jyear / 100;
3670 if (sgn < 0 && jyear % 100 != 0)
3671 icty -=1;
3672 while(icty % ncties != 0)
3673 icty--;
3674 #if 0
3675 if (icty < BEG_YEAR / 100)
3676 icty = BEG_YEAR / 100;
3677 if (icty >= END_YEAR / 100)
3678 icty = END_YEAR / 100 - ncties;
3679 #endif
3680 /* B.C. or A.D. */
3681 if (icty < 0)
3682 strcat(fname, "m");
3683 else
3684 strcat(fname, "_");
3685 icty = abs(icty);
3686 sprintf(fname + strlen(fname), "%02d.%s", icty, SE_FILE_SUFFIX);
3687 #if 0
3688 printf("fname %s\n", fname);
3689 fflush(stdout);
3690 #endif
3691 }
3692
3693 /**************************************************************
3694 cut the string s at any char in cutlist; put pointers to partial strings
3695 into cpos[0..n-1], return number of partial strings;
3696 if less than nmax fields are found, the first empty pointer is
3697 set to NULL.
3698 More than one character of cutlist in direct sequence count as one
3699 separator only! cut_str_any("word,,,word2",","..) cuts only two parts,
3700 cpos[0] = "word" and cpos[1] = "word2".
3701 If more than nmax fields are found, nmax is returned and the
3702 last field nmax-1 rmains un-cut.
3703 **************************************************************/
swi_cutstr(char * s,char * cutlist,char * cpos[],int nmax)3704 int swi_cutstr(char *s, char *cutlist, char *cpos[], int nmax)
3705 {
3706 int n = 1;
3707 cpos [0] = s;
3708 while (*s != '\0') {
3709 if ((strchr(cutlist, (int) *s) != NULL) && n < nmax) {
3710 *s = '\0';
3711 while (*(s + 1) != '\0' && strchr (cutlist, (int) *(s + 1)) != NULL) s++;
3712 cpos[n++] = s + 1;
3713 }
3714 if (*s == '\n' || *s == '\r') { /* treat nl or cr like end of string */
3715 *s = '\0';
3716 break;
3717 }
3718 s++;
3719 }
3720 if (n < nmax) cpos[n] = NULL;
3721 return (n);
3722 } /* cutstr */
3723
swi_right_trim(char * s)3724 char *swi_right_trim(char *s)
3725 {
3726 char *sp = s + strlen(s) - 1;
3727 // while (isspace((int)(unsigned char) *sp) && sp >= s)
3728 while (sp >= s && isspace((int)(unsigned char) *sp))
3729 *sp-- = '\0';
3730 return s;
3731 }
3732
swi_strnlen(const char * str,size_t n)3733 size_t swi_strnlen(const char *str, size_t n) {
3734 const char * stop = (char *)memchr(str, '\0', n);
3735 return stop ? stop - str : n;
3736 }
3737
3738 /*
3739 * The following C code (by Rob Warnock rpw3@sgi.com) does CRC-32 in
3740 * BigEndian/BigEndian byte/bit order. That is, the data is sent most
3741 * significant byte first, and each of the bits within a byte is sent most
3742 * significant bit first, as in FDDI. You will need to twiddle with it to do
3743 * Ethernet CRC, i.e., BigEndian/LittleEndian byte/bit order.
3744 *
3745 * The CRCs this code generates agree with the vendor-supplied Verilog models
3746 * of several of the popular FDDI "MAC" chips.
3747 */
3748 static TLS uint32 crc32_table[256];
3749 /* Initialized first time "crc32()" is called. If you prefer, you can
3750 * statically initialize it at compile time. [Another exercise.]
3751 */
3752
swi_crc32(unsigned char * buf,int len)3753 uint32 swi_crc32(unsigned char *buf, int len)
3754 {
3755 unsigned char *p;
3756 uint32 crc;
3757 if (!crc32_table[1]) /* if not already done, */
3758 init_crc32(); /* build table */
3759 crc = 0xffffffff; /* preload shift register, per CRC-32 spec */
3760 for (p = buf; len > 0; ++p, --len)
3761 crc = (crc << 8) ^ crc32_table[(crc >> 24) ^ *p];
3762 return ~crc; /* transmit complement, per CRC-32 spec */
3763 }
3764
3765 /*
3766 * Build auxiliary table for parallel byte-at-a-time CRC-32.
3767 */
3768 #define CRC32_POLY 0x04c11db7 /* AUTODIN II, Ethernet, & FDDI */
3769
init_crc32(void)3770 static void init_crc32(void)
3771 {
3772 int32 i, j;
3773 uint32 c;
3774 for (i = 0; i < 256; ++i) {
3775 for (c = i << 24, j = 8; j > 0; --j)
3776 c = c & 0x80000000 ? (c << 1) ^ CRC32_POLY : (c << 1);
3777 crc32_table[i] = c;
3778 }
3779 }
3780
3781 /*******************************************************
3782 * other functions from swephlib.c;
3783 * they are not needed for Swiss Ephemeris,
3784 * but may be useful to former Placalc users.
3785 ********************************************************/
3786
3787 /************************************
3788 normalize argument into interval [0..DEG360]
3789 *************************************/
swe_csnorm(centisec p)3790 centisec CALL_CONV swe_csnorm(centisec p)
3791 {
3792 if (p < 0)
3793 do { p += DEG360; } while (p < 0);
3794 else if (p >= DEG360)
3795 do { p -= DEG360; } while (p >= DEG360);
3796 return (p);
3797 }
3798
3799 /************************************
3800 distance in centisecs p1 - p2
3801 normalized to [0..360[
3802 **************************************/
swe_difcsn(centisec p1,centisec p2)3803 centisec CALL_CONV swe_difcsn (centisec p1, centisec p2)
3804 {
3805 return (swe_csnorm(p1 - p2));
3806 }
3807
swe_difdegn(double p1,double p2)3808 double CALL_CONV swe_difdegn (double p1, double p2)
3809 {
3810 return (swe_degnorm(p1 - p2));
3811 }
3812
3813 /************************************
3814 distance in centisecs p1 - p2
3815 normalized to [-180..180[
3816 **************************************/
swe_difcs2n(centisec p1,centisec p2)3817 centisec CALL_CONV swe_difcs2n(centisec p1, centisec p2)
3818 { centisec dif;
3819 dif = swe_csnorm(p1 - p2);
3820 if (dif >= DEG180) return (dif - DEG360);
3821 return (dif);
3822 }
3823
swe_difdeg2n(double p1,double p2)3824 double CALL_CONV swe_difdeg2n(double p1, double p2)
3825 { double dif;
3826 dif = swe_degnorm(p1 - p2);
3827 if (dif >= 180.0) return (dif - 360.0);
3828 return (dif);
3829 }
3830
swe_difrad2n(double p1,double p2)3831 double CALL_CONV swe_difrad2n(double p1, double p2)
3832 { double dif;
3833 dif = swe_radnorm(p1 - p2);
3834 if (dif >= TWOPI / 2) return (dif - TWOPI);
3835 return (dif);
3836 }
3837
3838 /*************************************
3839 round second, but at 29.5959 always down
3840 *************************************/
swe_csroundsec(centisec x)3841 centisec CALL_CONV swe_csroundsec(centisec x)
3842 {
3843 centisec t;
3844 t = (x + 50) / 100 *100L; /* round to seconds */
3845 if (t > x && t % DEG30 == 0) /* was rounded up to next sign */
3846 t = x / 100 * 100L; /* round last second of sign downwards */
3847 return (t);
3848 }
3849
3850 /*************************************
3851 double to int32 with rounding, no overflow check
3852 *************************************/
swe_d2l(double x)3853 int32 CALL_CONV swe_d2l(double x)
3854 {
3855 if (x >=0)
3856 return ((int32) (x + 0.5));
3857 else
3858 return (- (int32) (0.5 - x));
3859 }
3860
3861 /*
3862 * monday = 0, ... sunday = 6
3863 */
swe_day_of_week(double jd)3864 int CALL_CONV swe_day_of_week(double jd)
3865 {
3866 return (((int) floor (jd - 2433282 - 1.5) %7) + 7) % 7;
3867 }
3868
swe_cs2timestr(CSEC t,int sep,AS_BOOL suppressZero,char * a)3869 char *CALL_CONV swe_cs2timestr(CSEC t, int sep, AS_BOOL suppressZero, char *a)
3870 /* does not suppress zeros in hours or minutes */
3871 {
3872 /* static char a[9];*/
3873 centisec h,m,s;
3874 strcpy (a, " ");
3875 a[2] = a [5] = sep;
3876 t = ((t + 50) / 100) % (24L *3600L); /* round to seconds */
3877 s = t % 60L;
3878 m = (t / 60) % 60L;
3879 h = t / 3600 % 100L;
3880 if (s == 0 && suppressZero)
3881 a[5] = '\0';
3882 else {
3883 a [6] = (char) (s / 10 + '0');
3884 a [7] = (char) (s % 10 + '0');
3885 };
3886 a [0] = (char) (h / 10 + '0');
3887 a [1] = (char) (h % 10 + '0');
3888 a [3] = (char) (m / 10 + '0');
3889 a [4] = (char) (m % 10 + '0');
3890 return (a);
3891 } /* swe_cs2timestr() */
3892
swe_cs2lonlatstr(CSEC t,char pchar,char mchar,char * sp)3893 char *CALL_CONV swe_cs2lonlatstr(CSEC t, char pchar, char mchar, char *sp)
3894 {
3895 char a[10]; /* must be initialized at each call */
3896 char *aa;
3897 centisec h,m,s;
3898 strcpy (a, " ' ");
3899 /* mask dddEmm'ss" */
3900 if (t < 0 ) pchar = mchar;
3901 t = (ABS4 (t) + 50) / 100; /* round to seconds */
3902 s = t % 60L;
3903 m = t / 60 % 60L;
3904 h = t / 3600 % 1000L;
3905 if (s == 0)
3906 a[6] = '\0'; /* cut off seconds */
3907 else {
3908 a [7] = (char) (s / 10 + '0');
3909 a [8] = (char) (s % 10 + '0');
3910 }
3911 a [3] = pchar;
3912 if (h > 99) a [0] = (char) (h / 100 + '0');
3913 if (h > 9) a [1] = (char) (h % 100 / 10 + '0');
3914 a [2] = (char) (h % 10 + '0');
3915 a [4] = (char) (m / 10 + '0');
3916 a [5] = (char) (m % 10 + '0');
3917 aa = a;
3918 while (*aa == ' ') aa++;
3919 strcpy(sp, aa);
3920 return (sp);
3921 } /* swe_cs2lonlatstr() */
3922
swe_cs2degstr(CSEC t,char * a)3923 char *CALL_CONV swe_cs2degstr(CSEC t, char *a)
3924 /* does suppress leading zeros in degrees */
3925 {
3926 /* char a[9]; must be initialized at each call */
3927 centisec h,m,s;
3928 t = t / 100 % (30L*3600L); /* truncate to seconds */
3929 s = t % 60L;
3930 m = t / 60 % 60L;
3931 h = t / 3600 % 100L; /* only 0..99 degrees */
3932 sprintf(a, "%2d%s%02d'%02d", h, ODEGREE_STRING, m, s);
3933 return (a);
3934 } /* swe_cs2degstr() */
3935
3936 /******************************************************************
3937 * decimal degrees in zodiac to nakshatra position, deg, min, sec *
3938 * for definition of input see function swe_split_deg().
3939 * output:
3940 * ideg degrees,
3941 * imin minutes,
3942 * isec seconds,
3943 * dsecfr fraction of seconds
3944 * inak nakshatra number;
3945 ******************************************************************/
split_deg_nakshatra(double ddeg,int32 roundflag,int32 * ideg,int32 * imin,int32 * isec,double * dsecfr,int32 * inak)3946 static void split_deg_nakshatra(double ddeg, int32 roundflag, int32 *ideg, int32 *imin, int32 *isec, double *dsecfr, int32 *inak)
3947 {
3948 double dadd = 0;
3949 double dnakshsize = 13.33333333333333;
3950 double ddeghelp = fmod(ddeg, dnakshsize);
3951 *inak = 1;
3952 if (ddeg < 0) {
3953 *inak = -1;
3954 ddeg = 0;
3955 }
3956 // Sheoran "Vedic" ayanamsha: 0 Aries = 3°20 Ashvini
3957 if ((swed.sidd.sid_mode & SE_SIDM_TRUE_SHEORAN) == SE_SIDM_TRUE_SHEORAN)
3958 ddeg = swe_degnorm(ddeg + 3.33333333333333);
3959 if (roundflag & SE_SPLIT_DEG_ROUND_DEG) {
3960 dadd = 0.5;
3961 } else if (roundflag & SE_SPLIT_DEG_ROUND_MIN) {
3962 dadd = 0.5 / 60;
3963 } else if (roundflag & SE_SPLIT_DEG_ROUND_SEC) {
3964 dadd = 0.5 / 3600;
3965 }
3966 if (roundflag & SE_SPLIT_DEG_KEEP_DEG) {
3967 if ((int32) (ddeghelp + dadd) - (int32) ddeghelp > 0)
3968 dadd = 0;
3969 } else if (roundflag & SE_SPLIT_DEG_KEEP_SIGN) {
3970 if (ddeghelp + dadd >= dnakshsize)
3971 dadd = 0;
3972 }
3973 ddeg += dadd;
3974 *inak = (int32) (ddeg / dnakshsize);
3975 if (*inak == 27) *inak = 0; // with rounding up from 359.9999
3976 ddeg = fmod(ddeg, dnakshsize);
3977 *ideg = (int32) ddeg;
3978 ddeg -= *ideg;
3979 *imin = (int32) (ddeg * 60);
3980 ddeg -= *imin / 60.0;
3981 *isec = (int32) (ddeg * 3600);
3982 if (!(roundflag & (SE_SPLIT_DEG_ROUND_DEG | SE_SPLIT_DEG_ROUND_MIN | SE_SPLIT_DEG_ROUND_SEC))) {
3983 *dsecfr = ddeg * 3600 - *isec;
3984 }
3985 } /* end split_deg_nakshtra */
3986
3987 /************************************************************
3988 * splitting decimal degrees into (zod.sign,) deg, min, sec. *
3989 * input:
3990 * ddeg decimal degrees, ecliptic longitude
3991 * roundflag by default there is no rounding. if rounding is
3992 * required, the following bits can be set:
3993 # define SE_SPLIT_DEG_ROUND_SEC 1
3994 # define SE_SPLIT_DEG_ROUND_MIN 2
3995 # define SE_SPLIT_DEG_ROUND_DEG 4
3996 # define SE_SPLIT_DEG_ZODIACAL 8 * split into zodiac signs
3997 # define SE_SPLIT_DEG_NAKSHATRA 1024 * split into nakshatras *
3998 # define SE_SPLIT_DEG_KEEP_SIGN 16 * don't round to next zodiac sign,
3999 * e.g. 29.9999998 will be rounded
4000 * to 29°59'59" (or 29°59' or 29°) *
4001 * or next nakshatra:
4002 * e.g. 13.3333332 will be rounded
4003 * to 13°19'59" (or 13°19' or 13°) *
4004 # define SE_SPLIT_DEG_KEEP_DEG 32 * don't round to next degree
4005 * e.g. 10.9999999 will be rounded
4006 * to 10d59'59" (or 10d59' or 10d) *
4007 * output:
4008 * ideg degrees,
4009 * imin minutes,
4010 * isec seconds,
4011 * dsecfr fraction of seconds
4012 * isgn zodiac sign number;
4013 * or +/- sign
4014 *
4015 *********************************************************/
swe_split_deg(double ddeg,int32 roundflag,int32 * ideg,int32 * imin,int32 * isec,double * dsecfr,int32 * isgn)4016 void CALL_CONV swe_split_deg(double ddeg, int32 roundflag, int32 *ideg, int32 *imin, int32 *isec, double *dsecfr, int32 *isgn)
4017 {
4018 double dadd = 0;
4019 *isgn = 1;
4020 if (ddeg < 0) {
4021 *isgn = -1;
4022 ddeg = -ddeg;
4023 } else if (roundflag & SE_SPLIT_DEG_NAKSHATRA) {
4024 split_deg_nakshatra(ddeg, roundflag, ideg, imin, isec, dsecfr, isgn);
4025 return;
4026 }
4027 if (roundflag & SE_SPLIT_DEG_ROUND_DEG) {
4028 dadd = 0.5;
4029 } else if (roundflag & SE_SPLIT_DEG_ROUND_MIN) {
4030 dadd = 0.5 / 60.0;
4031 } else if (roundflag & SE_SPLIT_DEG_ROUND_SEC) {
4032 dadd = 0.5 / 3600.0;
4033 }
4034 if (roundflag & SE_SPLIT_DEG_KEEP_DEG) {
4035 if ((int32) (ddeg + dadd) - (int32) ddeg > 0)
4036 dadd = 0;
4037 } else if (roundflag & SE_SPLIT_DEG_KEEP_SIGN) {
4038 if (fmod(ddeg, 30) + dadd >= 30)
4039 dadd = 0;
4040 }
4041 ddeg += dadd;
4042 if (roundflag & SE_SPLIT_DEG_ZODIACAL) {
4043 *isgn = (int32) (ddeg / 30);
4044 if (*isgn == 12) // 360° = 0°
4045 *isgn = 0;
4046 ddeg = fmod(ddeg, 30);
4047 }
4048 *ideg = (int32) ddeg;
4049 ddeg -= *ideg;
4050 *imin = (int32) (ddeg * 60);
4051 ddeg -= *imin / 60.0;
4052 *isec = (int32) (ddeg * 3600);
4053 if (!(roundflag & (SE_SPLIT_DEG_ROUND_DEG | SE_SPLIT_DEG_ROUND_MIN | SE_SPLIT_DEG_ROUND_SEC))) {
4054 *dsecfr = ddeg * 3600 - *isec;
4055 }
4056 } /* end split_deg */
4057
swi_kepler(double E,double M,double ecce)4058 double swi_kepler(double E, double M, double ecce)
4059 {
4060 double dE = 1, E0;
4061 double x;
4062 /* simple formula for small eccentricities */
4063 if (ecce < 0.4) {
4064 while(dE > 1e-12) {
4065 E0 = E;
4066 E = M + ecce * sin(E0);
4067 dE = fabs(E - E0);
4068 }
4069 /* complicated formula for high eccentricities */
4070 } else {
4071 while(dE > 1e-12) {
4072 E0 = E;
4073 /*
4074 * Alois 21-jul-2000: workaround an optimizer problem in gcc
4075 * swi_mod2PI sees very small negative argument e-322 and returns +2PI;
4076 * we avoid swi_mod2PI for small x.
4077 */
4078 x = (M + ecce * sin(E0) - E0) / (1 - ecce * cos(E0));
4079 dE = fabs(x);
4080 if (dE < 1e-2) {
4081 E = E0 + x;
4082 } else {
4083 E = swi_mod2PI(E0 + x);
4084 dE = fabs(E - E0);
4085 }
4086 }
4087 }
4088 return E;
4089 }
4090
swi_FK4_FK5(double * xp,double tjd)4091 void swi_FK4_FK5(double *xp, double tjd)
4092 {
4093 AS_BOOL correct_speed = TRUE;
4094 if (xp[0] == 0 && xp[1] == 0 && xp[2] == 0)
4095 return;
4096 /* with zero speed, we assume that it should be really zero */
4097 if (xp[3] == 0)
4098 correct_speed = FALSE;
4099 swi_cartpol_sp(xp, xp);
4100 /* according to Expl.Suppl., p. 167f. */
4101 xp[0] += (0.035 + 0.085 * (tjd - B1950) / 36524.2198782) / 3600 * 15 * DEGTORAD;
4102 if (correct_speed)
4103 xp[3] += (0.085 / 36524.2198782) / 3600 * 15 * DEGTORAD;
4104 swi_polcart_sp(xp, xp);
4105 }
4106
swi_FK5_FK4(double * xp,double tjd)4107 void swi_FK5_FK4(double *xp, double tjd)
4108 {
4109 if (xp[0] == 0 && xp[1] == 0 && xp[2] == 0)
4110 return;
4111 swi_cartpol_sp(xp, xp);
4112 /* according to Expl.Suppl., p. 167f. */
4113 xp[0] -= (0.035 + 0.085 * (tjd - B1950) / 36524.2198782) / 3600 * 15 * DEGTORAD;
4114 xp[3] -= (0.085 / 36524.2198782) / 3600 * 15 * DEGTORAD;
4115 swi_polcart_sp(xp, xp);
4116 }
4117
4118 /* function for inhouse testing only */
set_astro_models(char * samod)4119 void set_astro_models(char *samod)
4120 {
4121 int *pmodel = &(swed.astro_models[0]);
4122 char *sp, *sp2;
4123 int i = 0;
4124 swi_init_swed_if_start();
4125 sp = samod;
4126 pmodel[0] = atoi(sp);
4127 i++;
4128 while((sp2 = strchr(sp, ',')) != NULL && i < NSE_MODELS) {
4129 sp = sp2 + 1;
4130 pmodel[i] = atoi(sp);
4131 i++;
4132 }
4133 }
4134
4135
4136 /*
4137 * Function for inhouse testing of old SE versions.
4138 *
4139 * Values of the following macros are defined in swephexp.h
4140 * and must be identical to the ones used in the defines below
4141 D1 SEMOD_DELTAT_STEPHENSON_MORRISON_1984
4142 D2 SEMOD_DELTAT_STEPHENSON_1997
4143 D3 SEMOD_DELTAT_STEPHENSON_MORRISON_2004
4144 D4 SEMOD_DELTAT_ESPENAK_MEEUS_2006
4145 D5 SEMOD_DELTAT_STEPHENSON_ETC_2016
4146
4147 P1 SEMOD_PREC_IAU_1976
4148 P2 SEMOD_PREC_LASKAR_1986
4149 P3 SEMOD_PREC_WILL_EPS_LASK
4150 P4 SEMOD_PREC_WILLIAMS_1994
4151 P5 SEMOD_PREC_SIMON_1994
4152 P6 SEMOD_PREC_IAU_2000
4153 P7 SEMOD_PREC_BRETAGNON_2003
4154 P8 SEMOD_PREC_IAU_2006
4155 P9 SEMOD_PREC_VONDRAK_2011
4156 P10 SEMOD_PREC_OWEN_1990
4157 P11 SEMOD_PREC_NEWCOMB
4158
4159 N1 SEMOD_NUT_IAU_1980
4160 N2 SEMOD_NUT_IAU_CORR_1987
4161 N3 SEMOD_NUT_IAU_2000A
4162 N4 SEMOD_NUT_IAU_2000B
4163 N5 SEMOD_NUT_WOOLARD
4164
4165 B1 SEMOD_BIAS_NONE
4166 B2 SEMOD_BIAS_IAU2000
4167 B3 SEMOD_BIAS_IAU2006
4168
4169 S1 SEMOD_SIDT_IAU_1976
4170 S2 SEMOD_SIDT_IAU_2006
4171 S3 SEMOD_SIDT_IERS_CONV_2010
4172 S4 SEMOD_SIDT_LONGTERM
4173 * D P P N B J J S
4174 */
4175 # define AMODELS_SE_1_00 "1,3,1,1,1,0,0,1"
4176 # define AMODELS_SE_1_64 "2,3,1,1,1,0,0,1"
4177 # define AMODELS_SE_1_70 "2,8,8,4,2,0,0,2"
4178 # define AMODELS_SE_1_72 "3,8,8,4,2,0,0,2"
4179 # define AMODELS_SE_1_77 "4,8,8,4,2,0,0,2"
4180 # define AMODELS_SE_1_78 "4,9,9,4,2,0,0,2"
4181 # define AMODELS_SE_1_80 "4,9,9,4,3,0,0,1" /* note sid. time (S)! */
4182 # define AMODELS_SE_2_00 "4,9,9,4,3,0,0,4"
4183 # define AMODELS_SE_2_06 "5,9,9,4,3,0,0,4"
swe_set_astro_models(char * samod,int32 iflag)4184 void CALL_CONV swe_set_astro_models(char *samod, int32 iflag)
4185 {
4186 double dversion;
4187 char s[30], *sp;
4188 swi_init_swed_if_start();
4189 if (*samod != '\0' && isdigit(*samod)) {
4190 set_astro_models(samod);
4191 } else if (*samod == '\0' || strncmp(samod, "SE", 2) == 0) {
4192 strncpy(s, samod, 20);
4193 s[20] = '\0';
4194 if ((sp = strchr(s + 5, '.')) != NULL) // remove second '.' in "SE2.05.01"
4195 swi_strcpy(sp, sp+1);
4196 if ((sp = strchr(s + 5, 'b')) != NULL) // remove 'b' in "SE2.05.02b04"
4197 swi_strcpy(sp, sp+1);
4198 dversion = atof(s + 2);
4199 if (dversion == 0)
4200 dversion = atof(SE_VERSION);
4201 if (dversion >= 2.06) {
4202 set_astro_models(AMODELS_SE_2_06);
4203 } else if (dversion >= 2.01) {
4204 set_astro_models(AMODELS_SE_2_00);
4205 } else if (dversion >= 2.00) {
4206 set_astro_models(AMODELS_SE_2_00);
4207 if (swi_get_denum(SEI_SUN, iflag) == 431)
4208 swe_set_tid_acc(SE_TIDAL_DE406);
4209 } else if (dversion >= 1.80) {
4210 set_astro_models(AMODELS_SE_1_80);
4211 swe_set_tid_acc(SE_TIDAL_DE406);
4212 } else if (dversion >= 1.78) {
4213 set_astro_models(AMODELS_SE_1_78);
4214 swe_set_tid_acc(SE_TIDAL_DE406);
4215 } else if (dversion >= 1.77) {
4216 set_astro_models(AMODELS_SE_1_77);
4217 swe_set_tid_acc(SE_TIDAL_DE406);
4218 } else if (dversion >= 1.72) {
4219 set_astro_models(AMODELS_SE_1_72);
4220 swe_set_tid_acc(-25.7376);
4221 } else if (dversion >= 1.70) {
4222 set_astro_models(AMODELS_SE_1_70);
4223 swe_set_tid_acc(-25.7376);
4224 } else if (dversion >= 1.64) {
4225 set_astro_models(AMODELS_SE_1_64);
4226 swe_set_tid_acc(-25.7376);
4227 } else {
4228 set_astro_models(AMODELS_SE_1_00);
4229 swe_set_tid_acc(-25.7376);
4230 }
4231 }
4232 }
4233
4234 /* function for inhouse testing only */
get_precession_model(int precmod,int32 iflag,char * s)4235 static void get_precession_model(int precmod, int32 iflag, char *s)
4236 {
4237 if (precmod == 0)
4238 precmod = SEMOD_PREC_DEFAULT;
4239 if (iflag & SEFLG_JPLEPH) {
4240 if (iflag & SEFLG_JPLHOR) {
4241 strcpy(s, "IAU 1976 (Lieske) / Owen 1990 before 1799");
4242 return;
4243 }
4244 if (iflag & SEFLG_JPLHOR_APPROX) {
4245 strcpy(s, "Vondrak 2011 / IAU 1976 (Lieske) before 1962 / Owen 1990 before 1799");
4246 return;
4247 }
4248 }
4249 switch(precmod) {
4250 case SEMOD_PREC_IAU_1976:
4251 strcpy(s, "IAU 1976 (Lieske)");
4252 break;
4253 case SEMOD_PREC_IAU_2000:
4254 strcpy(s, "IAU 2000 (Lieske 1976, Mathews 2002)");
4255 break;
4256 case SEMOD_PREC_IAU_2006:
4257 strcpy(s, "IAU 2006 (Capitaine & alii)");
4258 break;
4259 case SEMOD_PREC_BRETAGNON_2003:
4260 strcpy(s, "Bretagnon 2003");
4261 break;
4262 case SEMOD_PREC_LASKAR_1986:
4263 strcpy(s, "Laskar 1986");
4264 break;
4265 case SEMOD_PREC_SIMON_1994:
4266 strcpy(s, "Simon 1994");
4267 break;
4268 case SEMOD_PREC_WILLIAMS_1994:
4269 strcpy(s, "Williams 1994");
4270 break;
4271 case SEMOD_PREC_WILL_EPS_LASK:
4272 strcpy(s, "Williams 1994 / Epsilon Laskar 1986");
4273 break;
4274 case SEMOD_PREC_OWEN_1990:
4275 strcpy(s, "Owen 1990");
4276 break;
4277 case SEMOD_PREC_NEWCOMB:
4278 strcpy(s, "Newcomb 1895");
4279 break;
4280 case SEMOD_PREC_VONDRAK_2011:
4281 strcpy(s, "Vondrák 2011");
4282 break;
4283 default:
4284 break;
4285 }
4286 }
4287
4288 /* function for inhouse testing only */
get_deltat_model(int dtmod,char * s)4289 static void get_deltat_model(int dtmod, char *s)
4290 {
4291 if (dtmod == 0)
4292 dtmod = SEMOD_DELTAT_DEFAULT;
4293 switch(dtmod) {
4294 case SEMOD_DELTAT_ESPENAK_MEEUS_2006:
4295 strcpy(s, "Espenak/Meeus 2006 (before 1633)");
4296 break;
4297 case SEMOD_DELTAT_STEPHENSON_MORRISON_2004:
4298 strcpy(s, "Stephenson/Morrison 2004 (before 1600)");
4299 break;
4300 case SEMOD_DELTAT_STEPHENSON_1997:
4301 strcpy(s, "Stephenson 1997 (before 1600)");
4302 break;
4303 case SEMOD_DELTAT_STEPHENSON_MORRISON_1984:
4304 strcpy(s, "Stephenson/Morrison 1984 (before 1600)");
4305 break;
4306 case SEMOD_DELTAT_STEPHENSON_ETC_2016:
4307 strcpy(s, "Stephenson/Morrison/Hohenkerk 2016 (before 1955)");
4308 break;
4309 }
4310 }
4311
4312 /* function for inhouse testing only */
get_nutation_model(int nutmod,int32 iflag,char * s)4313 static void get_nutation_model(int nutmod, int32 iflag, char *s)
4314 {
4315 int jplhormod = swed.astro_models[SE_MODEL_JPLHOR_MODE];
4316 int jplhoramod = swed.astro_models[SE_MODEL_JPLHORA_MODE];
4317 if (jplhormod == 0)
4318 jplhormod = SEMOD_JPLHOR_DEFAULT;
4319 if (jplhoramod == 0)
4320 jplhoramod = SEMOD_JPLHORA_DEFAULT;
4321 if (nutmod == 0)
4322 nutmod = SEMOD_NUT_DEFAULT;
4323 switch(nutmod) {
4324 case SEMOD_NUT_WOOLARD:
4325 strcpy(s, "Woolard 1953");
4326 break;
4327 case SEMOD_NUT_IAU_1980:
4328 strcpy(s, "IAU 1980 (Wahr)");
4329 break;
4330 case SEMOD_NUT_IAU_CORR_1987:
4331 strcpy(s, "Herring 1986");
4332 break;
4333 case SEMOD_NUT_IAU_2000A:
4334 strcpy(s, "IAU 2000A (Mathews)");
4335 break;
4336 case SEMOD_NUT_IAU_2000B:
4337 strcpy(s, "IAU 2000B (Mathews)");
4338 break;
4339 }
4340 if (iflag & SEFLG_JPLEPH) {
4341 if (iflag & SEFLG_JPLHOR)
4342 strcpy(s, "IAU 1980 (Wahr)");
4343 if (iflag & SEFLG_JPLHOR) {
4344 strcat(s, "\n+ daily corrections to dpsi/deps 1962-today");
4345 if (jplhormod == SEMOD_JPLHOR_LONG_AGREEMENT)
4346 strcat(s, "\n good agreement with JPL Horizons between 1800 and today");
4347 else
4348 strcat(s, "\n defaults to SEFLG_JPLEPH_APPROX before 1962");
4349 } else if (iflag & SEFLG_JPLHOR_APPROX){
4350 strcat(s, "\n+ some corrections, approximating JPL Horizons");
4351 if (jplhoramod == SEMOD_JPLHORA_1)
4352 strcat(s, " (SEMOD_JPLHORA_1)");
4353 else if (jplhoramod == SEMOD_JPLHORA_2)
4354 strcat(s, " (SEMOD_JPLHORA_2)");
4355 else
4356 strcat(s, " (SEMOD_JPLHORA_3)");
4357 }
4358 }
4359 }
4360
4361 /* function for inhouse testing only */
get_frame_bias_model(int biasmod,char * s)4362 static void get_frame_bias_model(int biasmod, char *s)
4363 {
4364 if (biasmod == 0)
4365 biasmod = SEMOD_BIAS_DEFAULT;
4366 switch(biasmod) {
4367 case SEMOD_BIAS_IAU2000:
4368 strcpy(s, "IAU 2000");
4369 break;
4370 case SEMOD_BIAS_IAU2006:
4371 strcpy(s, "IAU 2006");
4372 break;
4373 case SEMOD_BIAS_NONE:
4374 strcpy(s, "none");
4375 break;
4376 }
4377 }
4378
4379 /* function for inhouse testing only */
get_sidt_model(int sidtmod,char * s)4380 static void get_sidt_model(int sidtmod, char *s)
4381 {
4382 if (sidtmod == 0)
4383 sidtmod = SEMOD_SIDT_DEFAULT;
4384 switch(sidtmod) {
4385 case SEMOD_SIDT_IAU_1976:
4386 strcpy(s, "IAU 1976");
4387 break;
4388 case SEMOD_SIDT_IAU_2006:
4389 strcpy(s, "IAU 2006 (Capitaine 2003)");
4390 break;
4391 case SEMOD_SIDT_IERS_CONV_2010:
4392 strcpy(s, "IERS Convention 2010");
4393 break;
4394 case SEMOD_SIDT_LONGTERM:
4395 strcpy(s, "IERS Convention 2010 + long-term extension by Astrodienst");
4396 break;
4397 }
4398 }
4399
4400 /* function for inhouse testing only */
swe_get_astro_models(char * samod,char * sdet,int32 iflag)4401 void CALL_CONV swe_get_astro_models(char *samod, char *sdet, int32 iflag)
4402 {
4403 int i, imod;
4404 int *pmodel = &(swed.astro_models[0]);
4405 char s[AS_MAXCH], samod0[AS_MAXCH];
4406 AS_BOOL list_all_models = FALSE;
4407 if (samod != NULL) {
4408 if (strchr(samod, '+') != NULL)
4409 list_all_models = TRUE;
4410 swe_set_astro_models(samod, iflag);
4411 }
4412 *samod0 = '\0';
4413 for (i = 0; i < NSE_MODELS; i++) {
4414 imod = pmodel[i];
4415 switch(i) {
4416 case SE_MODEL_PREC_LONGTERM:
4417 if (imod == SEMOD_PREC_DEFAULT) imod = 0; break;
4418 case SE_MODEL_PREC_SHORTTERM:
4419 if (imod == SEMOD_PREC_DEFAULT_SHORT) imod = 0; break;
4420 case SE_MODEL_NUT:
4421 if (imod == SEMOD_NUT_DEFAULT) imod = 0; break;
4422 case SE_MODEL_SIDT:
4423 if (imod == SEMOD_SIDT_DEFAULT) imod = 0; break;
4424 case SE_MODEL_BIAS:
4425 if (imod == SEMOD_BIAS_DEFAULT) imod = 0; break;
4426 case SE_MODEL_JPLHOR_MODE:
4427 if (imod == SEMOD_JPLHOR_DEFAULT) imod = 0; break;
4428 case SE_MODEL_JPLHORA_MODE:
4429 if (imod == SEMOD_JPLHORA_DEFAULT) imod = 0; break;
4430 case SE_MODEL_DELTAT:
4431 if (imod == SEMOD_DELTAT_DEFAULT) imod = 0; break;
4432 }
4433 sprintf(samod0 + strlen(samod0), "%d,", imod);
4434 }
4435 /*if (samod != NULL)
4436 strcpy(samod, samod0);*/
4437 *sdet = '\0';
4438 if (sdet != NULL) {
4439 /* JPL ephemeris number and tidal acceleration used with it */
4440 sprintf(sdet + strlen(sdet), "JPL eph. %d; tidal acc. Moon used by SE: %.4f\n",
4441 swi_get_denum(SEI_SUN, iflag), swe_get_tid_acc());
4442 if (iflag & SEFLG_JPLEPH) {
4443 if (iflag & SEFLG_JPLHOR)
4444 strcat(sdet, "JPL Horizons method:\n");
4445 if (iflag & SEFLG_JPLHOR_APPROX)
4446 strcat(sdet, "JPL Horizons method (approximation):\n");
4447 } else if (iflag & SEFLG_SWIEPH) {
4448 strcat(sdet, "Swiss Ephemeris compressed files sepl*/semo*\n");
4449 } else {
4450 strcat(sdet, "Moshier semi-analytical approximation\n");
4451 }
4452 /* long-term Delta T calculation */
4453 get_deltat_model(pmodel[SE_MODEL_DELTAT], s);
4454 sprintf(sdet + strlen(sdet), "Delta T (long-term): %s\n", s);
4455 /* precession model */
4456 get_precession_model(pmodel[SE_MODEL_PREC_LONGTERM], iflag, s);
4457 sprintf(sdet + strlen(sdet), "Precession: %s\n", s);
4458 if (pmodel[SE_MODEL_PREC_LONGTERM] != pmodel[SE_MODEL_PREC_SHORTTERM] && !(iflag & (SEFLG_JPLHOR | SEFLG_JPLHOR_APPROX))) {
4459 get_precession_model(pmodel[SE_MODEL_PREC_SHORTTERM], iflag, s);
4460 sprintf(sdet + strlen(sdet), "+ short-term model: %s\n", s);
4461 }
4462 /* nutation */
4463 get_nutation_model(pmodel[SE_MODEL_NUT], iflag, s);
4464 sprintf(sdet + strlen(sdet), "Nutation: %s\n", s);
4465 /* frame bias */
4466 get_frame_bias_model(pmodel[SE_MODEL_BIAS], s);
4467 sprintf(sdet + strlen(sdet), "Frame bias: %s\n", s);
4468 /* sidereal time */
4469 get_sidt_model(pmodel[SE_MODEL_SIDT], s);
4470 sprintf(sdet + strlen(sdet), "Sid. time: %s\n", s);
4471 /* swetest parameters */
4472 sprintf(sdet + strlen(sdet), "swetest parameters: D P P N B J J S\n");
4473 sprintf(sdet + strlen(sdet), " -amod%s", samod0);
4474 sprintf(sdet + strlen(sdet), " -tidacc%f", swe_get_tid_acc());
4475 strcat(sdet, "\n");
4476 /* list all available astronomical models */
4477 if (!list_all_models) {
4478 sprintf(sdet + strlen(sdet), "For list of all available astronomical models, add a '+' to the version string\n(swetest parameter -amod%s+ or -amod%s+)\n", samod, samod0);
4479 } else {
4480 strcat(sdet, "DELTA T MODELS (D)\n");
4481 for (i = 0; i <= SEMOD_NDELTAT; i++) {
4482 if (i == SEMOD_DELTAT_DEFAULT) continue;
4483 sprintf(sdet + strlen(sdet), " (%d)", i);
4484 if (i == 0) sprintf(sdet + strlen(sdet), " (=%d)", SEMOD_DELTAT_DEFAULT);
4485 get_deltat_model(i, s);
4486 sprintf(sdet + strlen(sdet), ": %s\n", s);
4487 }
4488 strcat(sdet, "PRECESSION MODELS (P P) (long-term/short-term)\n");
4489 for (i = 0; i <= SEMOD_NPREC; i++) {
4490 if (i == SEMOD_PREC_DEFAULT) continue;
4491 sprintf(sdet + strlen(sdet), " (%d)", i);
4492 if (i == 0) sprintf(sdet + strlen(sdet), " (=%d)", SEMOD_PREC_DEFAULT);
4493 get_precession_model(i, iflag, s);
4494 sprintf(sdet + strlen(sdet), ": %s\n", s);
4495 }
4496 strcat(sdet, "NUTATION MODELS (N)\n");
4497 for (i = 0; i <= SEMOD_NNUT; i++) {
4498 if (i == SEMOD_NUT_DEFAULT) continue;
4499 sprintf(sdet + strlen(sdet), " (%d)", i);
4500 if (i == 0) sprintf(sdet + strlen(sdet), " (=%d)", SEMOD_NUT_DEFAULT);
4501 get_nutation_model(i, iflag, s);
4502 sprintf(sdet + strlen(sdet), ": %s\n", s);
4503 }
4504 strcat(sdet, "FRAME BIAS MODELS (B)\n");
4505 for (i = 0; i <= SEMOD_NBIAS; i++) {
4506 if (i == SEMOD_BIAS_DEFAULT) continue;
4507 sprintf(sdet + strlen(sdet), " (%d)", i);
4508 if (i == 0) sprintf(sdet + strlen(sdet), " (=%d)", SEMOD_BIAS_DEFAULT);
4509 get_frame_bias_model(i, s);
4510 sprintf(sdet + strlen(sdet), ": %s\n", s);
4511 }
4512 strcat(sdet, "JPL HORIZONS MODELS (J) (with SEFLG_JPLEPH|SEFLG_JPLHOR).\n");
4513 strcat(sdet, " IAU 1980 (Wahr) + daily corrections to dpsi/deps 1962-today.\n");
4514 strcat(sdet, " (0 (=1): between 1799 and 1962, dpsi/deps of 20-jan-1962 are used.\n");
4515 strcat(sdet, " For times beyond the dpsi/deps table, the last tabulated values are used.\n");
4516 strcat(sdet, " Beyond 1799 and 2201, precession Owen 1990 is used..\n");
4517 strcat(sdet, " Documentation in swephexp.h under 'methods of JPL Horizons'\n");
4518 strcat(sdet, "JPL HORIZONS APPROXIMATION (J) (with SEFLG_JPLEPH|SEFLG_JPLHORA)\n");
4519 strcat(sdet, " Documentation in swephexp.h under 'methods of JPL Horizons'\n");
4520 strcat(sdet, "SIDEREAL TIME MODELS (S)\n");
4521 for (i = 0; i <= SEMOD_NSIDT; i++) {
4522 if (i == SEMOD_SIDT_DEFAULT) continue;
4523 sprintf(sdet + strlen(sdet), " (%d)", i);
4524 if (i == 0) sprintf(sdet + strlen(sdet), " (=%d)", SEMOD_SIDT_DEFAULT);
4525 get_sidt_model(i, s);
4526 sprintf(sdet + strlen(sdet), ": %s\n", s);
4527 }
4528 }
4529 }
4530 }
4531
swi_strcpy(char * to,char * from)4532 char *swi_strcpy(char *to, char *from)
4533 {
4534 char *sp, s[AS_MAXCH];
4535 if (*from == '\0') {
4536 *to = '\0';
4537 return to;
4538 }
4539 if (strlen(from) < AS_MAXCH) {
4540 strcpy(s, from);
4541 strcpy(to, s);
4542 } else {
4543 sp = strdup(from);
4544 if (sp == NULL) {
4545 strcpy(to, from);
4546 } else {
4547 strcpy(to, sp);
4548 free(sp);
4549 }
4550 }
4551 return to;
4552 }
4553
4554 #ifdef TRACE
swi_open_trace(char * serr)4555 void swi_open_trace(char *serr)
4556 {
4557 swi_trace_count++;
4558 if (swi_trace_count >= TRACE_COUNT_MAX) {
4559 if (swi_trace_count == TRACE_COUNT_MAX) {
4560 if (serr != NULL)
4561 sprintf(serr, "trace stopped, %d calls exceeded.", TRACE_COUNT_MAX);
4562 if (swi_fp_trace_out != NULL)
4563 fprintf(swi_fp_trace_out, "trace stopped, %d calls exceeded.\n", TRACE_COUNT_MAX);
4564 if (swi_fp_trace_c != NULL)
4565 fprintf(swi_fp_trace_c, "/* trace stopped, %d calls exceeded. */\n", TRACE_COUNT_MAX);
4566 }
4567 return;
4568 }
4569 if (swi_fp_trace_c == NULL) {
4570 char fname[AS_MAXCH];
4571 #if TRACE == 2
4572 char *sp, *sp1;
4573 int ipid;
4574 #endif
4575 /* remove(fname_trace_c); */
4576 strcpy(fname, fname_trace_c);
4577 #if TRACE == 2
4578 sp = strchr(fname_trace_c, '.');
4579 sp1 = strchr(fname, '.');
4580 # if MSDOS
4581 ipid = _getpid();
4582 # else
4583 ipid = getpid();
4584 # endif
4585 sprintf(sp1, "_%d%s", ipid, sp);
4586 #endif
4587 if ((swi_fp_trace_c = fopen(fname, FILE_A_ACCESS)) == NULL) {
4588 if (serr != NULL) {
4589 sprintf(serr, "could not open trace output file '%s'", fname);
4590 }
4591 } else {
4592 fputs("#include \"sweodef.h\"\n", swi_fp_trace_c);
4593 fputs("#include \"swephexp.h\"\n\n", swi_fp_trace_c);
4594 fputs("void main()\n{\n", swi_fp_trace_c);
4595 fputs(" double tjd, t, nut, eps; int i, ipl, retc; int32 iflag;\n", swi_fp_trace_c);
4596 fputs(" double armc, geolat, cusp[12], ascmc[10]; int hsys;\n", swi_fp_trace_c);
4597 fputs(" double xx[6]; int32 iflgret;\n", swi_fp_trace_c);
4598 fputs(" char s[AS_MAXCH], star[AS_MAXCH], serr[AS_MAXCH];\n", swi_fp_trace_c);
4599 fflush(swi_fp_trace_c);
4600 }
4601 }
4602 if (swi_fp_trace_out == NULL) {
4603 char fname[AS_MAXCH];
4604 #if TRACE == 2
4605 char *sp, *sp1;
4606 int ipid;
4607 #endif
4608 /* remove(fname_trace_out); */
4609 strcpy(fname, fname_trace_out);
4610 #if TRACE == 2
4611 sp = strchr(fname_trace_out, '.');
4612 sp1 = strchr(fname, '.');
4613 # if MSDOS
4614 ipid = _getpid();
4615 # else
4616 ipid = getpid();
4617 # endif
4618 sprintf(sp1, "_%d%s", ipid, sp);
4619 #endif
4620 if ((swi_fp_trace_out = fopen(fname, FILE_A_ACCESS)) == NULL) {
4621 if (serr != NULL) {
4622 sprintf(serr, "could not open trace output file '%s'", fname);
4623 }
4624 }
4625 }
4626 }
4627 #endif
4628