1 /* SWISSEPH
2 Moshier planet routines
3
4 modified for SWISSEPH by Dieter Koch
5
6 **************************************************************/
7 /* Copyright (C) 1997 - 2021 Astrodienst AG, Switzerland. All rights reserved.
8
9 License conditions
10 ------------------
11
12 This file is part of Swiss Ephemeris.
13
14 Swiss Ephemeris is distributed with NO WARRANTY OF ANY KIND. No author
15 or distributor accepts any responsibility for the consequences of using it,
16 or for whether it serves any particular purpose or works at all, unless he
17 or she says so in writing.
18
19 Swiss Ephemeris is made available by its authors under a dual licensing
20 system. The software developer, who uses any part of Swiss Ephemeris
21 in his or her software, must choose between one of the two license models,
22 which are
23 a) GNU Affero General Public License (AGPL)
24 b) Swiss Ephemeris Professional License
25
26 The choice must be made before the software developer distributes software
27 containing parts of Swiss Ephemeris to others, and before any public
28 service using the developed software is activated.
29
30 If the developer choses the AGPL software license, he or she must fulfill
31 the conditions of that license, which includes the obligation to place his
32 or her whole software project under the AGPL or a compatible license.
33 See https://www.gnu.org/licenses/agpl-3.0.html
34
35 If the developer choses the Swiss Ephemeris Professional license,
36 he must follow the instructions as found in http://www.astro.com/swisseph/
37 and purchase the Swiss Ephemeris Professional Edition from Astrodienst
38 and sign the corresponding license contract.
39
40 The License grants you the right to use, copy, modify and redistribute
41 Swiss Ephemeris, but only under certain conditions described in the License.
42 Among other things, the License requires that the copyright notices and
43 this notice be preserved on all copies.
44
45 Authors of the Swiss Ephemeris: Dieter Koch and Alois Treindl
46
47 The authors of Swiss Ephemeris have no control or influence over any of
48 the derived works, i.e. over software or services created by other
49 programmers which use Swiss Ephemeris functions.
50
51 The names of the authors or of the copyright holder (Astrodienst) must not
52 be used for promoting any software, product or service which uses or contains
53 the Swiss Ephemeris. This copyright notice is the ONLY place where the
54 names of the authors can legally appear, except in cases where they have
55 given special permission in writing.
56
57 The trademarks 'Swiss Ephemeris' and 'Swiss Ephemeris inside' may be used
58 for promoting such software, products or services.
59 */
60
61 #include <string.h>
62 #include "swephexp.h"
63 #include "sweph.h"
64 #include "swephlib.h"
65 #include "swemptab.h"
66
67 #define TIMESCALE 3652500.0
68
69 #define mods3600(x) ((x) - 1.296e6 * floor ((x)/1.296e6))
70
71 #define FICT_GEO 1
72 #define KGAUSS_GEO 0.0000298122353216 /* Earth only */
73 /* #define KGAUSS_GEO 0.00002999502129737 Earth + Moon */
74
75 static void embofs_mosh(double J, double *xemb);
76 static int check_t_terms(double t, char *sinp, double *doutp);
77
78 static int read_elements_file(int32 ipl, double tjd,
79 double *tjd0, double *tequ,
80 double *mano, double *sema, double *ecce,
81 double *parg, double *node, double *incl,
82 char *pname, int32 *fict_ifl, char *serr);
83
84 static const int pnoint2msh[] = {2, 2, 0, 1, 3, 4, 5, 6, 7, 8, };
85
86
87 /* From Simon et al (1994) */
88 static const double freqs[] =
89 {
90 /* Arc sec per 10000 Julian years. */
91 53810162868.8982,
92 21066413643.3548,
93 12959774228.3429,
94 6890507749.3988,
95 1092566037.7991,
96 439960985.5372,
97 154248119.3933,
98 78655032.0744,
99 52272245.1795
100 };
101
102 static const double phases[] =
103 {
104 /* Arc sec. */
105 252.25090552 * 3600.,
106 181.97980085 * 3600.,
107 100.46645683 * 3600.,
108 355.43299958 * 3600.,
109 34.35151874 * 3600.,
110 50.07744430 * 3600.,
111 314.05500511 * 3600.,
112 304.34866548 * 3600.,
113 860492.1546,
114 };
115
116 static const struct plantbl *planets[] =
117 {
118 &mer404,
119 &ven404,
120 &ear404,
121 &mar404,
122 &jup404,
123 &sat404,
124 &ura404,
125 &nep404,
126 &plu404
127 };
128
129 static TLS double ss[9][24];
130 static TLS double cc[9][24];
131
132 static void sscc (int k, double arg, int n);
133
swi_moshplan2(double J,int iplm,double * pobj)134 int swi_moshplan2 (double J, int iplm, double *pobj)
135 {
136 int i, j, k, m, k1, ip, np, nt;
137 signed char *p;
138 double *pl, *pb, *pr;
139 double su, cu, sv, cv, T;
140 double t, sl, sb, sr;
141 const struct plantbl *plan = planets[iplm];
142
143 T = (J - J2000) / TIMESCALE;
144 /* Calculate sin( i*MM ), etc. for needed multiple angles. */
145 for (i = 0; i < 9; i++)
146 {
147 if ((j = plan->max_harmonic[i]) > 0)
148 {
149 sr = (mods3600 (freqs[i] * T) + phases[i]) * STR;
150 sscc (i, sr, j);
151 }
152 }
153
154 /* Point to start of table of arguments. */
155 p = plan->arg_tbl;
156 /* Point to tabulated cosine and sine amplitudes. */
157 pl = plan->lon_tbl;
158 pb = plan->lat_tbl;
159 pr = plan->rad_tbl;
160 sl = 0.0;
161 sb = 0.0;
162 sr = 0.0;
163
164 for (;;)
165 {
166 /* argument of sine and cosine */
167 /* Number of periodic arguments. */
168 np = *p++;
169 if (np < 0)
170 break;
171 if (np == 0)
172 { /* It is a polynomial term. */
173 nt = *p++;
174 /* Longitude polynomial. */
175 cu = *pl++;
176 for (ip = 0; ip < nt; ip++)
177 {
178 cu = cu * T + *pl++;
179 }
180 sl += mods3600 (cu);
181 /* Latitude polynomial. */
182 cu = *pb++;
183 for (ip = 0; ip < nt; ip++)
184 {
185 cu = cu * T + *pb++;
186 }
187 sb += cu;
188 /* Radius polynomial. */
189 cu = *pr++;
190 for (ip = 0; ip < nt; ip++)
191 {
192 cu = cu * T + *pr++;
193 }
194 sr += cu;
195 continue;
196 }
197 k1 = 0;
198 cv = 0.0;
199 sv = 0.0;
200 for (ip = 0; ip < np; ip++)
201 {
202 /* What harmonic. */
203 j = *p++;
204 /* Which planet. */
205 m = *p++ - 1;
206 if (j)
207 {
208 k = j;
209 if (j < 0)
210 k = -k;
211 k -= 1;
212 su = ss[m][k]; /* sin(k*angle) */
213 if (j < 0)
214 su = -su;
215 cu = cc[m][k];
216 if (k1 == 0)
217 { /* set first angle */
218 sv = su;
219 cv = cu;
220 k1 = 1;
221 }
222 else
223 { /* combine angles */
224 t = su * cv + cu * sv;
225 cv = cu * cv - su * sv;
226 sv = t;
227 }
228 }
229 }
230 /* Highest power of T. */
231 nt = *p++;
232 /* Longitude. */
233 cu = *pl++;
234 su = *pl++;
235 for (ip = 0; ip < nt; ip++)
236 {
237 cu = cu * T + *pl++;
238 su = su * T + *pl++;
239 }
240 sl += cu * cv + su * sv;
241 /* Latitiude. */
242 cu = *pb++;
243 su = *pb++;
244 for (ip = 0; ip < nt; ip++)
245 {
246 cu = cu * T + *pb++;
247 su = su * T + *pb++;
248 }
249 sb += cu * cv + su * sv;
250 /* Radius. */
251 cu = *pr++;
252 su = *pr++;
253 for (ip = 0; ip < nt; ip++)
254 {
255 cu = cu * T + *pr++;
256 su = su * T + *pr++;
257 }
258 sr += cu * cv + su * sv;
259 }
260 pobj[0] = STR * sl;
261 pobj[1] = STR * sb;
262 pobj[2] = STR * plan->distance * sr + plan->distance;
263 return OK;
264 }
265
266 /* Moshier ephemeris.
267 * computes heliocentric cartesian equatorial coordinates of
268 * equinox 2000
269 * for earth and a planet
270 * tjd julian day
271 * ipli internal SWEPH planet number
272 * xp array of 6 doubles for planet's position and speed
273 * xe earth's
274 * serr error string
275 */
swi_moshplan(double tjd,int ipli,AS_BOOL do_save,double * xpret,double * xeret,char * serr)276 int swi_moshplan(double tjd, int ipli, AS_BOOL do_save, double *xpret, double *xeret, char *serr)
277 {
278 int i;
279 int do_earth = FALSE;
280 double dx[3], x2[3], xxe[6], xxp[6];
281 double *xp, *xe;
282 double dt;
283 char s[AS_MAXCH];
284 int iplm = pnoint2msh[ipli];
285 struct plan_data *pdp = &swed.pldat[ipli];
286 struct plan_data *pedp = &swed.pldat[SEI_EARTH];
287 double seps2000 = swed.oec2000.seps;
288 double ceps2000 = swed.oec2000.ceps;
289 if (do_save) {
290 xp = pdp->x;
291 xe = pedp->x;
292 } else {
293 xp = xxp;
294 xe = xxe;
295 }
296 if (do_save || ipli == SEI_EARTH || xeret != NULL)
297 do_earth = TRUE;
298 /* tjd beyond ephemeris limits, give some margin for spped at edge */
299 if (tjd < MOSHPLEPH_START - 0.3 || tjd > MOSHPLEPH_END + 0.3) {
300 if (serr != NULL) {
301 sprintf(s, "jd %f outside Moshier planet range %.2f .. %.2f ",
302 tjd, MOSHPLEPH_START, MOSHPLEPH_END);
303 if (strlen(serr) + strlen(s) < AS_MAXCH)
304 strcat(serr, s);
305 }
306 return(ERR);
307 }
308 /* earth, for geocentric position */
309 if (do_earth) {
310 if (tjd == pedp->teval
311 && pedp->iephe == SEFLG_MOSEPH) {
312 xe = pedp->x;
313 } else {
314 /* emb */
315 swi_moshplan2(tjd, pnoint2msh[SEI_EMB], xe); /* emb hel. ecl. 2000 polar */
316 swi_polcart(xe, xe); /* to cartesian */
317 swi_coortrf2(xe, xe, -seps2000, ceps2000);/* and equator 2000 */
318 embofs_mosh(tjd, xe); /* emb -> earth */
319 if (do_save) {
320 pedp->teval = tjd;
321 pedp->xflgs = -1;
322 pedp->iephe = SEFLG_MOSEPH;
323 }
324 /* one more position for speed. */
325 swi_moshplan2(tjd - PLAN_SPEED_INTV, pnoint2msh[SEI_EMB], x2);
326 swi_polcart(x2, x2);
327 swi_coortrf2(x2, x2, -seps2000, ceps2000);
328 embofs_mosh(tjd - PLAN_SPEED_INTV, x2);/**/
329 for (i = 0; i <= 2; i++)
330 dx[i] = (xe[i] - x2[i]) / PLAN_SPEED_INTV;
331 /* store speed */
332 for (i = 0; i <= 2; i++) {
333 xe[i+3] = dx[i];
334 }
335 }
336 if (xeret != NULL)
337 for (i = 0; i <= 5; i++)
338 xeret[i] = xe[i];
339 }
340 /* earth is the planet wanted */
341 if (ipli == SEI_EARTH) {
342 xp = xe;
343 } else {
344 /* other planet */
345 /* if planet has already been computed, return */
346 if (tjd == pdp->teval && pdp->iephe == SEFLG_MOSEPH) {
347 xp = pdp->x;
348 } else {
349 swi_moshplan2(tjd, iplm, xp);
350 swi_polcart(xp, xp);
351 swi_coortrf2(xp, xp, -seps2000, ceps2000);
352 if (do_save) {
353 pdp->teval = tjd;/**/
354 pdp->xflgs = -1;
355 pdp->iephe = SEFLG_MOSEPH;
356 }
357 /* one more position for speed.
358 * the following dt gives good speed for light-time correction
359 */
360 #if 0
361 for (i = 0; i <= 2; i++)
362 dx[i] = xp[i] - pedp->x[i];
363 dt = LIGHTTIME_AUNIT * sqrt(square_sum(dx));
364 #endif
365 dt = PLAN_SPEED_INTV;
366 swi_moshplan2(tjd - dt, iplm, x2);
367 swi_polcart(x2, x2);
368 swi_coortrf2(x2, x2, -seps2000, ceps2000);
369 for (i = 0; i <= 2; i++)
370 dx[i] = (xp[i] - x2[i]) / dt;
371 /* store speed */
372 for (i = 0; i <= 2; i++) {
373 xp[i+3] = dx[i];
374 }
375 }
376 if (xpret != NULL)
377 for (i = 0; i <= 5; i++)
378 xpret[i] = xp[i];
379 }
380 return(OK);
381 }
382
383
384 /* Prepare lookup table of sin and cos ( i*Lj )
385 * for required multiple angles
386 */
sscc(int k,double arg,int n)387 static void sscc (int k, double arg, int n)
388 {
389 double cu, su, cv, sv, s;
390 int i;
391
392 su = sin (arg);
393 cu = cos (arg);
394 ss[k][0] = su; /* sin(L) */
395 cc[k][0] = cu; /* cos(L) */
396 sv = 2.0 * su * cu;
397 cv = cu * cu - su * su;
398 ss[k][1] = sv; /* sin(2L) */
399 cc[k][1] = cv;
400 for (i = 2; i < n; i++)
401 {
402 s = su * cv + cu * sv;
403 cv = cu * cv - su * sv;
404 sv = s;
405 ss[k][i] = sv; /* sin( i+1 L ) */
406 cc[k][i] = cv;
407 }
408 }
409
410
411 /* Adjust position from Earth-Moon barycenter to Earth
412 *
413 * J = Julian day number
414 * xemb = rectangular equatorial coordinates of Earth
415 */
embofs_mosh(double tjd,double * xemb)416 static void embofs_mosh(double tjd, double *xemb)
417 {
418 double T, M, a, L, B, p;
419 double smp, cmp, s2mp, c2mp, s2d, c2d, sf, cf;
420 double s2f, sx, cx, xyz[6];
421 double seps = swed.oec.seps;
422 double ceps = swed.oec.ceps;
423 int i;
424 /* Short series for position of the Moon
425 */
426 T = (tjd-J1900)/36525.0;
427 /* Mean anomaly of moon (MP) */
428 a = swe_degnorm(((1.44e-5*T + 0.009192)*T + 477198.8491)*T + 296.104608);
429 a *= DEGTORAD;
430 smp = sin(a);
431 cmp = cos(a);
432 s2mp = 2.0*smp*cmp; /* sin(2MP) */
433 c2mp = cmp*cmp - smp*smp; /* cos(2MP) */
434 /* Mean elongation of moon (D) */
435 a = swe_degnorm(((1.9e-6*T - 0.001436)*T + 445267.1142)*T + 350.737486);
436 a = 2.0 * DEGTORAD * a;
437 s2d = sin(a);
438 c2d = cos(a);
439 /* Mean distance of moon from its ascending node (F) */
440 a = swe_degnorm((( -3.e-7*T - 0.003211)*T + 483202.0251)*T + 11.250889);
441 a *= DEGTORAD;
442 sf = sin(a);
443 cf = cos(a);
444 s2f = 2.0*sf*cf; /* sin(2F) */
445 sx = s2d*cmp - c2d*smp; /* sin(2D - MP) */
446 cx = c2d*cmp + s2d*smp; /* cos(2D - MP) */
447 /* Mean longitude of moon (LP) */
448 L = ((1.9e-6*T - 0.001133)*T + 481267.8831)*T + 270.434164;
449 /* Mean anomaly of sun (M) */
450 M = swe_degnorm((( -3.3e-6*T - 1.50e-4)*T + 35999.0498)*T + 358.475833);
451 /* Ecliptic longitude of the moon */
452 L = L
453 + 6.288750*smp
454 + 1.274018*sx
455 + 0.658309*s2d
456 + 0.213616*s2mp
457 - 0.185596*sin( DEGTORAD * M )
458 - 0.114336*s2f;
459 /* Ecliptic latitude of the moon */
460 a = smp*cf;
461 sx = cmp*sf;
462 B = 5.128189*sf
463 + 0.280606*(a+sx) /* sin(MP+F) */
464 + 0.277693*(a-sx) /* sin(MP-F) */
465 + 0.173238*(s2d*cf - c2d*sf); /* sin(2D-F) */
466 B *= DEGTORAD;
467 /* Parallax of the moon */
468 p = 0.950724
469 +0.051818*cmp
470 +0.009531*cx
471 +0.007843*c2d
472 +0.002824*c2mp;
473 p *= DEGTORAD;
474 /* Elongation of Moon from Sun
475 */
476 L = swe_degnorm(L);
477 L *= DEGTORAD;
478 /* Distance in au */
479 a = 4.263523e-5/sin(p);
480 /* Convert to rectangular ecliptic coordinates */
481 xyz[0] = L;
482 xyz[1] = B;
483 xyz[2] = a;
484 swi_polcart(xyz, xyz);
485 /* Convert to equatorial */
486 swi_coortrf2(xyz, xyz, -seps, ceps);
487 /* Precess to equinox of J2000.0 */
488 swi_precess(xyz, tjd, 0, J_TO_J2000);/**/
489 /* now emb -> earth */
490 for (i = 0; i <= 2; i++)
491 xemb[i] -= xyz[i] / (EARTH_MOON_MRAT + 1.0);
492 }
493
494 /* orbital elements of planets that are computed from osculating elements
495 * epoch
496 * equinox
497 * mean anomaly,
498 * semi axis,
499 * eccentricity,
500 * argument of perihelion,
501 * ascending node
502 * inclination
503 */
504 #define SE_NEELY /* use James Neely's revised elements
505 * of Uranian planets*/
506 static const char *plan_fict_nam[SE_NFICT_ELEM] =
507 {"Cupido", "Hades", "Zeus", "Kronos",
508 "Apollon", "Admetos", "Vulkanus", "Poseidon",
509 "Isis-Transpluto", "Nibiru", "Harrington",
510 "Leverrier", "Adams",
511 "Lowell", "Pickering",};
512
swi_get_fict_name(int32 ipl,char * snam)513 char *swi_get_fict_name(int32 ipl, char *snam)
514 {
515 if (read_elements_file(ipl, 0, NULL, NULL,
516 NULL, NULL, NULL, NULL, NULL, NULL,
517 snam, NULL, NULL) == ERR)
518 strcpy(snam, "name not found");
519 return snam;
520 }
521
522 static const double plan_oscu_elem[SE_NFICT_ELEM][8] = {
523 #ifdef SE_NEELY
524 {J1900, J1900, 163.7409, 40.99837, 0.00460, 171.4333, 129.8325, 1.0833},/* Cupido Neely */
525 {J1900, J1900, 27.6496, 50.66744, 0.00245, 148.1796, 161.3339, 1.0500},/* Hades Neely */
526 {J1900, J1900, 165.1232, 59.21436, 0.00120, 299.0440, 0.0000, 0.0000},/* Zeus Neely */
527 {J1900, J1900, 169.0193, 64.81960, 0.00305, 208.8801, 0.0000, 0.0000},/* Kronos Neely */
528 {J1900, J1900, 138.0533, 70.29949, 0.00000, 0.0000, 0.0000, 0.0000},/* Apollon Neely */
529 {J1900, J1900, 351.3350, 73.62765, 0.00000, 0.0000, 0.0000, 0.0000},/* Admetos Neely */
530 {J1900, J1900, 55.8983, 77.25568, 0.00000, 0.0000, 0.0000, 0.0000},/* Vulcanus Neely */
531 {J1900, J1900, 165.5163, 83.66907, 0.00000, 0.0000, 0.0000, 0.0000},/* Poseidon Neely */
532 #else
533 {J1900, J1900, 104.5959, 40.99837, 0, 0, 0, 0}, /* Cupido */
534 {J1900, J1900, 337.4517, 50.667443, 0, 0, 0, 0}, /* Hades */
535 {J1900, J1900, 104.0904, 59.214362, 0, 0, 0, 0}, /* Zeus */
536 {J1900, J1900, 17.7346, 64.816896, 0, 0, 0, 0}, /* Kronos */
537 {J1900, J1900, 138.0354, 70.361652, 0, 0, 0, 0}, /* Apollon */
538 {J1900, J1900, -8.678, 73.736476, 0, 0, 0, 0}, /* Admetos */
539 {J1900, J1900, 55.9826, 77.445895, 0, 0, 0, 0}, /* Vulkanus */
540 {J1900, J1900, 165.3595, 83.493733, 0, 0, 0, 0}, /* Poseidon */
541 #endif
542 /* Isis-Transpluto; elements from "Die Sterne" 3/1952, p. 70ff.
543 * Strubell does not give an equinox. 1945 is taken to best reproduce
544 * ASTRON ephemeris. (This is a strange choice, though.)
545 * The epoch is 1772.76. The year is understood to have 366 days.
546 * The fraction is counted from 1 Jan. 1772 */
547 {2368547.66, 2431456.5, 0.0, 77.775, 0.3, 0.7, 0, 0},
548 /* Nibiru, elements from Christian Woeltge, Hannover */
549 {1856113.380954, 1856113.380954, 0.0, 234.8921, 0.981092, 103.966, -44.567, 158.708},
550 /* Harrington, elements from Astronomical Journal 96(4), Oct. 1988 */
551 {2374696.5, J2000, 0.0, 101.2, 0.411, 208.5, 275.4, 32.4},
552 /* Leverrier's Neptune,
553 according to W.G. Hoyt, "Planets X and Pluto", Tucson 1980, p. 63 */
554 {2395662.5, 2395662.5, 34.05, 36.15, 0.10761, 284.75, 0, 0},
555 /* Adam's Neptune */
556 {2395662.5, 2395662.5, 24.28, 37.25, 0.12062, 299.11, 0, 0},
557 /* Lowell's Pluto */
558 {2425977.5, 2425977.5, 281, 43.0, 0.202, 204.9, 0, 0},
559 /* Pickering's Pluto */
560 {2425977.5, 2425977.5, 48.95, 55.1, 0.31, 280.1, 100, 15}, /**/
561 #if 0 /* Ceres JPL 1600, without perturbations from other minor planets,
562 * from following initial elements:
563 * 2450600.5 2000 0 1 164.7073602 73.0340746 80.5995101
564 * 10.5840296 0.07652422 0.0 2.770176095 */
565 {2305447.5, J2000, 0.5874558977449977e+02, 0.2766536058742327e+01,
566 0.7870946565779195e-01, 0.5809199028919189e+02,
567 0.8650119410725021e+02, 0.1066835622280712e+02},
568 /* Chiron, Bowell database 18-mar-1997 */
569 {2450500.5, J2000, 7.258191, 13.67387471, 0.38174778, 339.558345, 209.379239, 6.933360}, /**/
570 #endif
571 };
572
573 /* computes a planet from osculating elements *
574 * tjd julian day
575 * ipl body number
576 * ipli body number in planetary data structure
577 * iflag flags
578 */
swi_osc_el_plan(double tjd,double * xp,int ipl,int ipli,double * xearth,double * xsun,char * serr)579 int swi_osc_el_plan(double tjd, double *xp, int ipl, int ipli, double *xearth, double *xsun, char *serr)
580 {
581 double pqr[9], x[6];
582 double eps, K, fac, rho, cose, sine;
583 double alpha, beta, zeta, sigma, M2, Msgn, M_180_or_0;
584 double tjd0, tequ, mano, sema, ecce, parg, node, incl, dmot;
585 double cosnode, sinnode, cosincl, sinincl, cosparg, sinparg;
586 double M, E;
587 struct plan_data *pedp = &swed.pldat[SEI_EARTH];
588 struct plan_data *pdp = &swed.pldat[ipli];
589 int32 fict_ifl = 0;
590 int i;
591 /* orbital elements, either from file or, if file not found,
592 * from above built-in set
593 */
594 if (read_elements_file(ipl, tjd, &tjd0, &tequ,
595 &mano, &sema, &ecce, &parg, &node, &incl,
596 NULL, &fict_ifl, serr) == ERR)
597 return ERR;
598 dmot = 0.9856076686 * DEGTORAD / sema / sqrt(sema); /* daily motion */
599 if (fict_ifl & FICT_GEO)
600 dmot /= sqrt(SUN_EARTH_MRAT);
601 cosnode = cos(node);
602 sinnode = sin(node);
603 cosincl = cos(incl);
604 sinincl = sin(incl);
605 cosparg = cos(parg);
606 sinparg = sin(parg);
607 /* Gaussian vector */
608 pqr[0] = cosparg * cosnode - sinparg * cosincl * sinnode;
609 pqr[1] = -sinparg * cosnode - cosparg * cosincl * sinnode;
610 pqr[2] = sinincl * sinnode;
611 pqr[3] = cosparg * sinnode + sinparg * cosincl * cosnode;
612 pqr[4] = -sinparg * sinnode + cosparg * cosincl * cosnode;
613 pqr[5] = -sinincl * cosnode;
614 pqr[6] = sinparg * sinincl;
615 pqr[7] = cosparg * sinincl;
616 pqr[8] = cosincl;
617 /* Kepler problem */
618 E = M = swi_mod2PI(mano + (tjd - tjd0) * dmot); /* mean anomaly of date */
619 /* better E for very high eccentricity and small M */
620 if (ecce > 0.975) {
621 M2 = M * RADTODEG;
622 if (M2 > 150 && M2 < 210) {
623 M2 -= 180;
624 M_180_or_0 = 180;
625 } else
626 M_180_or_0 = 0;
627 if (M2 > 330)
628 M2 -= 360;
629 if (M2 < 0) {
630 M2 = -M2;
631 Msgn = -1;
632 } else
633 Msgn = 1;
634 if (M2 < 30) {
635 M2 *= DEGTORAD;
636 alpha = (1 - ecce) / (4 * ecce + 0.5);
637 beta = M2 / (8 * ecce + 1);
638 zeta = pow(beta + sqrt(beta * beta + alpha * alpha), 1/3);
639 sigma = zeta - alpha / 2;
640 sigma = sigma - 0.078 * sigma * sigma * sigma * sigma * sigma / (1 + ecce);
641 E = Msgn * (M2 + ecce * (3 * sigma - 4 * sigma * sigma * sigma))
642 + M_180_or_0;
643 }
644 }
645 E = swi_kepler(E, M, ecce);
646 /* position and speed, referred to orbital plane */
647 if (fict_ifl & FICT_GEO)
648 K = KGAUSS_GEO / sqrt(sema);
649 else
650 K = KGAUSS / sqrt(sema);
651 cose = cos(E);
652 sine = sin(E);
653 fac = sqrt((1 - ecce) * (1 + ecce));
654 rho = 1 - ecce * cose;
655 x[0] = sema * (cose - ecce);
656 x[1] = sema * fac * sine;
657 x[3] = -K * sine / rho;
658 x[4] = K * fac * cose / rho;
659 /* transformation to ecliptic */
660 xp[0] = pqr[0] * x[0] + pqr[1] * x[1];
661 xp[1] = pqr[3] * x[0] + pqr[4] * x[1];
662 xp[2] = pqr[6] * x[0] + pqr[7] * x[1];
663 xp[3] = pqr[0] * x[3] + pqr[1] * x[4];
664 xp[4] = pqr[3] * x[3] + pqr[4] * x[4];
665 xp[5] = pqr[6] * x[3] + pqr[7] * x[4];
666 /* transformation to equator */
667 eps = swi_epsiln(tequ, 0);
668 swi_coortrf(xp, xp, -eps);
669 swi_coortrf(xp+3, xp+3, -eps);
670 /* precess to J2000 */
671 if (tequ != J2000) {
672 swi_precess(xp, tequ, 0, J_TO_J2000);
673 swi_precess(xp+3, tequ, 0, J_TO_J2000);
674 }
675 /* to solar system barycentre */
676 if (fict_ifl & FICT_GEO) {
677 for (i = 0; i <= 5; i++) {
678 xp[i] += xearth[i];
679 }
680 } else {
681 for (i = 0; i <= 5; i++) {
682 xp[i] += xsun[i];
683 }
684 }
685 if (pdp->x == xp) {
686 pdp->teval = tjd; /* for precession! */
687 pdp->iephe = pedp->iephe;
688 }
689 return OK;
690 }
691
692 #if 1
693 /* note: input parameter tjd is required for T terms in elements */
read_elements_file(int32 ipl,double tjd,double * tjd0,double * tequ,double * mano,double * sema,double * ecce,double * parg,double * node,double * incl,char * pname,int32 * fict_ifl,char * serr)694 static int read_elements_file(int32 ipl, double tjd,
695 double *tjd0, double *tequ,
696 double *mano, double *sema, double *ecce,
697 double *parg, double *node, double *incl,
698 char *pname, int32 *fict_ifl, char *serr)
699 {
700 int i, iline, iplan, retc, ncpos;
701 FILE *fp = NULL;
702 char s[AS_MAXCH], *sp;
703 char *cpos[20], serri[AS_MAXCH];
704 AS_BOOL elem_found = FALSE;
705 double tt = 0;
706 /* -1, because file information is not saved, file is always closed */
707 if ((fp = swi_fopen(-1, SE_FICTFILE, swed.ephepath, serr)) == NULL) {
708 /* file does not exist, use built-in bodies */
709 if (ipl >= SE_NFICT_ELEM) {
710 if (serr != NULL)
711 sprintf(serr, "error no elements for fictitious body no %7.0f", (double) ipl);
712 return ERR;
713 }
714 if (tjd0 != NULL)
715 *tjd0 = plan_oscu_elem[ipl][0]; /* epoch */
716 if (tequ != NULL)
717 *tequ = plan_oscu_elem[ipl][1]; /* equinox */
718 if (mano != NULL)
719 *mano = plan_oscu_elem[ipl][2] * DEGTORAD; /* mean anomaly */
720 if (sema != NULL)
721 *sema = plan_oscu_elem[ipl][3]; /* semi-axis */
722 if (ecce != NULL)
723 *ecce = plan_oscu_elem[ipl][4]; /* eccentricity */
724 if (parg != NULL)
725 *parg = plan_oscu_elem[ipl][5] * DEGTORAD; /* arg. of peri. */
726 if (node != NULL)
727 *node = plan_oscu_elem[ipl][6] * DEGTORAD; /* asc. node */
728 if (incl != NULL)
729 *incl = plan_oscu_elem[ipl][7] * DEGTORAD; /* inclination */
730 if (pname != NULL)
731 strcpy(pname, plan_fict_nam[ipl]);
732 return OK;
733 }
734 /*
735 * find elements in file
736 */
737 iline = 0;
738 iplan = -1;
739 while (fgets(s, AS_MAXCH, fp) != NULL) {
740 iline++;
741 sp = s;
742 while(*sp == ' ' || *sp == '\t')
743 sp++;
744 swi_strcpy(s, sp);
745 if (*s == '#')
746 continue;
747 if (*s == '\r')
748 continue;
749 if (*s == '\n')
750 continue;
751 if (*s == '\0')
752 continue;
753 if ((sp = strchr(s, '#')) != NULL)
754 *sp = '\0';
755 ncpos = swi_cutstr(s, ",", cpos, 20);
756 sprintf(serri, "error in file %s, line %7.0f:", SE_FICTFILE, (double) iline);
757 if (ncpos < 9) {
758 if (serr != NULL) {
759 sprintf(serr, "%s nine elements required", serri);
760 }
761 goto return_err;
762 }
763 iplan++;
764 if (iplan != ipl)
765 continue;
766 elem_found = TRUE;
767 /* epoch of elements */
768 if (tjd0 != NULL) {
769 sp = cpos[0];
770 for (i = 0; i < 5; i++)
771 sp[i] = tolower(sp[i]);
772 if (strncmp(sp, "j2000", 5) == OK)
773 *tjd0 = J2000;
774 else if (strncmp(sp, "b1950", 5) == OK)
775 *tjd0 = B1950;
776 else if (strncmp(sp, "j1900", 5) == OK)
777 *tjd0 = J1900;
778 else if (*sp == 'j' || *sp == 'b') {
779 if (serr != NULL) {
780 sprintf(serr, "%s invalid epoch", serri);
781 }
782 goto return_err;
783 } else
784 *tjd0 = atof(sp);
785 tt = tjd - *tjd0;
786 }
787 /* equinox */
788 if (tequ != NULL) {
789 sp = cpos[1];
790 while(*sp == ' ' || *sp == '\t')
791 sp++;
792 for (i = 0; i < 5; i++)
793 sp[i] = tolower(sp[i]);
794 if (strncmp(sp, "j2000", 5) == OK)
795 *tequ = J2000;
796 else if (strncmp(sp, "b1950", 5) == OK)
797 *tequ = B1950;
798 else if (strncmp(sp, "j1900", 5) == OK)
799 *tequ = J1900;
800 else if (strncmp(sp, "jdate", 5) == OK)
801 *tequ = tjd;
802 else if (*sp == 'j' || *sp == 'b') {
803 if (serr != NULL) {
804 sprintf(serr, "%s invalid equinox", serri);
805 }
806 goto return_err;
807 } else
808 *tequ = atof(sp);
809 }
810 /* mean anomaly t0 */
811 if (mano != NULL) {
812 retc = check_t_terms(tt, cpos[2], mano);
813 *mano = swe_degnorm(*mano);
814 if (retc == ERR) {
815 if (serr != NULL) {
816 sprintf(serr, "%s mean anomaly value invalid", serri);
817 }
818 goto return_err;
819 }
820 /* if mean anomaly has t terms (which happens with fictitious
821 * planet Vulcan), we set
822 * epoch = tjd, so that no motion will be added anymore
823 * equinox = tjd */
824 if (retc == 1) {
825 *tjd0 = tjd;
826 }
827 *mano *= DEGTORAD;
828 }
829 /* semi-axis */
830 if (sema != NULL) {
831 retc = check_t_terms(tt, cpos[3], sema);
832 if (*sema <= 0 || retc == ERR) {
833 if (serr != NULL) {
834 sprintf(serr, "%s semi-axis value invalid", serri);
835 }
836 goto return_err;
837 }
838 }
839 /* eccentricity */
840 if (ecce != NULL) {
841 retc = check_t_terms(tt, cpos[4], ecce);
842 if (*ecce >= 1 || *ecce < 0 || retc == ERR) {
843 if (serr != NULL) {
844 sprintf(serr, "%s eccentricity invalid (no parabolic or hyperbolic orbits allowed)", serri);
845 }
846 goto return_err;
847 }
848 }
849 /* perihelion argument */
850 if (parg != NULL) {
851 retc = check_t_terms(tt, cpos[5], parg);
852 *parg = swe_degnorm(*parg);
853 if (retc == ERR) {
854 if (serr != NULL) {
855 sprintf(serr, "%s perihelion argument value invalid", serri);
856 }
857 goto return_err;
858 }
859 *parg *= DEGTORAD;
860 }
861 /* node */
862 if (node != NULL) {
863 retc = check_t_terms(tt, cpos[6], node);
864 *node = swe_degnorm(*node);
865 if (retc == ERR) {
866 if (serr != NULL) {
867 sprintf(serr, "%s node value invalid", serri);
868 }
869 goto return_err;
870 }
871 *node *= DEGTORAD;
872 }
873 /* inclination */
874 if (incl != NULL) {
875 retc = check_t_terms(tt, cpos[7], incl);
876 *incl = swe_degnorm(*incl);
877 if (retc == ERR) {
878 if (serr != NULL) {
879 sprintf(serr, "%s inclination value invalid", serri);
880 }
881 goto return_err;
882 }
883 *incl *= DEGTORAD;
884 }
885 /* planet name */
886 if (pname != NULL) {
887 sp = cpos[8];
888 while(*sp == ' ' || *sp == '\t')
889 sp++;
890 swi_right_trim(sp);
891 strcpy(pname, sp);
892 }
893 /* geocentric */
894 if (fict_ifl != NULL && ncpos > 9) {
895 for (sp = cpos[9]; *sp != '\0'; sp++)
896 *sp = tolower(*sp);
897 if (strstr(cpos[9], "geo") != NULL)
898 *fict_ifl |= FICT_GEO;
899 }
900 break;
901 }
902 if (!elem_found) {
903 if (serr != NULL) {
904 sprintf(serr, "%s elements for planet %7.0f not found", serri, (double) ipl);
905 }
906 goto return_err;
907 }
908 fclose(fp);
909 return OK;
910 return_err:
911 fclose(fp);
912 return ERR;
913 }
914 #endif
915
check_t_terms(double t,char * sinp,double * doutp)916 static int check_t_terms(double t, char *sinp, double *doutp)
917 {
918 int i, isgn = 1, z;
919 int retc = 0;
920 char *sp;
921 double tt[5], fac;
922 tt[0] = t / 36525;
923 tt[1] = tt[0];
924 tt[2] = tt[1] * tt[1];
925 tt[3] = tt[2] * tt[1];
926 tt[4] = tt[3] * tt[1];
927 if ((sp = strpbrk(sinp, "+-")) != NULL)
928 retc = 1; /* with additional terms */
929 sp = sinp;
930 *doutp = 0;
931 fac = 1;
932 z = 0;
933 while (1) {
934 while(*sp != '\0' && strchr(" \t", *sp) != NULL)
935 sp++;
936 if (strchr("+-", *sp) || *sp == '\0') {
937 if (z > 0)
938 *doutp += fac;
939 isgn = 1;
940 if (*sp == '-')
941 isgn = -1;
942 fac = 1 * isgn;
943 if (*sp == '\0')
944 return retc;
945 sp++;
946 } else {
947 while(*sp != '\0' && strchr("* \t", *sp) != NULL)
948 sp++;
949 if (*sp != '\0' && strchr("tT", *sp) != NULL) {
950 /* a T */
951 sp++;
952 if (*sp != '\0' && strchr("+-", *sp))
953 fac *= tt[0];
954 else if ((i = atoi(sp)) <= 4 && i >= 0)
955 fac *= tt[i];
956 } else {
957 /* a number */
958 if (atof(sp) != 0 || *sp == '0')
959 fac *= atof(sp);
960 }
961 while (*sp != '\0' && strchr("0123456789.", *sp))
962 sp++;
963 }
964 z++;
965 }
966 return retc; /* there have been additional terms */
967 }
968