1 /* jpleph.cpp: JPL ephemeris functions
2
3 Copyright (C) 2011, Project Pluto
4
5 This program is free software; you can redistribute it and/or
6 modify it under the terms of the GNU General Public License
7 as published by the Free Software Foundation; either version 2
8 of the License, or (at your option) any later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program; if not, write to the Free Software
17 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
18 02110-1301, USA. */
19
20 /*****************************************************************************
21 * ***** jpl planetary and lunar ephemerides ***** C ver.1.2 *
22 ******************************************************************************
23 * *
24 * This program was written in standard fortran-77 and it was manually *
25 * translated to C language by Piotr A. Dybczynski (dybol@phys.amu.edu.pl), *
26 * subsequently revised heavily by Bill J Gray (pluto@gwi.net), just short *
27 * of a total re-write. *
28 * *
29 ******************************************************************************
30 * Last modified: July 23, 1997 by PAD *
31 ******************************************************************************
32 21 Apr 2010: Revised by Bill J. Gray. The code now determines the kernel
33 size, then allocates memory accordingly. This should 'future-proof' us in
34 case JPL (or someone else) creates kernels that are larger than the previously
35 arbitrary MAX_KERNEL_SIZE parameter. 'swap_long' and 'swap_double' have
36 been replaced with 'swap_32_bit_val' and 'swap_64_bit_val'. It also works
37 on 64-bit compiles now.
38
39 16 Mar 2001: Revised by Bill J. Gray. You can now use binary
40 ephemerides with either byte order ('big-endian' or 'small-endian');
41 the code checks to see if the data is in the "wrong" order for the
42 current platform, and swaps bytes on-the-fly if needed. (Yes, this
43 can result in a slowdown... sometimes as much as 1%. The function is
44 so mathematically intensive that the byte-swapping is the least of our
45 troubles.) You can also use DE-200, 403, 404, 405, or 406 without
46 recompiling (the constan() function now determines which ephemeris is
47 in use and its byte order).
48
49 Also, I did some minor optimization of the interp() (Chebyshev
50 interpolation) function, resulting in a bit of a speedup.
51
52 The code has been modified to be a separately linkable component, with
53 details of the implementation encapsulated.
54 *****************************************************************************/
55 // Make Large File Support work also on ARM boards, explicitly.
56 #define _FILE_OFFSET_BITS 64
57
58 #include <stdio.h>
59 #include <assert.h>
60 #include <math.h>
61 #include <string.h>
62 #include <stdlib.h>
63 #include <stdint.h>
64
65 #include "StelUtils.hpp"
66 /**** include variable and type definitions, specific for this C version */
67
68 #include "jpleph.h"
69 #include "jpl_int.h"
70
71 #define TRUE 1
72 #define FALSE 0
73
74
75 // GZ patches for Large File Support for DE431 past AD10100...
76 #if defined(Q_OS_WIN)
77 #define FSeek(__FILE, __OFFSET, _MODE) _fseeki64(__FILE, __OFFSET, _MODE)
78 #else
79 #define FSeek(__FILE, __OFFSET, _MODE) fseeko(__FILE, __OFFSET, _MODE)
80 #endif
81
82
jpl_get_double(const void * ephem,const int value)83 double DLL_FUNC jpl_get_double(const void *ephem, const int value)
84 {
85 return(*reinterpret_cast<const double *>(static_cast<const char *>(ephem) + value));
86 }
87
jpl_get_long(const void * ephem,const int value)88 long DLL_FUNC jpl_get_long(const void *ephem, const int value)
89 {
90 return(*reinterpret_cast<const int32_t *>(static_cast<const char *>(ephem) + value));
91 }
92
93
94 /*****************************************************************************
95 ** jpl_pleph(ephem,et,ntar,ncent,rrd,calc_velocity) **
96 ******************************************************************************
97 ** **
98 ** This subroutine reads the jpl planetary ephemeris **
99 ** and gives the position and velocity of the point 'ntarg' **
100 ** with respect to 'ncent'. **
101 ** **
102 ** Calling sequence parameters: **
103 ** **
104 ** et = (double) julian ephemeris date at which interpolation **
105 ** is wanted. **
106 ** **
107 ** ntarg = integer number of 'target' point. **
108 ** **
109 ** ncent = integer number of center point. **
110 ** **
111 ** The numbering convention for 'ntarg' and 'ncent' is: **
112 ** **
113 ** 1 = mercury 8 = neptune **
114 ** 2 = venus 9 = pluto **
115 ** 3 = earth 10 = moon **
116 ** 4 = mars 11 = sun **
117 ** 5 = jupiter 12 = solar-system barycenter **
118 ** 6 = saturn 13 = earth-moon barycenter **
119 ** 7 = uranus 14 = nutations (longitude and obliq) **
120 ** 15 = librations, if on eph. file **
121 ** 16 = lunar mantle omega_x,omega_y,omega_z**
122 ** 17 = TT-TDB, if on eph. file **
123 ** **
124 ** (If nutations are wanted, set ntarg = 14. **
125 ** For librations, set ntarg = 15. set ncent= 0. **
126 ** For TT-TDB, set ntarg = 17. I've not actually **
127 ** seen an ntarg = 16 case yet.) **
128 ** **
129 ** rrd = output 6-element, double array of position and velocity **
130 ** of point 'ntarg' relative to 'ncent'. The units are au and **
131 ** au/day. For librations the units are radians and radians **
132 ** per day. In the case of nutations the first four words of **
133 ** rrd will be set to nutations and rates, having units of **
134 ** radians and radians/day. **
135 ** **
136 ** The option is available to have the units in km and km/sec. **
137 ** for this, set km=TRUE at the beginning of the program. **
138 ** **
139 ** calc_velocity = integer flag; if nonzero, velocities will be **
140 ** computed, otherwise not. **
141 ** **
142 *****************************************************************************/
jpl_pleph(void * ephem,const double et,const int ntarg,const int ncent,double rrd[],const int calc_velocity)143 int DLL_FUNC jpl_pleph(void *ephem, const double et, const int ntarg,
144 const int ncent, double rrd[], const int calc_velocity)
145 {
146 struct jpl_eph_data *eph = static_cast<struct jpl_eph_data *>(ephem);
147 double pv[13][6]={{0.}};/* pv is the position/velocity array
148 NUMBERED FROM ZERO: 0=Mercury,1=Venus,...
149 8=Pluto,9=Moon,10=Sun,11=SSBary,12=EMBary
150 First 10 elements (0-9) are affected by
151 jpl_state(), all are adjusted here. */
152
153
154 int rval = 0;
155 const int list_val = (calc_velocity ? 2 : 1);
156 unsigned i;
157 int list[14]; /* list is a vector denoting, for which "body"
158 ephemeris values should be calculated by
159 jpl_state(): 0=Mercury,1=Venus,2=EMBary,...,
160 8=Pluto, 9=geocentric Moon, 10=nutations in
161 long. & obliq. 11= lunar librations;
162 12 = TT-TDB, 13=lunar mantle omegas */
163
164 for(i = 0; i < 6; ++i) rrd[i] = 0.0;
165
166 if(ntarg == ncent) return(0);
167
168 for(i = 0; i < sizeof(list) / sizeof(list[0]); i++)
169 list[i] = 0;
170
171 /* Because of the whacko indexing in JPL ephemerides, we need */
172 /* to work out way through the following indexing schemes : */
173 /* ntarg ipt list */
174 /* 14 11 10 Nutations */
175 /* 15 12 11 Librations */
176 /* 16 13 12 Lunar mantle angular vel */
177 /* 17 14 13 TT - TDB */
178
179 for(i = 0; i < 4; i++)
180 if(ntarg == static_cast<int>(i) + 14)
181 {
182 if(eph->ipt[i + 11][1] > 0) /* quantity is in ephemeris */
183 {
184 list[i + 10] = list_val;
185
186 // GZ: Coverity Scan (travis) automatic test chokes on next line with:
187 // CID 134920: Memory - corruptions (OVERRUN)
188 // Overrunning array "list" of 56 bytes by passing it to a function
189 // which accesses it at byte offset 56.
190 // I see it does explicitly NOT access list[14].
191 // TODO: check again after next round of travis.
192 rval = jpl_state(ephem, et, list, pv, rrd, 0);
193 }
194 else /* quantity doesn't exist in the ephemeris file */
195 rval = JPL_EPH_QUANTITY_NOT_IN_EPHEMERIS;
196 return(rval);
197 }
198 if(ntarg > 13 || ncent > 13 || ntarg < 1 || ncent < 1)
199 return(JPL_EPH_INVALID_INDEX);
200
201 /* force barycentric output by 'state' */
202
203 /* set up proper entries in 'list' array for state call */
204
205 for(i = 0; i < 2; i++) /* list[] IS NUMBERED FROM ZERO ! */
206 {
207 const unsigned k = static_cast<unsigned int>(i ? ncent : ntarg) - 1;
208
209 if(k <= 9) list[k] = list_val; /* Major planets */
210 if(k == 9) list[2] = list_val; /* for moon, earth state is needed */
211 if(k == 2) list[9] = list_val; /* for earth, moon state is needed */
212 if(k == 12) list[2] = list_val; /* EMBary state additionally */
213 }
214
215 /* make call to state */
216 rval = jpl_state(eph, et, list, pv, rrd, 1);
217 /* Solar System barycentric Sun state goes to pv[10][] */
218 if(ntarg == 11 || ncent == 11)
219 for(i = 0; i < 6; i++)
220 pv[10][i] = eph->pvsun[i];
221
222 /* Solar System Barycenter coordinates & velocities equal to zero */
223 if(ntarg == 12 || ncent == 12)
224 for(i = 0; i < 6; i++)
225 pv[11][i] = 0.0;
226
227 /* Solar System barycentric EMBary state: */
228 if(ntarg == 13 || ncent == 13)
229 for(i = 0; i < 6; i++)
230 pv[12][i] = pv[2][i];
231
232 /* if moon from earth or earth from moon ..... */
233 if((ntarg*ncent) == 30 && (ntarg+ncent) == 13)
234 for(i = 0; i < 6; ++i) pv[2][i]=0.0;
235 else
236 {
237 if(list[2]) /* calculate earth state from EMBary */
238 for(i = 0; i < static_cast<unsigned int>(list[2]) * 3u; ++i)
239 pv[2][i] -= pv[9][i]/(1.0+eph->emrat);
240
241 if(list[9]) /* calculate Solar System barycentric moon state */
242 for(i = 0; i < static_cast<unsigned int>(list[9]) * 3u; ++i)
243 pv[9][i] += pv[2][i];
244 }
245
246 for(i = 0; i < static_cast<unsigned int>(list_val) * 3u; ++i)
247 rrd[i] = pv[ntarg-1][i] - pv[ncent-1][i];
248
249 return(rval);
250 }
251
252 /* Some notes about the information stored in 'iinfo': the posn_coeff[]
253 array contains the Chebyshev polynomials for tc,
254
255 posn_coeff[i]=T (tc).
256 i
257
258 The vel_coeff[] array contains the derivatives of the same polynomials,
259
260 vel_coeff[i]=T'(tc).
261 i
262
263 Evaluating these polynomials is a little expensive, and we don't want
264 to evaluate any more than we have to. (Some planets require many more
265 Chebyshev polynomials than others.) So if 'tc' is unchanged, we can
266 rest assured that 'n_posn_avail' Chebyshev polynomials, and 'n_vel_avail'
267 derivatives of Chebyshev polynomials, have already been evaluated, and
268 we start from there, using the recurrence formulae
269
270 T (x) = 1 T (x) = x T (x) = 2xT (x) - T (x)
271 0 1 n+1 n n-1
272
273 T'(x) = 0 T'(x) = 1 T' (x) = 2xT'(x) + 2T (x) - T' (x)
274 0 1 n+1 n n n-1
275
276 (the second set being just the derivatives of the first). To get the
277 _acceleration_ of an object, we just keep going and get the second
278 derivatives as
279
280 T"(x) = 0 T"(x) = 1 T" (x) = 2xT"(x) + 4T'(x) - T" (x)
281 0 1 n+1 n n n-1
282
283 At present, i can range from 0 to 17. If future JPL ephems require
284 Chebyshev polynomials beyond T , those arrays may need to be expanded.
285 17 */
286
287 /*****************************************************************************
288 ** interp(buf,t,ncf,ncm,na,ifl,pv) **
289 ******************************************************************************
290 ** **
291 ** this subroutine differentiates and interpolates a **
292 ** set of chebyshev coefficients to give position and velocity **
293 ** **
294 ** calling sequence parameters: **
295 ** **
296 ** input: **
297 ** **
298 ** iinfo stores certain chunks of interpolation info, in hopes **
299 ** that if you call again with similar parameters, the **
300 ** function won't have to re-compute all coefficients/data. **
301 ** **
302 ** coef 1st location of array of d.p. chebyshev coefficients **
303 ** of position **
304 ** **
305 ** t t[0] is double fractional time in interval covered by **
306 ** coefficients at which interpolation is wanted **
307 ** (0 <= t[0] <= 1). t[1] is dp length of whole **
308 ** interval in input time units. **
309 ** **
310 ** ncf # of coefficients per component **
311 ** **
312 ** ncm # of components per set of coefficients **
313 ** **
314 ** na # of sets of coefficients in full array **
315 ** (i.e., # of sub-intervals in full interval) **
316 ** **
317 ** ifl integer flag: =1 for positions only **
318 ** =2 for pos and vel **
319 ** =3 for pos, vel, accel (currently used for **
320 ** pvsun only) **
321 ** **
322 ** output: **
323 ** **
324 ** posvel interpolated quantities requested. dimension **
325 ** expected is posvel[ncm*ifl], double precision. **
326 ** **
327 *****************************************************************************/
interp(struct interpolation_info * iinfo,const double coef[],const double t[2],const unsigned ncf,const unsigned ncm,const unsigned na,const int velocity_flag,double posvel[])328 static void interp(struct interpolation_info *iinfo,
329 const double coef[], const double t[2], const unsigned ncf, const unsigned ncm,
330 const unsigned na, const int velocity_flag, double posvel[])
331 {
332 const double dna = static_cast<double>(na);
333 const double temp = dna * t[0];
334 unsigned l = static_cast<unsigned>(temp);
335 double vfac, unused_temp1;
336 double tc = 2.0 * modf(temp, &unused_temp1) - 1.0;
337 unsigned i, j;
338
339 assert(ncf < MAX_CHEBY);
340 if(l == na)
341 {
342 l--;
343 tc = 1.;
344 }
345 assert(tc >= -1.);
346 assert(tc <= 1.);
347
348 /* check to see whether chebyshev time has changed, and compute new
349 polynomial values if it has.
350 (the element iinfo->posn_coeff[1] is the value of t1[tc] and hence
351 contains the value of tc on the previous call.) */
352
353
354 if(tc != iinfo->posn_coeff[1])
355 {
356 iinfo->n_posn_avail = 2;
357 iinfo->n_vel_avail = 2;
358 iinfo->posn_coeff[1] = tc;
359 iinfo->twot = tc+tc;
360 }
361
362 /* be sure that at least 'ncf' polynomials have been evaluated and are
363 stored in the array 'iinfo->posn_coeff'. Note that we start out with
364 posn_coeff[0] = 1. and posn_coeff[1] = tc (see 'jpl_init_ephemeris'
365 below), and vel_coeff[0] and [1] are similarly preset. We do that
366 because you need the first two coeffs of those series to start the
367 Chebyshev recurrence; see the comments above this function. */
368
369 if(iinfo->n_posn_avail < ncf)
370 {
371 double *pc_ptr = iinfo->posn_coeff + iinfo->n_posn_avail;
372
373 for(i=ncf - iinfo->n_posn_avail; i; i--, pc_ptr++)
374 *pc_ptr = iinfo->twot * pc_ptr[-1] - pc_ptr[-2];
375 iinfo->n_posn_avail=ncf;
376 }
377
378 /* interpolate to get position for each component */
379
380 for(i = 0; i < ncm; ++i) /* ncm is a number of coordinates */
381 {
382 const double *coeff_ptr = coef + ncf * (i + l * ncm + 1);
383 const double *pc_ptr = iinfo->posn_coeff + ncf;
384
385 *posvel = 0.0;
386 for(j = ncf; j; j--)
387 *posvel += (*--pc_ptr) * (*--coeff_ptr);
388 posvel++;
389 }
390
391 if(velocity_flag <= 1) return;
392
393 /* if velocity interpolation is wanted, be sure enough
394 derivative polynomials have been generated and stored. */
395
396 if(iinfo->n_vel_avail < ncf)
397 {
398 double *vc_ptr = iinfo->vel_coeff + iinfo->n_vel_avail;
399 const double *pc_ptr = iinfo->posn_coeff + iinfo->n_vel_avail - 1;
400
401 for(i = ncf - iinfo->n_vel_avail; i; i--, vc_ptr++, pc_ptr++)
402 *vc_ptr = iinfo->twot * vc_ptr[-1] + *pc_ptr + *pc_ptr - vc_ptr[-2];
403 iinfo->n_vel_avail = ncf;
404 }
405
406 /* interpolate to get velocity for each component */
407
408 vfac = (dna + dna) / t[1];
409 for(i = 0; i < ncm; ++i)
410 {
411 double tval = 0.;
412 const double *coeff_ptr = coef + ncf * (i + l * ncm + 1);
413 const double *vc_ptr = iinfo->vel_coeff + ncf;
414
415 for(j = ncf - 1; j; j--)
416 tval += (*--vc_ptr) * (*--coeff_ptr);
417 *posvel++ = tval * vfac;
418 }
419
420 /* Accelerations are rarely computed -- at present, only */
421 /* for pvsun -- so we don't get so tricky in optimizing. */
422 /* The accel_coeffs (the second derivative of the Chebyshev */
423 /* polynomials) are not stored for repeated use, for example. */
424 if(velocity_flag == 3)
425 {
426 double accel_coeffs[MAX_CHEBY];
427
428 accel_coeffs[0] = accel_coeffs[1] = 0.;
429 for(i = 2; i < ncf; i++) /* recurrence for T"(x) */
430 accel_coeffs[i] = 4. * iinfo->vel_coeff[i - 1]
431 + iinfo->twot * accel_coeffs[i - 1]
432 - accel_coeffs[i - 2];
433
434 for(i = 0; i < ncm; ++i) /* ncm is a number of coordinates */
435 {
436 double tval = 0.;
437 const double *coeff_ptr = coef + ncf * (i + l * ncm + 1);
438 const double *ac_ptr = accel_coeffs + ncf;
439
440 for(j = ncf; j; j--)
441 tval += (*--ac_ptr) * (*--coeff_ptr);
442 *posvel++ = tval * vfac * vfac;
443 }
444 }
445
446 return;
447 }
448
449 /* swap_32_bit_val() and swap_64_bit_val() are used when reading a binary
450 ephemeris that was created on a machine with 'opposite' byte order to
451 the currently-used machine (signalled by the 'swap_bytes' flag in the
452 jpl_eph_data structure). In such cases, every double and integer
453 value read from the ephemeris must be byte-swapped by these two functions. */
454
455 #define SWAP_MACRO(A, B, TEMP) { TEMP = A; A = B; B = TEMP; }
456
swap_32_bit_val(void * ptr)457 static void swap_32_bit_val(void *ptr)
458 {
459 char *tptr = static_cast<char *>(ptr), tchar;
460
461 SWAP_MACRO(tptr[0], tptr[3], tchar);
462 SWAP_MACRO(tptr[1], tptr[2], tchar);
463 }
464
swap_64_bit_val(void * ptr,long count)465 static void swap_64_bit_val(void *ptr, long count)
466 {
467 char *tptr = static_cast<char *>(ptr), tchar;
468
469 while(count--)
470 {
471 SWAP_MACRO(tptr[0], tptr[7], tchar);
472 SWAP_MACRO(tptr[1], tptr[6], tchar);
473 SWAP_MACRO(tptr[2], tptr[5], tchar);
474 SWAP_MACRO(tptr[3], tptr[4], tchar);
475
476 tptr += 8;
477 }
478 }
479
480 /* Most ephemeris quantities have a dimension of three. Planet positions
481 have an x, y, and z; librations and lunar mantle angles have three Euler
482 angles. But TDT-TT is a single quantity, and nutation is expressed as
483 two angles. */
484
dimension(const unsigned int idx)485 static unsigned int dimension(const unsigned int idx)
486 {
487 unsigned int rval;
488
489 if(idx == 11) /* Nutations */
490 rval = 2;
491 else if(idx == 14) /* TDT - TT */
492 rval = 1;
493 else /* planets, lunar mantle angles, librations */
494 rval = 3;
495 return(rval);
496 }
497
498 /*****************************************************************************
499 ** jpl_state(ephem,et2,list,pv,nut,bary) **
500 ******************************************************************************
501 ** This subroutine reads and interpolates the jpl planetary ephemeris file **
502 ** **
503 ** Calling sequence parameters: **
504 ** **
505 ** Input: **
506 ** **
507 ** et2[] double, 2-element JED epoch at which interpolation **
508 ** is wanted. Any combination of et2[0]+et2[1] which falls **
509 ** within the time span on the file is a permissible epoch. **
510 ** **
511 ** a. for ease in programming, the user may put the **
512 ** entire epoch in et2[0] and set et2[1]=0.0 **
513 ** **
514 ** b. for maximum interpolation accuracy, set et2[0] = **
515 ** the most recent midnight at or before interpolation **
516 ** epoch and set et2[1] = fractional part of a day **
517 ** elapsed between et2[0] and epoch. **
518 ** **
519 ** c. as an alternative, it may prove convenient to set **
520 ** et2[0] = some fixed epoch, such as start of integration,**
521 ** and et2[1] = elapsed interval between then and epoch. **
522 ** **
523 ** list 13-element integer array specifying what interpolation **
524 ** is wanted for each of the "bodies" on the file. **
525 ** **
526 ** list[i]=0, no interpolation for body i **
527 ** =1, position only **
528 ** =2, position and velocity **
529 ** **
530 ** the designation of the astronomical bodies by i is: **
531 ** **
532 ** i = 0: mercury **
533 ** = 1: venus **
534 ** = 2: earth-moon barycenter **
535 ** = 3: mars **
536 ** = 4: jupiter **
537 ** = 5: saturn **
538 ** = 6: uranus **
539 ** = 7: neptune **
540 ** = 8: pluto **
541 ** = 9: geocentric moon **
542 ** =10: nutations in lon & obliq (if on file) **
543 ** =11: lunar librations (if on file) **
544 ** =12: lunar mantle omegas **
545 ** =13: TT-TDB (if on file) **
546 ** **
547 ** Note that I've not actually seen case 12 yet. It probably doesn't work. **
548 ** **
549 ** output: **
550 ** **
551 ** pv[][6] double array that will contain requested interpolated **
552 ** quantities. The body specified by list[i] will have its **
553 ** state in the array starting at pv[i][0] (on any given **
554 ** call, only those words in 'pv' which are affected by the **
555 ** first 10 'list' entries (and by list(11) if librations are **
556 ** on the file) are set. The rest of the 'pv' array **
557 ** is untouched.) The order of components in pv[][] is: **
558 ** pv[][0]=x,....pv[][5]=dz. **
559 ** **
560 ** All output vectors are referenced to the earth mean **
561 ** equator and equinox of epoch. The moon state is always **
562 ** geocentric; the other nine states are either heliocentric **
563 ** or solar-system barycentric, depending on the setting of **
564 ** global variables (see below). **
565 ** **
566 ** Lunar librations, if on file, are put into pv[10][k] if **
567 ** list[11] is 1 or 2. **
568 ** **
569 ** nut dp 4-word array that will contain nutations and rates, **
570 ** depending on the setting of list[10]. the order of **
571 ** quantities in nut is: **
572 ** **
573 ** d psi (nutation in longitude) **
574 ** d epsilon (nutation in obliquity) **
575 ** d psi dot **
576 ** d epsilon dot **
577 ** **
578 *****************************************************************************/
jpl_state(void * ephem,const double et,const int list[14],double pv[][6],double nut[4],const int bary)579 int DLL_FUNC jpl_state(void *ephem, const double et, const int list[14],
580 double pv[][6], double nut[4], const int bary)
581 {
582 struct jpl_eph_data *eph = static_cast<struct jpl_eph_data *>(ephem);
583 unsigned i, j, n_intervals;
584 uint32_t nr;
585 double *buf = eph->cache;
586 double t[2];
587 const double block_loc = (et - eph->ephem_start) / eph->ephem_step;
588 bool recompute_pvsun;
589 const double aufac = 1.0 / eph->au;
590
591 /* error return for epoch out of range */
592 if(et < eph->ephem_start || et > eph->ephem_end)
593 return(JPL_EPH_OUTSIDE_RANGE);
594
595 /* calculate record # and relative time in interval */
596
597 nr = static_cast<uint32_t>(block_loc);
598 t[0] = block_loc - static_cast<double>(nr);
599 if(t[0]==0.0 && nr)
600 {
601 t[0] = 1.;
602 nr--;
603 }
604
605 /* read correct record if not in core (static vector buf[]) */
606 if(nr != eph->curr_cache_loc)
607 {
608 eph->curr_cache_loc = nr;
609 /* Read two blocks ahead to account for header: */
610 if(FSeek(eph->ifile, (nr + 2) * eph->recsize, SEEK_SET)) // lgtm [cpp/integer-multiplication-cast-to-long]
611 {
612 // GZ: Make sure we will try again on next call...
613 eph->curr_cache_loc=0;
614 return(JPL_EPH_FSEEK_ERROR);
615 }
616 if(fread(buf, sizeof(double), static_cast<size_t>(eph->ncoeff), eph->ifile)
617 != static_cast<size_t>(eph->ncoeff))
618 return(JPL_EPH_READ_ERROR);
619
620 if(eph->swap_bytes)
621 swap_64_bit_val(buf, eph->ncoeff);
622 }
623 t[1] = eph->ephem_step;
624
625 if(eph->pvsun_t != et) /* If several calls are made for the same et, */
626 { /* don't recompute pvsun each time... only on */
627 recompute_pvsun = true; /* the first run through. */
628 eph->pvsun_t = et;
629 }
630 else
631 recompute_pvsun = false;
632
633 /* Here, i loops through the "traditional" 14 listed items -- 10
634 solar system objects, nutations, librations, lunar mantle angles,
635 and TT-TDT -- plus a fifteenth: the solar system barycenter. That
636 last is quite different: it's computed 'as needed', rather than
637 from list[]; the output goes to pvsun rather than the pv array;
638 and three quantities (position, velocity, acceleration) are
639 computed (nobody else gets accelerations at present.) */
640 for(n_intervals = 1; n_intervals <= 8; n_intervals *= 2)
641 for(i = 0; i < 15; i++)
642 {
643 unsigned quantities;
644 uint32_t *iptr = &eph->ipt[i + 1][0];
645
646 if(i == 14)
647 {
648 quantities = (recompute_pvsun ? 3 : 0);
649 iptr = &eph->ipt[10][0];
650 }
651 else
652 {
653 quantities = static_cast<unsigned int>(list[i]);
654 iptr = &eph->ipt[i < 10 ? i : i + 1][0];
655 }
656 if(n_intervals == iptr[2] && quantities)
657 {
658 double *dest = ((i == 10) ? eph->pvsun : pv[i]);
659
660 if(i < 10)
661 dest = pv[i];
662 else if(i == 14)
663 dest = eph->pvsun;
664 else
665 dest = nut;
666
667 interp(&eph->iinfo, &buf[iptr[0]-1], t, static_cast<unsigned int>(iptr[1]),
668 dimension(i + 1),
669 n_intervals, static_cast<int>(quantities), dest);
670
671 if(i < 10 || i == 14) /* convert km to AU */
672 for(j = 0; j < quantities * 3; j++)
673 dest[j] *= aufac;
674 }
675 }
676 if(!bary) /* gotta correct everybody for */
677 for(i = 0; i < 9; i++) /* the solar system barycenter */
678 for(j = 0; j < static_cast<unsigned>(list[i]) * 3; j++)
679 pv[i][j] -= eph->pvsun[j];
680 return(0);
681 }
682
683 static int init_err_code = JPL_INIT_NOT_CALLED;
684
jpl_init_error_code(void)685 int DLL_FUNC jpl_init_error_code(void)
686 {
687 return(init_err_code);
688 }
689
jpl_init_error_message(void)690 const char * jpl_init_error_message(void)
691 {
692 switch(init_err_code)
693 {
694 case 0:
695 return "JPL_INIT_NO_ERROR";
696 case -1:
697 return "JPL_INIT_FILE_NOT_FOUND";
698 case -2:
699 return "JPL_INIT_FSEEK_FAILED";
700 case -3:
701 return "JPL_INIT_FREAD_FAILED";
702 case -4:
703 return "JPL_INIT_FREAD2_FAILED";
704 case -5:
705 return "JPL_INIT_FILE_CORRUPT";
706 case -6:
707 return "JPL_INIT_MEMORY_FAILURE";
708 case -7:
709 return "JPL_INIT_FREAD3_FAILED";
710 case -8:
711 return "JPL_INIT_FREAD4_FAILED";
712 case -9:
713 return "JPL_INIT_NOT_CALLED";
714 case -10:
715 return "JPL_INIT_FREAD5_FAILED";
716 default:
717 return "ERROR_NOT_RECOGNIZED";
718 }
719 }
720
721 /* DE-430 has 572 constants. That's more than the 400 constants */
722 /* originally expected. The remaining 172 are stored after the */
723 /* other header data : */
724
725 #define START_400TH_CONSTANT_NAME (84 * 3 + 400 * 6 + 5 * sizeof(double) \
726 + 41 * sizeof(int32_t))
727
728 /* ...which comes out to 2856. See comments in 'jpl_int.h'. */
729
730 /****************************************************************************
731 ** jpl_init_ephemeris(ephemeris_filename, nam, val, n_constants) **
732 *****************************************************************************
733 ** **
734 ** this function does the initial prep work for use of binary JPL **
735 ** ephemerides. **
736 ** const char *ephemeris_filename = full path/filename of the binary **
737 ** ephemeris (on the Willmann-Bell CDs, this is UNIX.200, 405, **
738 ** or 406) **
739 ** char nam[][6] = array of constant names (max 6 characters each) **
740 ** You can pass nam=NULL if you don't care about the names **
741 ** double *val = array of values of constants **
742 ** You can pass val=NULL if you don't care about the constants **
743 ** Return value is a pointer to the jpl_eph_data structure **
744 ** NULL is returned if the file isn't opened or memory isn't alloced **
745 ** Errors can be determined with the above jpl_init_error_code() **
746 ****************************************************************************/
747
jpl_init_ephemeris(const char * ephemeris_filename,char nam[][6],double * val)748 void * DLL_FUNC jpl_init_ephemeris(const char *ephemeris_filename,
749 char nam[][6], double *val)
750 {
751 unsigned i, j;
752 unsigned long de_version;
753 char title[84];
754 FILE *ifile = fopen(ephemeris_filename, "rb");
755
756 struct jpl_eph_data *rval;
757 struct jpl_eph_data temp_data;
758
759 init_err_code = 0;
760 temp_data.ifile = ifile;
761 if(!ifile)
762 init_err_code = JPL_INIT_FILE_NOT_FOUND;
763 else if(fread(title, 84, 1, ifile) != 1)
764 init_err_code = JPL_INIT_FREAD_FAILED;
765 else if(FSeek(ifile, 2652L, SEEK_SET))
766 init_err_code = JPL_INIT_FSEEK_FAILED;
767 else if(fread(&temp_data, JPL_HEADER_SIZE, 1, ifile) != 1)
768 init_err_code = JPL_INIT_FREAD2_FAILED;
769
770 if(init_err_code)
771 {
772 if(ifile)
773 fclose(ifile);
774 return(Q_NULLPTR);
775 }
776
777 de_version = strtoul(title+26, Q_NULLPTR, 10);
778
779 /* A small piece of trickery: in the binary file, data is stored */
780 /* for ipt[0...11], then the ephemeris version, then the */
781 /* remaining ipt[12] data. A little switching is required to get */
782 /* the correct order. */
783 temp_data.ipt[12][0] = temp_data.ipt[12][1];
784 temp_data.ipt[12][1] = temp_data.ipt[12][2];
785 temp_data.ipt[12][2] = temp_data.ipt[13][0];
786 temp_data.ephemeris_version = de_version;
787
788 //qDebug() << "DE_Version: " << de_version;
789
790
791 temp_data.swap_bytes = (temp_data.ncon > 65536L);
792 if(temp_data.swap_bytes) /* byte order is wrong for current platform */
793 {
794 qDebug() << "Byte order is wrong for current platform";
795
796 swap_64_bit_val(&temp_data.ephem_start, 1);
797 swap_64_bit_val(&temp_data.ephem_end, 1);
798 swap_64_bit_val(&temp_data.ephem_step, 1);
799 swap_32_bit_val(&temp_data.ncon);
800 swap_64_bit_val(&temp_data.au, 1);
801 swap_64_bit_val(&temp_data.emrat, 1);
802 }
803
804 /* It's a little tricky to tell if an ephemeris really has */
805 /* TT-TDB data (with offsets in ipt[13][] and ipt[14][]). */
806 /* Essentially, we read the data and sanity-check it, and */
807 /* zero it if it "doesn't add up" correctly. */
808 /* Also: certain ephems I've generated with ncon capped */
809 /* at 400 have no TT-TDB data. So if ncon == 400, don't */
810 /* try to read such data; you may get garbage. */
811 if(de_version >= 430 && temp_data.ncon != 400)
812 {
813 /* If there are 400 or fewer constants, data for ipt[13][0...2] */
814 /* immediately follows that for ipt[12][0..2]; i.e., we don't */
815 /* need to fseek(). Otherwise, we gotta skip 6*(n_constants-400) */
816 /* bytes. */
817 if(temp_data.ncon > 400)
818 {
819 if ( FSeek(ifile, static_cast<long long>(static_cast<size_t>(temp_data.ncon - 400) * 6), SEEK_CUR) != 0)
820 {
821 qWarning() << "jpl_init_ephemeris(): Cannot seek in file. Result will be undefined.";
822 }
823 }
824 if(fread(&temp_data.ipt[13][0], sizeof(int32_t), 6, ifile) != 6)
825 init_err_code = JPL_INIT_FREAD5_FAILED;
826 }
827 else /* mark header data as invalid */
828 temp_data.ipt[13][0] = static_cast<uint32_t>(-1);
829
830 if(temp_data.swap_bytes) /* byte order is wrong for current platform */
831 {
832 for(j = 0; j < 3; j++)
833 {
834 for(i = 0; i < 15; i++)
835 {
836 swap_32_bit_val(&temp_data.ipt[i][j]);
837 }
838 }
839 }
840
841 if(temp_data.ipt[13][0] != /* if these don't add up correctly, */
842 temp_data.ipt[12][0] + temp_data.ipt[12][1] * temp_data.ipt[12][2] * 3
843 || temp_data.ipt[14][0] != /* zero them out (they've garbage data) */
844 temp_data.ipt[13][0] + temp_data.ipt[13][1] * temp_data.ipt[13][2] * 3)
845 { /* not valid pointers to TT-TDB data */
846 memset(&temp_data.ipt[13][0], 0, 6 * sizeof(int32_t));
847 }
848
849 /* A sanity check: if the earth-moon ratio is outside reasonable */
850 /* bounds, we must be looking at a wrong or corrupted file. */
851 /* In DE-102, emrat = 81.3007; in DE-405/406, emrat = 81.30056. */
852 /* Those are the low and high ranges. We'll allow some slop in */
853 /* case the earth/moon mass ratio changes: */
854 if(temp_data.emrat > 81.3008 || temp_data.emrat < 81.30055)
855 {
856 init_err_code = JPL_INIT_FILE_CORRUPT;
857 qWarning() << "temp_data: " << temp_data.emrat << "(should have been =~81). JPL_INIT_FILE_CORRUPT!";
858 }
859
860 if(init_err_code)
861 {
862 fclose(ifile);
863 return(Q_NULLPTR);
864 }
865
866 /* Once upon a time, the kernel size was determined from the */
867 /* DE version. This was not a terrible idea, except that it */
868 /* meant that when the code faced a new version, it broke. */
869 /* Now we use some logic to compute the kernel size. */
870 temp_data.kernel_size = 4;
871 for(i = 0; i < 15; i++)
872 temp_data.kernel_size +=
873 2 * temp_data.ipt[i][1] * temp_data.ipt[i][2] * dimension(i);
874 // for(i = 0; i < 13; i++)
875 // temp_data.kernel_size +=
876 // temp_data.ipt[i][1] * temp_data.ipt[i][2] * ((i == 11) ? 4 : 6);
877 // /* ...and then add in space required for the TT-TDB data : */
878 // temp_data.kernel_size += temp_data.ipt[14][1] * temp_data.ipt[14][2] * 2;
879 temp_data.recsize = temp_data.kernel_size * 4L;
880 temp_data.ncoeff = temp_data.kernel_size / 2L;
881
882 /* Rather than do two separate allocations, everything */
883 /* we need is allocated in _one_ chunk, then parceled out. */
884 /* This looks a little weird, but it simplifies error */
885 /* handling and cleanup. */
886 rval = static_cast<struct jpl_eph_data *>(calloc(sizeof(struct jpl_eph_data)
887 + temp_data.recsize, 1));
888 if(!rval)
889 {
890 init_err_code = JPL_INIT_MEMORY_FAILURE;
891 fclose(ifile);
892 return(Q_NULLPTR);
893 }
894 memcpy(rval, &temp_data, sizeof(struct jpl_eph_data));
895 rval->iinfo.posn_coeff[0] = 1.0;
896 /* Seed a bogus value here. The first and subsequent calls to */
897 /* 'interp' will correct it to a value between -1 and +1. */
898 rval->iinfo.posn_coeff[1] = -2.0;
899 rval->iinfo.vel_coeff[0] = 0.0;
900 rval->iinfo.vel_coeff[1] = 1.0;
901 rval->curr_cache_loc = static_cast<uint32_t>(-1);
902
903 /* The 'cache' data is right after the 'jpl_eph_data' struct: */
904 rval->cache = reinterpret_cast<double *>(rval + 1);
905 /* If there are more than 400 constants, the names of */
906 /* the extra constants are stored in what would normally */
907 /* be zero-padding after the header record. However, */
908 /* older ephemeris-reading software won't know about that. */
909 /* So we store ncon=400, then actually check the names to */
910 /* see how many constants there _really_ are. Older readers */
911 /* will just see 400 names and won't know about the others. */
912 /* But on the upside, they won't crash. */
913
914 if(rval->ncon == 400)
915 {
916 char buff[7];
917
918 buff[6] = '\0';
919 FSeek(ifile, START_400TH_CONSTANT_NAME, SEEK_SET);
920 while(fread(buff, 6, 1, ifile) && strlen(buff) == 6)
921 {
922 rval->ncon++;
923 }
924 }
925
926 if(val)
927 {
928 FSeek(ifile, rval->recsize, SEEK_SET);
929 if(fread(val, sizeof(double), static_cast<size_t>(rval->ncon), ifile)
930 != static_cast<size_t>(rval->ncon))
931 init_err_code = JPL_INIT_FREAD3_FAILED;
932 else if(rval->swap_bytes) /* gotta swap the constants, too */
933 swap_64_bit_val(val, rval->ncon);
934 }
935
936 if(!init_err_code && nam)
937 {
938 FSeek(ifile, 84L * 3L, SEEK_SET); /* just after the 3 'title' lines */
939 for(i = 0; i < rval->ncon && !init_err_code; i++)
940 {
941 if(i == 400)
942 FSeek(ifile, START_400TH_CONSTANT_NAME, SEEK_SET);
943 if(fread(nam[i], 6, 1, ifile) != 1)
944 init_err_code = JPL_INIT_FREAD4_FAILED;
945 }
946 }
947 return(rval);
948 }
949
950 /****************************************************************************
951 ** jpl_close_ephemeris(ephem) **
952 *****************************************************************************
953 ** **
954 ** this function closes files and frees up memory allocated by the **
955 ** jpl_init_ephemeris() function. **
956 ****************************************************************************/
jpl_close_ephemeris(void * ephem)957 void DLL_FUNC jpl_close_ephemeris(void *ephem)
958 {
959 struct jpl_eph_data *eph = static_cast<struct jpl_eph_data *>(ephem);
960
961 fclose(eph->ifile);
962 free(ephem);
963 }
964
965 /* Added 2011 Jan 18: random access to any desired JPL constant */
966
967
jpl_get_constant(const int idx,void * ephem,char * constant_name)968 double DLL_FUNC jpl_get_constant(const int idx, void *ephem, char *constant_name)
969 {
970 struct jpl_eph_data *eph = static_cast<struct jpl_eph_data *>(ephem);
971 double rval = 0.;
972
973 *constant_name = '\0';
974 if(idx >= 0 && idx < static_cast<int>(eph->ncon))
975 {
976 // GZ extended from const long to const long long
977 const long long seek_loc = (idx < 400 ? 84L * 3L + static_cast<long>(idx) * 6 :
978 static_cast<long long>(START_400TH_CONSTANT_NAME) + (idx - 400) * 6);
979
980 if (FSeek(eph->ifile, seek_loc, SEEK_SET) != 0)
981 {
982 qWarning() << "jpl_get_constant(): Cannot seek in file. Result will be undefined.";
983 }
984 if(fread(constant_name, 1, 6, eph->ifile) == 6 )
985 {
986 constant_name[6] = '\0';
987 if (FSeek(eph->ifile, static_cast<long long>(eph->recsize) + static_cast<long long>(idx) * static_cast<long long>(sizeof(double)), SEEK_SET) != 0)
988 {
989 qWarning() << "jpl_get_constant(): Cannot seek in file (call2). Result will be undefined.";
990 }
991
992 // GZ added tests to make travis test suite (Coverity Scan) happier
993 if( fread(&rval, 1, sizeof(double), eph->ifile) == sizeof(double) )
994 {
995 if(eph->swap_bytes) /* gotta swap the constants, too */
996 swap_64_bit_val(&rval, 1);
997 }
998 else
999 qWarning() << "jpl_get_constant(): fread() failed to read a double. Result will be undefined.";
1000 }
1001 else
1002 qWarning() << "jpl_get_constant(): fread() failed to read a constant name. Result will be 0, likely wrong.";
1003 }
1004 return(rval);
1005 }
1006 /*************************** THE END ***************************************/
1007