1
2 /*
3 |
4 | Subroutines for reading JPL ephemerides.
5 | derived from testeph.f as contained in DE403 distribution July 1995.
6 | works with DE200, DE102, DE403, DE404, DE405, DE406, DE431
7 | (attention, these ephemerides do not have exactly the same reference frame)
8
9 Authors: Dieter Koch and Alois Treindl, Astrodienst Zurich
10
11 ************************************************************/
12 /* Copyright (C) 1997 - 2021 Astrodienst AG, Switzerland. All rights reserved.
13
14 License conditions
15 ------------------
16
17 This file is part of Swiss Ephemeris.
18
19 Swiss Ephemeris is distributed with NO WARRANTY OF ANY KIND. No author
20 or distributor accepts any responsibility for the consequences of using it,
21 or for whether it serves any particular purpose or works at all, unless he
22 or she says so in writing.
23
24 Swiss Ephemeris is made available by its authors under a dual licensing
25 system. The software developer, who uses any part of Swiss Ephemeris
26 in his or her software, must choose between one of the two license models,
27 which are
28 a) GNU Affero General Public License (AGPL)
29 b) Swiss Ephemeris Professional License
30
31 The choice must be made before the software developer distributes software
32 containing parts of Swiss Ephemeris to others, and before any public
33 service using the developed software is activated.
34
35 If the developer choses the AGPL software license, he or she must fulfill
36 the conditions of that license, which includes the obligation to place his
37 or her whole software project under the AGPL or a compatible license.
38 See https://www.gnu.org/licenses/agpl-3.0.html
39
40 If the developer choses the Swiss Ephemeris Professional license,
41 he must follow the instructions as found in http://www.astro.com/swisseph/
42 and purchase the Swiss Ephemeris Professional Edition from Astrodienst
43 and sign the corresponding license contract.
44
45 The License grants you the right to use, copy, modify and redistribute
46 Swiss Ephemeris, but only under certain conditions described in the License.
47 Among other things, the License requires that the copyright notices and
48 this notice be preserved on all copies.
49
50 Authors of the Swiss Ephemeris: Dieter Koch and Alois Treindl
51
52 The authors of Swiss Ephemeris have no control or influence over any of
53 the derived works, i.e. over software or services created by other
54 programmers which use Swiss Ephemeris functions.
55
56 The names of the authors or of the copyright holder (Astrodienst) must not
57 be used for promoting any software, product or service which uses or contains
58 the Swiss Ephemeris. This copyright notice is the ONLY place where the
59 names of the authors can legally appear, except in cases where they have
60 given special permission in writing.
61
62 The trademarks 'Swiss Ephemeris' and 'Swiss Ephemeris inside' may be used
63 for promoting such software, products or services.
64 */
65
66 #if MSDOS
67 #else
68 #define _FILE_OFFSET_BITS 64
69 #endif
70
71 #include <string.h>
72 #include "swephexp.h"
73 #include "sweph.h"
74 #include "swejpl.h"
75
76 #if MSDOS
77 typedef __int64 off_t64;
78 #define FSEEK _fseeki64
79 #define FTELL _ftelli64
80 #else
81 typedef off_t off_t64;
82 #define FSEEK fseeko
83 #define FTELL ftello
84 #endif
85
86 #define DEBUG_DO_SHOW FALSE
87
88 /*
89 * local globals
90 */
91 struct jpl_save {
92 char *jplfname;
93 char *jplfpath;
94 FILE *jplfptr;
95 short do_reorder;
96 double eh_cval[400];
97 double eh_ss[3], eh_au, eh_emrat;
98 int32 eh_denum, eh_ncon, eh_ipt[39];
99 char ch_cnam[6*400];
100 double pv[78];
101 double pvsun[6];
102 double buf[1500];
103 double pc[18], vc[18], ac[18], jc[18];
104 short do_km;
105 };
106
107 static TLS struct jpl_save *js;
108
109 static int state (double et, int32 *list, int do_bary,
110 double *pv, double *pvsun, double *nut, char *serr);
111 static int interp(double *buf, double t, double intv, int32 ncfin,
112 int32 ncmin, int32 nain, int32 ifl, double *pv);
113 static int32 fsizer(char *serr);
114 static void reorder(char *x, int size, int number);
115 static int read_const_jpl(double *ss, char *serr);
116
117 /* information about eh_ipt[] and buf[]
118 DE200 DE102 DE403
119 3 3 ipt[0] 3 body 0 (mercury) starts at buf[2]
120 12 15 ipt[1] 14 body 0, ncf = coefficients per component
121 4 2 ipt[2] 4 na = nintervals, tot 14*4*3=168
122 147 93 ipt[3] 171 body 1 (venus) starts at buf[170]
123 12 15 ipt[4] 10 ncf = coefficients per component
124 1 1 ipt[5] 2 total 10*2*3=60
125 183 138 ipt[6] 231 body 2 (earth) starts at buf[230]
126 15 15 ipt[7] 13 ncf = coefficients per component
127 2 2 ipt[8] 2 total 13*2*3=78
128 273 228 ipt[9] 309 body 3 (mars) starts at buf[308]
129 10 10 ipt[10] 11 ncf = coefficients per component
130 1 1 ipt[11] 1 total 11*1*3=33
131 303 258 ipt[12] 342 body 4 (jupiter) at buf[341]
132 9 9 ipt[13] 8 total 8 * 1 * 3 = 24
133 1 1 ipt[14] 1
134 330 285 ipt[15] 366 body 5 (saturn) at buf[365]
135 8 8 ipt[16] 7 total 7 * 1 * 3 = 21
136 1 1 ipt[17] 1
137 354 309 ipt[18] 387 body 6 (uranus) at buf[386]
138 8 8 ipt[19] 6 total 6 * 1 * 3 = 18
139 1 1 ipt[20] 1
140 378 333 ipt[21] 405 body 7 (neptune) at buf[404]
141 6 6 ipt[22] 6 total 18
142 1 1 ipt[23] 1
143 396 351 ipt[24] 423 body 8 (pluto) at buf[422]
144 6 6 ipt[25] 6 total 18
145 1 1 ipt[26] 1
146 414 369 ipt[27] 441 body 9 (moon) at buf[440]
147 12 15 ipt[28] 13 total 13 * 8 * 3 = 312
148 8 8 ipt[29] 8
149 702 729 ipt[30] 753 SBARY SUN, starts at buf[752]
150 15 15 ipt[31] 11 SBARY SUN, ncf = coeff per component
151 1 1 ipt[32] 2 total 11*2*3=66
152 747 774 ipt[33] 819 nutations, starts at buf[818]
153 10 0 ipt[34] 10 total 10 * 4 * 2 = 80
154 4 0 ipt[35] 4 (nutation only two coordinates)
155 0 0 ipt[36] 899 librations, start at buf[898]
156 0 0 ipt[37] 10 total 10 * 4 * 3 = 120
157 0 0 ipt[38] 4
158
159 last element of buf[1017]
160 buf[0] contains start jd and buf[1] end jd of segment;
161 each segment is 32 days in de403, 64 days in DE102, 32 days in DE200
162
163 Length of blocks: DE406 = 1456*4=5824 bytes = 728 double
164 DE405 = 2036*4=8144 bytes = 1018 double
165 DE404 = 1456*4=5824 bytes = 728 double
166 DE403 = 2036*4=8144 bytes = 1018 double
167 DE200 = 1652*4=6608 bytes = 826 double
168 DE102 = 1546*4=6184 bytes = 773 double
169 each DE102 record has 53*8=424 fill bytes so that
170 the records have the same length as DE200.
171 */
172
173 /*
174 * This subroutine opens the file jplfname, with a phony record length,
175 * reads the first record, and uses the info to compute ksize,
176 * the number of single precision words in a record.
177 * RETURN: ksize (record size of ephemeris data)
178 * jplfptr is opened on return.
179 * note 26-aug-2008: now record size is computed by fsizer(), not
180 * set to a fixed value depending as in previous releases. The caller of
181 * fsizer() will verify by data comparison whether it computed correctly.
182 */
fsizer(char * serr)183 static int32 fsizer(char *serr)
184 {
185 /* Local variables */
186 int32 ncon;
187 double emrat;
188 int32 numde;
189 double au, ss[3];
190 int i, kmx, khi, nd;
191 int32 ksize, lpt[3];
192 char ttl[6*14*3];
193 size_t nrd; /* unused, removes compile warnings */
194 if ((js->jplfptr = swi_fopen(SEI_FILE_PLANET, js->jplfname, js->jplfpath, serr)) == NULL) {
195 return NOT_AVAILABLE;
196 }
197 /* ttl = ephemeris title, e.g.
198 * "JPL Planetary Ephemeris DE404/LE404
199 * Start Epoch: JED= 625296.5-3001 DEC 21 00:00:00
200 * Final Epoch: JED= 2817168.5 3001 JAN 17 00:00:00c */
201 nrd = fread((void *) &ttl[0], 1, 252, js->jplfptr);
202 if (nrd != 252) return NOT_AVAILABLE;
203 /* cnam = names of constants */
204 nrd = fread((void *) js->ch_cnam, 1, 6*400, js->jplfptr);
205 if (nrd != 6*400) return NOT_AVAILABLE;
206 /* ss[0] = start epoch of ephemeris
207 * ss[1] = end epoch
208 * ss[2] = segment size in days */
209 nrd = fread((void *) &ss[0], sizeof(double), 3, js->jplfptr);
210 if (nrd != 3) return NOT_AVAILABLE;
211 /* reorder ? */
212 if (ss[2] < 1 || ss[2] > 200)
213 js->do_reorder = TRUE;
214 else
215 js->do_reorder = 0;
216 for (i = 0; i < 3; i++)
217 js->eh_ss[i] = ss[i];
218 if (js->do_reorder)
219 reorder((char *) &js->eh_ss[0], sizeof(double), 3);
220 /* plausibility test of these constants. Start and end date must be
221 * between -20000 and +20000, segment size >= 1 and <= 200 */
222 if (js->eh_ss[0] < -5583942 || js->eh_ss[1] > 9025909 || js->eh_ss[2] < 1 || js->eh_ss[2] > 200) {
223 if (serr != NULL) {
224 strcpy(serr, "alleged ephemeris file has invalid format.");
225 if (strlen(serr) + strlen(js->jplfname) + 3 < AS_MAXCH) {
226 sprintf(serr, "alleged ephemeris file (%s) has invalid format.", js->jplfname);
227 }
228 }
229 return(NOT_AVAILABLE);
230 }
231 /* ncon = number of constants */
232 nrd = fread((void *) &ncon, sizeof(int32), 1, js->jplfptr);
233 if (nrd != 1) return NOT_AVAILABLE;
234 if (js->do_reorder)
235 reorder((char *) &ncon, sizeof(int32), 1);
236 /* au = astronomical unit */
237 nrd = fread((void *) &au, sizeof(double), 1, js->jplfptr);
238 if (nrd != 1) return NOT_AVAILABLE;
239 if (js->do_reorder)
240 reorder((char *) &au, sizeof(double), 1);
241 /* emrat = earth moon mass ratio */
242 nrd = fread((void *) &emrat, sizeof(double), 1, js->jplfptr);
243 if (nrd != 1) return NOT_AVAILABLE;
244 if (js->do_reorder)
245 reorder((char *) &emrat, sizeof(double), 1);
246 /* ipt[i+0]: coefficients of planet i start at buf[ipt[i+0]-1]
247 * ipt[i+1]: number of coefficients (interpolation order - 1)
248 * ipt[i+2]: number of intervals in segment */
249 nrd = fread((void *) &js->eh_ipt[0], sizeof(int32), 36, js->jplfptr);
250 if (nrd != 36) return NOT_AVAILABLE;
251 if (js->do_reorder)
252 reorder((char *) &js->eh_ipt[0], sizeof(int32), 36);
253 /* numde = number of jpl ephemeris "404" with de404 */
254 nrd = fread((void *) &numde, sizeof(int32), 1, js->jplfptr);
255 if (nrd != 1) return NOT_AVAILABLE;
256 if (js->do_reorder)
257 reorder((char *) &numde, sizeof(int32), 1);
258 /* read librations */
259 nrd = fread(&lpt[0], sizeof(int32), 3, js->jplfptr);
260 if (nrd != 3) return NOT_AVAILABLE;
261 if (js->do_reorder)
262 reorder((char *) &lpt[0], sizeof(int32), 3);
263 /* fill librations into eh_ipt[36]..[38] */
264 for (i = 0; i < 3; ++i)
265 js->eh_ipt[i + 36] = lpt[i];
266 rewind(js->jplfptr);
267 /* find the number of ephemeris coefficients from the pointers */
268 /* re-activated this code on 26-aug-2008 */
269 kmx = 0;
270 khi = 0;
271 for (i = 0; i < 13; i++) {
272 if (js->eh_ipt[i * 3] > kmx) {
273 kmx = js->eh_ipt[i * 3];
274 khi = i + 1;
275 }
276 }
277 if (khi == 12)
278 nd = 2;
279 else
280 nd = 3;
281 ksize = (js->eh_ipt[khi * 3 - 3] + nd * js->eh_ipt[khi * 3 - 2] * js->eh_ipt[khi * 3 - 1] - 1L) * 2L;
282 /*
283 * de102 files give wrong ksize, because they contain 424 empty bytes
284 * per record. Fixed by hand!
285 */
286 if (ksize == 1546)
287 ksize = 1652;
288 #if 0 /* we prefer to compute ksize to be comaptible
289 with new DE releases */
290 switch (numde) {
291 case 403:
292 case 405:
293 case 410:
294 case 413:
295 case 414:
296 case 418:
297 case 421:
298 ksize = 2036;
299 break;
300 case 404:
301 case 406:
302 ksize = 1456;
303 break;
304 case 200:
305 ksize = 1652;
306 break;
307 case 102:
308 ksize = 1652; /* de102 is filled with blanks to length of de200 */
309 break;
310 default:
311 if (serr != NULL)
312 sprintf(serr,"unknown numde value %d;", numde);
313 return ERR;
314 }
315 #endif
316 if (ksize < 1000 || ksize > 5000) {
317 if (serr != NULL)
318 sprintf(serr, "JPL ephemeris file does not provide valid ksize (%d)", ksize);/**/
319 return NOT_AVAILABLE;
320 }
321 return ksize;
322 }
323
324 /*
325 * This subroutine reads the jpl planetary ephemeris
326 * and gives the position and velocity of the point 'ntarg'
327 * with respect to 'ncent'.
328 * calling sequence parameters:
329 * et = d.p. julian ephemeris date at which interpolation
330 * is wanted.
331 * ** note the entry dpleph for a doubly-dimensioned time **
332 * the reason for this option is discussed in the
333 * subroutine state
334 * ntarg = integer number of 'target' point.
335 * ncent = integer number of center point.
336 * the numbering convention for 'ntarg' and 'ncent' is:
337 * 0 = mercury 7 = neptune
338 * 1 = venus 8 = pluto
339 * 2 = earth 9 = moon
340 * 3 = mars 10 = sun
341 * 4 = jupiter 11 = solar-system barycenter
342 * 5 = saturn 12 = earth-moon barycenter
343 * 6 = uranus 13 = nutations (longitude and obliq)
344 * 14 = librations, if on eph file
345 * (if nutations are wanted, set ntarg = 13. for librations,
346 * set ntarg = 14. set ncent=0.)
347 * rrd = output 6-word d.p. array containing position and velocity
348 * of point 'ntarg' relative to 'ncent'. the units are au and
349 * au/day. for librations the units are radians and radians
350 * per day. in the case of nutations the first four words of
351 * rrd will be set to nutations and rates, having units of
352 * radians and radians/day.
353 * The option is available to have the units in km and km/sec.
354 * For this, set do_km=TRUE (default FALSE).
355 */
swi_pleph(double et,int ntarg,int ncent,double * rrd,char * serr)356 int swi_pleph(double et, int ntarg, int ncent, double *rrd, char *serr)
357 {
358 int i, retc;
359 int32 list[12];
360 double *pv = js->pv;
361 double *pvsun = js->pvsun;
362 for (i = 0; i < 6; ++i)
363 rrd[i] = 0.0;
364 if (ntarg == ncent)
365 return 0;
366 for (i = 0; i < 12; ++i)
367 list[i] = 0;
368 /* check for nutation call */
369 if (ntarg == J_NUT) {
370 if (js->eh_ipt[34] > 0) {
371 list[10] = 2;
372 return(state(et, list, FALSE, pv, pvsun, rrd, serr));
373 } else {
374 if (serr != NULL)
375 sprintf(serr,"No nutations on the JPL ephemeris file;");
376 return (NOT_AVAILABLE);
377 }
378 }
379 if (ntarg == J_LIB) {
380 if (js->eh_ipt[37] > 0) {
381 list[11] = 2;
382 if ((retc = state(et, list, FALSE, pv, pvsun, rrd, serr)) != OK)
383 return (retc);
384 for (i = 0; i < 6; ++i)
385 rrd[i] = pv[i + 60];
386 return 0;
387 } else {
388 if (serr != NULL)
389 sprintf(serr,"No librations on the ephemeris file;");
390 return (NOT_AVAILABLE);
391 }
392 }
393 /* set up proper entries in 'list' array for state call */
394 if (ntarg < J_SUN)
395 list[ntarg] = 2;
396 if (ntarg == J_MOON) /* Mooon needs Earth */
397 list[J_EARTH] = 2;
398 if (ntarg == J_EARTH) /* Earth needs Moon */
399 list[J_MOON] = 2;
400 if (ntarg == J_EMB) /* EMB needs Earth */
401 list[J_EARTH] = 2;
402 if (ncent < J_SUN)
403 list[ncent] = 2;
404 if (ncent == J_MOON) /* Mooon needs Earth */
405 list[J_EARTH] = 2;
406 if (ncent == J_EARTH) /* Earth needs Moon */
407 list[J_MOON] = 2;
408 if (ncent == J_EMB) /* EMB needs Earth */
409 list[J_EARTH] = 2;
410 if ((retc = state(et, list, TRUE, pv, pvsun, rrd, serr)) != OK)
411 return (retc);
412 if (ntarg == J_SUN || ncent == J_SUN) {
413 for (i = 0; i < 6; ++i)
414 pv[i + 6*J_SUN] = pvsun[i];
415 }
416 if (ntarg == J_SBARY || ncent == J_SBARY) {
417 for (i = 0; i < 6; ++i) {
418 pv[i + 6*J_SBARY] = 0.;
419 }
420 }
421 if (ntarg == J_EMB || ncent == J_EMB) {
422 for (i = 0; i < 6; ++i)
423 pv[i + 6*J_EMB] = pv[i + 6*J_EARTH];
424 }
425 if ((ntarg==J_EARTH && ncent==J_MOON) || (ntarg == J_MOON && ncent==J_EARTH)){
426 for (i = 0; i < 6; ++i)
427 pv[i + 6*J_EARTH] = 0.;
428
429 } else {
430 if (list[J_EARTH] == 2) {
431 for (i = 0; i < 6; ++i)
432 pv[i + 6*J_EARTH] -= pv[i + 6*J_MOON] / (js->eh_emrat + 1.);
433 }
434 if (list[J_MOON] == 2) {
435 for (i = 0; i < 6; ++i) {
436 pv[i + 6*J_MOON] += pv[i + 6*J_EARTH];
437 }
438 }
439 }
440 for (i = 0; i < 6; ++i)
441 rrd[i] = pv[i + ntarg * 6] - pv[i + ncent * 6];
442 return OK;
443 }
444
445 /*
446 * This subroutine differentiates and interpolates a
447 * set of chebyshev coefficients to give pos, vel, acc, and jerk
448 * calling sequence parameters:
449 * input:
450 * buf 1st location of array of d.p. chebyshev coefficients of position
451 * t is dp fractional time in interval covered by
452 * coefficients at which interpolation is wanted, 0 <= t <= 1
453 * intv is dp length of whole interval in input time units.
454 * ncf number of coefficients per component
455 * ncm number of components per set of coefficients
456 * na number of sets of coefficients in full array
457 * (i.e., number of sub-intervals in full interval)
458 * ifl int flag: =1 for positions only
459 * =2 for pos and vel
460 * =3 for pos, vel, and acc
461 * =4 for pos, vel, acc, and jerk
462 * output:
463 * pv d.p. interpolated quantities requested.
464 * assumed dimension is pv(ncm,fl).
465 */
interp(double * buf,double t,double intv,int32 ncfin,int32 ncmin,int32 nain,int32 ifl,double * pv)466 static int interp(double *buf, double t, double intv, int32 ncfin,
467 int32 ncmin, int32 nain, int32 ifl, double *pv)
468 {
469 /* Initialized data */
470 static TLS int np, nv;
471 static TLS int nac;
472 static TLS int njk;
473 static TLS double twot = 0.;
474 double *pc = js->pc;
475 double *vc = js->vc;
476 double *ac = js->ac;
477 double *jc = js->jc;
478 int ncf = (int) ncfin;
479 int ncm = (int) ncmin;
480 int na = (int) nain;
481 /* Local variables */
482 double temp;
483 int i, j, ni;
484 double tc;
485 double dt1, bma;
486 double bma2, bma3;
487 /*
488 | get correct sub-interval number for this set of coefficients and then
489 | get normalized chebyshev time within that subinterval.
490 */
491 if (t >= 0)
492 dt1 = floor(t);
493 else
494 dt1 = -floor(-t);
495 temp = na * t;
496 ni = (int) (temp - dt1);
497 /* tc is the normalized chebyshev time (-1 <= tc <= 1) */
498 tc = (fmod(temp, 1.0) + dt1) * 2. - 1.;
499 /*
500 * check to see whether chebyshev time has changed,
501 * and compute new polynomial values if it has.
502 * (the element pc(2) is the value of t1(tc) and hence
503 * contains the value of tc on the previous call.)
504 */
505 if (tc != pc[1]) {
506 np = 2;
507 nv = 3;
508 nac = 4;
509 njk = 5;
510 pc[1] = tc;
511 twot = tc + tc;
512 }
513 /*
514 * be sure that at least 'ncf' polynomials have been evaluated
515 * and are stored in the array 'pc'.
516 */
517 if (np < ncf) {
518 for (i = np; i < ncf; ++i)
519 pc[i] = twot * pc[i - 1] - pc[i - 2];
520 np = ncf;
521 }
522 /* interpolate to get position for each component */
523 for (i = 0; i < ncm; ++i) {
524 pv[i] = 0.;
525 for (j = ncf-1; j >= 0; --j)
526 pv[i] += pc[j] * buf[j + (i + ni * ncm) * ncf];
527 }
528 if (ifl <= 1)
529 return 0;
530 /*
531 * if velocity interpolation is wanted, be sure enough
532 * derivative polynomials have been generated and stored.
533 */
534 bma = (na + na) / intv;
535 vc[2] = twot + twot;
536 if (nv < ncf) {
537 for (i = nv; i < ncf; ++i)
538 vc[i] = twot * vc[i - 1] + pc[i - 1] + pc[i - 1] - vc[i - 2];
539 nv = ncf;
540 }
541 /* interpolate to get velocity for each component */
542 for (i = 0; i < ncm; ++i) {
543 pv[i + ncm] = 0.;
544 for (j = ncf-1; j >= 1; --j)
545 pv[i + ncm] += vc[j] * buf[j + (i + ni * ncm) * ncf];
546 pv[i + ncm] *= bma;
547 }
548 if (ifl == 2)
549 return 0;
550 /* check acceleration polynomial values, and */
551 /* re-do if necessary */
552 bma2 = bma * bma;
553 ac[3] = pc[1] * 24.;
554 if (nac < ncf) {
555 nac = ncf;
556 for (i = nac; i < ncf; ++i)
557 ac[i] = twot * ac[i - 1] + vc[i - 1] * 4. - ac[i - 2];
558 }
559 /* get acceleration for each component */
560 for (i = 0; i < ncm; ++i) {
561 pv[i + ncm * 2] = 0.;
562 for (j = ncf-1; j >= 2; --j)
563 pv[i + ncm * 2] += ac[j] * buf[j + (i + ni * ncm) * ncf];
564 pv[i + ncm * 2] *= bma2;
565 }
566 if (ifl == 3)
567 return 0;
568 /* check jerk polynomial values, and */
569 /* re-do if necessary */
570 bma3 = bma * bma2;
571 jc[4] = pc[1] * 192.;
572 if (njk < ncf) {
573 njk = ncf;
574 for (i = njk; i < ncf; ++i)
575 jc[i] = twot * jc[i - 1] + ac[i - 1] * 6. - jc[i - 2];
576 }
577 /* get jerk for each component */
578 for (i = 0; i < ncm; ++i) {
579 pv[i + ncm * 3] = 0.;
580 for (j = ncf-1; j >= 3; --j)
581 pv[i + ncm * 3] += jc[j] * buf[j + (i + ni * ncm) * ncf];
582 pv[i + ncm * 3] *= bma3;
583 }
584 return 0;
585 }
586
587 /*
588 | ********** state ********************
589 | this subroutine reads and interpolates the jpl planetary ephemeris file
590 | calling sequence parameters:
591 | input:
592 | et dp julian ephemeris epoch at which interpolation is wanted.
593 | list 12-word integer array specifying what interpolation
594 | is wanted for each of the bodies on the file.
595 | list(i)=0, no interpolation for body i
596 | =1, position only
597 | =2, position and velocity
598 | the designation of the astronomical bodies by i is:
599 | i = 0: mercury
600 | = 1: venus
601 | = 2: earth-moon barycenter, NOT earth!
602 | = 3: mars
603 | = 4: jupiter
604 | = 5: saturn
605 | = 6: uranus
606 | = 7: neptune
607 | = 8: pluto
608 | = 9: geocentric moon
609 | =10: nutations in longitude and obliquity
610 | =11: lunar librations (if on file)
611 | If called with list = NULL, only the header records are read and
612 | stored in the global areas.
613 | do_bary short, if true, barycentric, if false, heliocentric.
614 | only the 9 planets 0..8 are affected by it.
615 | output:
616 | pv dp 6 x 11 array that will contain requested interpolated
617 | quantities. the body specified by list(i) will have its
618 | state in the array starting at pv(1,i). (on any given
619 | call, only those words in 'pv' which are affected by the
620 | first 10 'list' entries (and by list(11) if librations are
621 | on the file) are set. the rest of the 'pv' array
622 | is untouched.) the order of components starting in
623 | pv is: x,y,z,dx,dy,dz.
624 | all output vectors are referenced to the earth mean
625 | equator and equinox of epoch. the moon state is always
626 | geocentric; the other nine states are either heliocentric
627 | or solar-system barycentric, depending on the setting of
628 | common flags (see below).
629 | lunar librations, if on file, are put into pv(k,10) if
630 | list(11) is 1 or 2.
631 | pvsun dp 6-word array containing the barycentric position and
632 | velocity of the sun.
633 | nut dp 4-word array that will contain nutations and rates,
634 | depending on the setting of list(10). the order of
635 | quantities in nut is:
636 | d psi (nutation in longitude)
637 | d epsilon (nutation in obliquity)
638 | d psi dot
639 | d epsilon dot
640 | globals used:
641 | do_km logical flag defining physical units of the output states.
642 | TRUE = return km and km/sec, FALSE = return au and au/day
643 | default value = FALSE (km determines time unit
644 | for nutations and librations. angle unit is always radians.)
645 */
state(double et,int32 * list,int do_bary,double * pv,double * pvsun,double * nut,char * serr)646 static int state(double et, int32 *list, int do_bary,
647 double *pv, double *pvsun, double *nut, char *serr)
648 {
649 int i, j, k;
650 int32 nseg;
651 off_t64 flen, nb;
652 double *buf = js->buf;
653 double aufac, s, t, intv, ts[4];
654 int32 nrecl, ksize;
655 int32 nr;
656 double et_mn, et_fr;
657 int32 *ipt = js->eh_ipt;
658 char ch_ttl[252];
659 static TLS int32 irecsz;
660 static TLS int32 nrl, lpt[3], ncoeffs;
661 size_t nrd; /* unused, removes compile warnings */
662 if (js->jplfptr == NULL) {
663 ksize = fsizer(serr); /* the number of single precision words in a record */
664 nrecl = 4;
665 if (ksize == NOT_AVAILABLE)
666 return NOT_AVAILABLE;
667 irecsz = nrecl * ksize; /* record size in bytes */
668 ncoeffs = ksize / 2; /* # of coefficients, doubles */
669 /* ttl = ephemeris title, e.g.
670 * "JPL Planetary Ephemeris DE404/LE404
671 * Start Epoch: JED= 625296.5-3001 DEC 21 00:00:00
672 * Final Epoch: JED= 2817168.5 3001 JAN 17 00:00:00c */
673 nrd = fread((void *) ch_ttl, 1, 252, js->jplfptr);
674 if (nrd != 252) return NOT_AVAILABLE;
675 /* cnam = names of constants */
676 nrd = fread((void *) js->ch_cnam, 1, 2400, js->jplfptr);
677 if (nrd != 2400) return NOT_AVAILABLE;
678 /* ss[0] = start epoch of ephemeris
679 * ss[1] = end epoch
680 * ss[2] = segment size in days */
681 nrd = fread((void *) &js->eh_ss[0], sizeof(double), 3, js->jplfptr);
682 if (nrd != 3) return NOT_AVAILABLE;
683 if (js->do_reorder)
684 reorder((char *) &js->eh_ss[0], sizeof(double), 3);
685 /* ncon = number of constants */
686 nrd = fread((void *) &js->eh_ncon, sizeof(int32), 1, js->jplfptr);
687 if (nrd != 1) return NOT_AVAILABLE;
688 if (js->do_reorder)
689 reorder((char *) &js->eh_ncon, sizeof(int32), 1);
690 /* au = astronomical unit */
691 nrd = fread((void *) &js->eh_au, sizeof(double), 1, js->jplfptr);
692 if (nrd != 1) return NOT_AVAILABLE;
693 if (js->do_reorder)
694 reorder((char *) &js->eh_au, sizeof(double), 1);
695 /* emrat = earth moon mass ratio */
696 nrd = fread((void *) &js->eh_emrat, sizeof(double), 1, js->jplfptr);
697 if (nrd != 1) return NOT_AVAILABLE;
698 if (js->do_reorder)
699 reorder((char *) &js->eh_emrat, sizeof(double), 1);
700 /* ipt[i+0]: coefficients of planet i start at buf[ipt[i+0]-1]
701 * ipt[i+1]: number of coefficients (interpolation order - 1)
702 * ipt[i+2]: number of intervals in segment */
703 nrd = fread((void *) &ipt[0], sizeof(int32), 36, js->jplfptr);
704 if (nrd != 36) return NOT_AVAILABLE;
705 if (js->do_reorder)
706 reorder((char *) &ipt[0], sizeof(int32), 36);
707 /* numde = number of jpl ephemeris "404" with de404 */
708 nrd = fread((void *) &js->eh_denum, sizeof(int32), 1, js->jplfptr);
709 if (nrd != 1) return NOT_AVAILABLE;
710 if (js->do_reorder)
711 reorder((char *) &js->eh_denum, sizeof(int32), 1);
712 nrd = fread((void *) &lpt[0], sizeof(int32), 3, js->jplfptr);
713 if (nrd != 3) return NOT_AVAILABLE;
714 if (js->do_reorder)
715 reorder((char *) &lpt[0], sizeof(int32), 3);
716 /* cval[]: other constants in next record */
717 FSEEK(js->jplfptr, (off_t64) (1L * irecsz), 0);
718 nrd = fread((void *) &js->eh_cval[0], sizeof(double), 400, js->jplfptr);
719 if (nrd != 400) return NOT_AVAILABLE;
720 if (js->do_reorder)
721 reorder((char *) &js->eh_cval[0], sizeof(double), 400);
722 /* new 26-aug-2008: verify correct block size */
723 for (i = 0; i < 3; ++i)
724 ipt[i + 36] = lpt[i];
725 nrl = 0;
726 /* is file length correct? */
727 /* file length */
728 FSEEK(js->jplfptr, (off_t64) 0L, SEEK_END);
729 flen = FTELL(js->jplfptr);
730 /* # of segments in file */
731 nseg = (int32) ((js->eh_ss[1] - js->eh_ss[0]) / js->eh_ss[2]);
732 /* sum of all cheby coeffs of all planets and segments */
733 for(i = 0, nb = 0; i < 13; i++) {
734 k = 3;
735 if (i == 11)
736 k = 2;
737 nb += (ipt[i*3+1] * ipt[i*3+2]) * k * nseg;
738 }
739 /* add start and end epochs of segments */
740 nb += 2 * nseg;
741 /* doubles to bytes */
742 nb *= 8;
743 /* add size of header and constants section */
744 nb += 2 * ksize * nrecl;
745 if (flen != nb
746 /* some of our files are one record too long */
747 && flen - nb != ksize * nrecl
748 ) {
749 if (serr != NULL) {
750 sprintf(serr, "JPL ephemeris file is mutilated; length = %d instead of %d.", (unsigned int) flen, (unsigned int) nb);
751 if (strlen(serr) + strlen(js->jplfname) < AS_MAXCH - 1) {
752 sprintf(serr, "JPL ephemeris file %s is mutilated; length = %d instead of %d.", js->jplfname, (unsigned int) flen, (unsigned int) nb);
753 }
754 }
755 return(NOT_AVAILABLE);
756 }
757 /* check if start and end dates in segments are the same as in
758 * file header */
759 FSEEK(js->jplfptr, (off_t64) (2L * irecsz), 0);
760 nrd = fread((void *) &ts[0], sizeof(double), 2, js->jplfptr);
761 if (nrd != 2) return NOT_AVAILABLE;
762 if (js->do_reorder)
763 reorder((char *) &ts[0], sizeof(double), 2);
764 FSEEK(js->jplfptr, (off_t64) ((nseg + 2 - 1) * ((off_t64) irecsz)), 0);
765 nrd = fread((void *) &ts[2], sizeof(double), 2, js->jplfptr);
766 if (nrd != 2) return NOT_AVAILABLE;
767 if (js->do_reorder)
768 reorder((char *) &ts[2], sizeof(double), 2);
769 if (ts[0] != js->eh_ss[0] || ts[3] != js->eh_ss[1]) {
770 if (serr != NULL)
771 sprintf(serr, "JPL ephemeris file is corrupt; start/end date check failed. %.1f != %.1f || %.1f != %.1f", ts[0],js->eh_ss[0],ts[3],js->eh_ss[1]);
772 return NOT_AVAILABLE;
773 }
774 }
775 if (list == NULL)
776 return 0;
777 s = et - .5;
778 et_mn = floor(s);
779 et_fr = s - et_mn; /* fraction of days since previous midnight */
780 et_mn += .5; /* midnight before epoch */
781 /* error return for epoch out of range */
782 if (et < js->eh_ss[0] || et > js->eh_ss[1]) {
783 if (serr != NULL)
784 sprintf(serr,"jd %f outside JPL eph. range %.2f .. %.2f;", et, js->eh_ss[0], js->eh_ss[1]);
785 return BEYOND_EPH_LIMITS;
786 }
787 /* calculate record # and relative time in interval */
788 nr = (int32) ((et_mn - js->eh_ss[0]) / js->eh_ss[2]) + 2;
789 if (et_mn == js->eh_ss[1])
790 --nr; /* end point of ephemeris, use last record */
791 t = (et_mn - ((nr - 2) * js->eh_ss[2] + js->eh_ss[0]) + et_fr) / js->eh_ss[2];
792 /* read correct record if not in core */
793 if (nr != nrl) {
794 nrl = nr;
795 if (FSEEK(js->jplfptr, (off_t64) (nr * ((off_t64) irecsz)), 0) != 0) {
796 if (serr != NULL)
797 sprintf(serr, "Read error in JPL eph. at %f\n", et);
798 return NOT_AVAILABLE;
799 }
800 for (k = 1; k <= ncoeffs; ++k) {
801 if ( fread((void *) &buf[k - 1], sizeof(double), 1, js->jplfptr) != 1) {
802 if (serr != NULL)
803 sprintf(serr, "Read error in JPL eph. at %f\n", et);
804 return NOT_AVAILABLE;
805 }
806 if (js->do_reorder)
807 reorder((char *) &buf[k-1], sizeof(double), 1);
808 }
809 }
810 if (js->do_km) {
811 intv = js->eh_ss[2] * 86400.;
812 aufac = 1.;
813 } else {
814 intv = js->eh_ss[2];
815 aufac = 1. / js->eh_au;
816 }
817 /* interpolate ssbary sun */
818 interp(&buf[(int) ipt[30] - 1], t, intv, ipt[31], 3L, ipt[32], 2L, pvsun);
819 for (i = 0; i < 6; ++i) {
820 pvsun[i] *= aufac;
821 }
822 /* check and interpolate whichever bodies are requested */
823 for (i = 0; i < 10; ++i) {
824 if (list[i] > 0) {
825 interp(&buf[(int) ipt[i * 3] - 1], t, intv, ipt[i * 3 + 1], 3L,
826 ipt[i * 3 + 2], list[i], &pv[i * 6]);
827 for (j = 0; j < 6; ++j) {
828 if (i < 9 && ! do_bary) {
829 pv[j + i * 6] = pv[j + i * 6] * aufac - pvsun[j];
830 } else {
831 pv[j + i * 6] *= aufac;
832 }
833 }
834 }
835 }
836 /* do nutations if requested (and if on file) */
837 if (list[10] > 0 && ipt[34] > 0) {
838 interp(&buf[(int) ipt[33] - 1], t, intv, ipt[34], 2L, ipt[35],
839 list[10], nut);
840 }
841 /* get librations if requested (and if on file) */
842 if (list[11] > 0 && ipt[37] > 0) {
843 interp(&buf[(int) ipt[36] - 1], t, intv, ipt[37], 3L, ipt[38], list[1],
844 &pv[60]);
845 }
846 return OK;
847 }
848
849 /*
850 * this entry obtains the constants from the ephemeris file
851 * call state to initialize the ephemeris and read in the constants
852 */
read_const_jpl(double * ss,char * serr)853 static int read_const_jpl(double *ss, char *serr)
854 {
855 int i, retc;
856 retc = state(0.0, NULL, FALSE, NULL, NULL, NULL, serr);
857 if (retc != OK)
858 return (retc);
859 for (i = 0; i < 3; i++)
860 ss[i] = js->eh_ss[i];
861 #if DEBUG_DO_SHOW
862 {
863 static const char *bname[] = {
864 "Mercury", "Venus", "EMB", "Mars", "Jupiter", "Saturn",
865 "Uranus", "Neptune", "Pluto", "Moon", "SunBary", "Nut", "Libr"};
866 int j, k;
867 int32 nb, nc;
868 printf(" JPL TEST-EPHEMERIS program. Version October 1995.\n");
869 for (i = 0; i < 13; i++) {
870 j = i * 3;
871 k = 3;
872 if (i == 11) k = 2;
873 nb = js->eh_ipt[j+1] * js->eh_ipt[j+2] * k;
874 nc = (int32) (nb * 36525L / js->eh_ss[2] * 8L);
875 printf("%s\t%d\tipt[%d]\t%3ld %2ld %2ld,\t",
876 bname[i], i, j, js->eh_ipt[j], js->eh_ipt[j+1], js->eh_ipt[j+2]);
877 printf("%3ld double, bytes per century = %6ld\n", nb, nc);
878 fflush(stdout);
879 }
880 printf("%16.2f %16.2f %16.2f\n", js->eh_ss[0], js->eh_ss[1], js->eh_ss[2]);
881 for (i = 0; i < js->eh_ncon; ++i)
882 printf("%.6s\t%24.16f\n", js->ch_cnam + i * 6, js->eh_cval[i]);
883 fflush(stdout);
884 }
885 #endif
886 return OK;
887 }
888
reorder(char * x,int size,int number)889 static void reorder(char *x, int size, int number)
890 {
891 int i, j;
892 char s[8];
893 char *sp1 = x;
894 char *sp2 = &s[0];
895 for (i = 0; i < number; i++) {
896 for (j = 0; j < size; j++)
897 *(sp2 + j) = *(sp1 + size - j - 1);
898 for (j = 0; j < size; j++)
899 *(sp1 + j) = *(sp2 + j);
900 sp1 += size;
901 }
902 }
903
swi_close_jpl_file(void)904 void swi_close_jpl_file(void)
905 {
906 if (js != NULL) {
907 if (js->jplfptr != NULL)
908 fclose(js->jplfptr);
909 if (js->jplfname != NULL)
910 FREE((void *) js->jplfname);
911 if (js->jplfpath != NULL)
912 FREE((void *) js->jplfpath);
913 FREE((void *) js);
914 js = NULL;
915 }
916 }
917
swi_open_jpl_file(double * ss,char * fname,char * fpath,char * serr)918 int swi_open_jpl_file(double *ss, char *fname, char *fpath, char *serr)
919 {
920 int retc = OK;
921 /* if open, return */
922 if (js != NULL && js->jplfptr != NULL)
923 return OK;
924 if ((js = (struct jpl_save *) CALLOC(1, sizeof(struct jpl_save))) == NULL
925 || (js->jplfname = (char *) MALLOC(strlen(fname)+1)) == NULL
926 || (js->jplfpath = (char *) MALLOC(strlen(fpath)+1)) == NULL
927 ) {
928 if (serr != NULL)
929 strcpy(serr, "error in malloc() with JPL ephemeris.");
930 return ERR;
931 }
932 strcpy(js->jplfname, fname);
933 strcpy(js->jplfpath, fpath);
934 retc = read_const_jpl(ss, serr);
935 if (retc != OK)
936 swi_close_jpl_file();
937 else {
938 /* intializations for function interpol() */
939 js->pc[0] = 1;
940 js->pc[1] = 2;
941 js->vc[1] = 1;
942 js->ac[2] = 4;
943 js->jc[3] = 24;
944 }
945 return retc;
946 }
947
swi_get_jpl_denum()948 int32 swi_get_jpl_denum()
949 {
950 return js->eh_denum;
951 }
952
953