1 /*****************************************************************************
2 * ***** jpl planetary and lunar ephemerides ***** C ver.1.2 *
3 ******************************************************************************
4 * *
5 * This program was written in standard fortran-77 and it was manually *
6 * translated to C language by Piotr A. Dybczynski (dybol@phys.amu.edu.pl), *
7 * subsequently revised heavily by Bill J Gray (pluto@gwi.net), just short *
8 * of a total re-write. *
9 * *
10 ******************************************************************************
11 * Last modified: July 23, 1997 by PAD *
12 ******************************************************************************
13 21 Apr 2010: Revised by Bill J. Gray. The code now determines the kernel
14 size, then allocates memory accordingly. This should 'future-proof' us in
15 case JPL (or someone else) creates kernels that are larger than the previously
16 arbitrary MAX_KERNEL_SIZE parameter. 'swap_long' and 'swap_double' have
17 been replaced with 'swap_32_bit_val' and 'swap_64_bit_val'. It also works
18 on 64-bit compiles now.
19
20 16 Mar 2001: Revised by Bill J. Gray. You can now use binary
21 ephemerides with either byte order ('big-endian' or 'small-endian');
22 the code checks to see if the data is in the "wrong" order for the
23 current platform, and swaps bytes on-the-fly if needed. (Yes, this
24 can result in a slowdown... sometimes as much as 1%. The function is
25 so mathematically intensive that the byte-swapping is the least of our
26 troubles.) You can also use DE-200, 403, 404, 405, or 406 without
27 recompiling (the constan( ) function now determines which ephemeris is
28 in use and its byte order).
29
30 Also, I did some minor optimization of the interp( ) (Chebyshev
31 interpolation) function, resulting in a bit of a speedup.
32
33 The code has been modified to be a separately linkable component, with
34 details of the implementation encapsulated.
35 *****************************************************************************/
36
37 #include <stdio.h>
38 #include <math.h>
39 #include <string.h>
40 #include <stdlib.h>
41 #ifdef _MSC_VER /* Microsoft Visual C/C++ lacks a 'stdint.h'; */
42 #include "stdintvc.h" /* 'stdintvc.h' is a replacement version */
43 #else
44 #include <stdint.h>
45 #endif
46
47 /**** include variable and type definitions, specific for this C version */
48
49 #include "jpleph.h"
50 #include "jpl_int.h"
51
52 #define TRUE 1
53 #define FALSE 0
54
jpl_get_double(const void * ephem,const int value)55 double DLL_FUNC jpl_get_double( const void *ephem, const int value)
56 {
57 return( *(double *)( (char *)ephem + value));
58 }
59
jpl_get_long(const void * ephem,const int value)60 long DLL_FUNC jpl_get_long( const void *ephem, const int value)
61 {
62 return( *(int32_t *)( (char *)ephem + value));
63 }
64
65
66 /*****************************************************************************
67 ** jpl_pleph( ephem,et,ntar,ncent,rrd,calc_velocity) **
68 ******************************************************************************
69 ** **
70 ** This subroutine reads the jpl planetary ephemeris **
71 ** and gives the position and velocity of the point 'ntarg' **
72 ** with respect to 'ncent'. **
73 ** **
74 ** Calling sequence parameters: **
75 ** **
76 ** et = (double) julian ephemeris date at which interpolation **
77 ** is wanted. **
78 ** **
79 ** ntarg = integer number of 'target' point. **
80 ** **
81 ** ncent = integer number of center point. **
82 ** **
83 ** The numbering convention for 'ntarg' and 'ncent' is: **
84 ** **
85 ** 1 = mercury 8 = neptune **
86 ** 2 = venus 9 = pluto **
87 ** 3 = earth 10 = moon **
88 ** 4 = mars 11 = sun **
89 ** 5 = jupiter 12 = solar-system barycenter **
90 ** 6 = saturn 13 = earth-moon barycenter **
91 ** 7 = uranus 14 = nutations (longitude and obliq) **
92 ** 15 = librations, if on eph. file **
93 ** **
94 ** (If nutations are wanted, set ntarg = 14. **
95 ** For librations, set ntarg = 15. set ncent= 0) **
96 ** **
97 ** rrd = output 6-element, double array of position and velocity **
98 ** of point 'ntarg' relative to 'ncent'. The units are au and **
99 ** au/day. For librations the units are radians and radians **
100 ** per day. In the case of nutations the first four words of **
101 ** rrd will be set to nutations and rates, having units of **
102 ** radians and radians/day. **
103 ** **
104 ** The option is available to have the units in km and km/sec. **
105 ** for this, set km=TRUE at the begining of the program. **
106 ** **
107 ** calc_velocity = integer flag; if nonzero, velocities will be **
108 ** computed, otherwise not. **
109 ** **
110 *****************************************************************************/
jpl_pleph(void * ephem,const double et,const int ntarg,const int ncent,double rrd[],const int calc_velocity)111 int DLL_FUNC jpl_pleph( void *ephem, const double et, const int ntarg,
112 const int ncent, double rrd[], const int calc_velocity)
113 {
114 struct jpl_eph_data *eph = (struct jpl_eph_data *)ephem;
115 double pv[13][6];/* pv is the position/velocity array
116 NUMBERED FROM ZERO: 0=Mercury,1=Venus,...
117 8=Pluto,9=Moon,10=Sun,11=SSBary,12=EMBary
118 First 10 elements (0-9) are affected by
119 jpl_state(), all are adjusted here. */
120
121
122 int rval = 0, list_val = (calc_velocity ? 2 : 1);
123 int i, k, list[12]; /* list is a vector denoting, for which "body"
124 ephemeris values should be calculated by
125 jpl_state(): 0=Mercury,1=Venus,2=EMBary,...,
126 8=Pluto, 9=geocentric Moon, 10=nutations in
127 long. & obliq. 11= lunar librations */
128
129 for( i = 0; i < 6; ++i) rrd[i] = 0.0;
130
131 if( ntarg == ncent) return( 0);
132
133 for( i = 0; i < 12; i++) list[i] = 0;
134
135 /* check for nutation call */
136
137 if( ntarg == 14)
138 {
139 if( eph->ipt[11][1] > 0) /* there is nutation on ephemeris */
140 {
141 list[10] = list_val;
142 if( jpl_state( ephem, et, list, pv, rrd, 0))
143 rval = -1;
144 }
145 else /* no nutations on the ephemeris file */
146 rval = -2;
147 return( rval);
148 }
149
150 /* check for librations */
151
152 if( ntarg == 15)
153 {
154 if( eph->ipt[12][1] > 0) /* there are librations on ephemeris file */
155 {
156 list[11] = list_val;
157 if( jpl_state( eph, et, list, pv, rrd, 0))
158 rval = -3;
159 for( i = 0; i < 6; ++i)
160 rrd[i] = pv[10][i]; /* librations */
161 }
162 else /* no librations on the ephemeris file */
163 rval = -4;
164 return( rval);
165 }
166
167 /* force barycentric output by 'state' */
168
169 /* set up proper entries in 'list' array for state call */
170
171 for( i = 0; i < 2; i++) /* list[] IS NUMBERED FROM ZERO ! */
172 {
173 k = ntarg-1;
174 if( i == 1) k=ncent-1; /* same for ntarg & ncent */
175 if( k <= 9) list[k] = list_val; /* Major planets */
176 if( k == 9) list[2] = list_val; /* for moon, earth state is needed */
177 if( k == 2) list[9] = list_val; /* for earth, moon state is needed */
178 if( k == 12) list[2] = list_val; /* EMBary state additionaly */
179 }
180
181 /* make call to state */
182
183 if( jpl_state( eph, et, list, pv, rrd, 1))
184 rval = -5;
185 /* Solar System barycentric Sun state goes to pv[10][] */
186 if( ntarg == 11 || ncent == 11)
187 for( i = 0; i < 6; i++)
188 pv[10][i] = eph->pvsun[i];
189
190 /* Solar System Barycenter coordinates & velocities equal to zero */
191 if( ntarg == 12 || ncent == 12)
192 for( i = 0; i < 6; i++)
193 pv[11][i] = 0.0;
194
195 /* Solar System barycentric EMBary state: */
196 if( ntarg == 13 || ncent == 13)
197 for( i = 0; i < 6; i++)
198 pv[12][i] = pv[2][i];
199
200 /* if moon from earth or earth from moon ..... */
201 if( (ntarg*ncent) == 30 && (ntarg+ncent) == 13)
202 for( i = 0; i < 6; ++i) pv[2][i]=0.0;
203 else
204 {
205 if( list[2]) /* calculate earth state from EMBary */
206 for( i = 0; i < list[2] * 3; ++i)
207 pv[2][i] -= pv[9][i]/(1.0+eph->emrat);
208
209 if(list[9]) /* calculate Solar System barycentric moon state */
210 for( i = 0; i < list[9] * 3; ++i)
211 pv[9][i] += pv[2][i];
212 }
213
214 for( i = 0; i < list_val * 3; ++i)
215 rrd[i] = pv[ntarg-1][i] - pv[ncent-1][i];
216
217 return( rval);
218 }
219
220 /*****************************************************************************
221 ** interp(buf,t,ncf,ncm,na,ifl,pv) **
222 ******************************************************************************
223 ** **
224 ** this subroutine differentiates and interpolates a **
225 ** set of chebyshev coefficients to give position and velocity **
226 ** **
227 ** calling sequence parameters: **
228 ** **
229 ** input: **
230 ** **
231 ** iinfo stores certain chunks of interpolation info, in hopes **
232 ** that if you call again with similar parameters, the **
233 ** function won't have to re-compute all coefficients/data. **
234 ** **
235 ** coef 1st location of array of d.p. chebyshev coefficients **
236 ** of position **
237 ** **
238 ** t t[0] is double fractional time in interval covered by **
239 ** coefficients at which interpolation is wanted **
240 ** (0 <= t[0] <= 1). t[1] is dp length of whole **
241 ** interval in input time units. **
242 ** **
243 ** ncf # of coefficients per component **
244 ** **
245 ** ncm # of components per set of coefficients **
246 ** **
247 ** na # of sets of coefficients in full array **
248 ** (i.e., # of sub-intervals in full interval) **
249 ** **
250 ** ifl integer flag: =1 for positions only **
251 ** =2 for pos and vel **
252 ** **
253 ** **
254 ** output: **
255 ** **
256 ** posvel interpolated quantities requested. dimension **
257 ** expected is posvel[ncm*ifl], double precision. **
258 ** **
259 *****************************************************************************/
interp(struct interpolation_info * iinfo,const double coef[],const double t[2],const int ncf,const int ncm,const int na,const int ifl,double posvel[])260 static void interp( struct interpolation_info *iinfo,
261 const double coef[], const double t[2], const int ncf,
262 const int ncm, const int na, const int ifl, double posvel[])
263 {
264 double dna, dt1, temp, tc, vfac, temp1;
265 int l, i, j;
266
267 /* entry point. get correct sub-interval number for this set
268 of coefficients and then get normalized chebyshev time
269 within that subinterval. */
270
271 dna = (double)na;
272 modf( t[0], &dt1);
273 temp = dna * t[0];
274 l = (int)(temp - dt1);
275
276 /* tc is the normalized chebyshev time (-1 <= tc <= 1) */
277
278 tc = 2.0 * (modf( temp, &temp1) + dt1) - 1.0;
279
280 /* check to see whether chebyshev time has changed,
281 and compute new polynomial values if it has.
282 (the element iinfo->pc[1] is the value of t1[tc] and hence
283 contains the value of tc on the previous call.) */
284
285 if(tc != iinfo->pc[1])
286 {
287 iinfo->np = 2;
288 iinfo->nv = 3;
289 iinfo->pc[1] = tc;
290 iinfo->twot = tc+tc;
291 }
292
293 /* be sure that at least 'ncf' polynomials have been evaluated
294 and are stored in the array 'iinfo->pc'. */
295
296 if(iinfo->np < ncf)
297 {
298 double *pc_ptr = iinfo->pc + iinfo->np;
299
300 for(i=ncf - iinfo->np; i; i--, pc_ptr++)
301 *pc_ptr = iinfo->twot * pc_ptr[-1] - pc_ptr[-2];
302 iinfo->np=ncf;
303 }
304
305 /* interpolate to get position for each component */
306
307 for( i = 0; i < ncm; ++i) /* ncm is a number of coordinates */
308 {
309 const double *coeff_ptr = coef + ncf * (i + l * ncm + 1);
310 const double *pc_ptr = iinfo->pc + ncf;
311
312 posvel[i]=0.0;
313 for( j = ncf; j; j--)
314 posvel[i] += (*--pc_ptr) * (*--coeff_ptr);
315 }
316
317 if(ifl <= 1) return;
318
319 /* if velocity interpolation is wanted, be sure enough
320 derivative polynomials have been generated and stored. */
321
322 vfac=(dna+dna)/t[1];
323 iinfo->vc[2] = iinfo->twot + iinfo->twot;
324 if( iinfo->nv < ncf)
325 {
326 double *vc_ptr = iinfo->vc + iinfo->nv;
327 const double *pc_ptr = iinfo->pc + iinfo->nv - 1;
328
329 for( i = ncf - iinfo->nv; i; i--, vc_ptr++, pc_ptr++)
330 *vc_ptr = iinfo->twot * vc_ptr[-1] + *pc_ptr + *pc_ptr - vc_ptr[-2];
331 iinfo->nv = ncf;
332 }
333
334 /* interpolate to get velocity for each component */
335
336 for( i = 0; i < ncm; ++i)
337 {
338 double tval = 0.;
339 const double *coeff_ptr = coef + ncf * (i + l * ncm + 1);
340 const double *vc_ptr = iinfo->vc + ncf;
341
342 for( j = ncf; j; j--)
343 tval += (*--vc_ptr) * (*--coeff_ptr);
344 posvel[ i + ncm] = tval * vfac;
345 }
346 return;
347 }
348
349 /* swap_32_bit_val() and swap_64_bit_val() are used when reading a binary
350 ephemeris that was created on a machine with 'opposite' byte order to
351 the currently-used machine (signalled by the 'swap_bytes' flag in the
352 jpl_eph_data structure). In such cases, every double and integer
353 value read from the ephemeris must be byte-swapped by these two functions. */
354
355 #define SWAP_MACRO( A, B, TEMP) { TEMP = A; A = B; B = TEMP; }
356
swap_32_bit_val(void * ptr)357 static void swap_32_bit_val( void *ptr)
358 {
359 char *tptr = (char *)ptr, tchar;
360
361 SWAP_MACRO( tptr[0], tptr[3], tchar);
362 SWAP_MACRO( tptr[1], tptr[2], tchar);
363 }
364
swap_64_bit_val(void * ptr,long count)365 static void swap_64_bit_val( void *ptr, long count)
366 {
367 char *tptr = (char *)ptr, tchar;
368
369 while( count--)
370 {
371 SWAP_MACRO( tptr[0], tptr[7], tchar);
372 SWAP_MACRO( tptr[1], tptr[6], tchar);
373 SWAP_MACRO( tptr[2], tptr[5], tchar);
374 SWAP_MACRO( tptr[3], tptr[4], tchar);
375 tptr += 8;
376 }
377 }
378
379 /*****************************************************************************
380 ** jpl_state(ephem,et2,list,pv,nut,bary) **
381 ******************************************************************************
382 ** This subroutine reads and interpolates the jpl planetary ephemeris file **
383 ** **
384 ** Calling sequence parameters: **
385 ** **
386 ** Input: **
387 ** **
388 ** et2[] double, 2-element JED epoch at which interpolation **
389 ** is wanted. Any combination of et2[0]+et2[1] which falls **
390 ** within the time span on the file is a permissible epoch. **
391 ** **
392 ** a. for ease in programming, the user may put the **
393 ** entire epoch in et2[0] and set et2[1]=0.0 **
394 ** **
395 ** b. for maximum interpolation accuracy, set et2[0] = **
396 ** the most recent midnight at or before interpolation **
397 ** epoch and set et2[1] = fractional part of a day **
398 ** elapsed between et2[0] and epoch. **
399 ** **
400 ** c. as an alternative, it may prove convenient to set **
401 ** et2[0] = some fixed epoch, such as start of integration,**
402 ** and et2[1] = elapsed interval between then and epoch. **
403 ** **
404 ** list 12-element integer array specifying what interpolation **
405 ** is wanted for each of the "bodies" on the file. **
406 ** **
407 ** list[i]=0, no interpolation for body i **
408 ** =1, position only **
409 ** =2, position and velocity **
410 ** **
411 ** the designation of the astronomical bodies by i is: **
412 ** **
413 ** i = 0: mercury **
414 ** = 1: venus **
415 ** = 2: earth-moon barycenter **
416 ** = 3: mars **
417 ** = 4: jupiter **
418 ** = 5: saturn **
419 ** = 6: uranus **
420 ** = 7: neptune **
421 ** = 8: pluto **
422 ** = 9: geocentric moon **
423 ** =10: nutations in longitude and obliquity **
424 ** =11: lunar librations (if on file) **
425 ** **
426 ** output: **
427 ** **
428 ** pv[][6] double array that will contain requested interpolated **
429 ** quantities. The body specified by list[i] will have its **
430 ** state in the array starting at pv[i][0] (on any given **
431 ** call, only those words in 'pv' which are affected by the **
432 ** first 10 'list' entries (and by list(11) if librations are **
433 ** on the file) are set. The rest of the 'pv' array **
434 ** is untouched.) The order of components in pv[][] is: **
435 ** pv[][0]=x,....pv[][5]=dz. **
436 ** **
437 ** All output vectors are referenced to the earth mean **
438 ** equator and equinox of epoch. The moon state is always **
439 ** geocentric; the other nine states are either heliocentric **
440 ** or solar-system barycentric, depending on the setting of **
441 ** global variables (see below). **
442 ** **
443 ** Lunar librations, if on file, are put into pv[10][k] if **
444 ** list[11] is 1 or 2. **
445 ** **
446 ** nut dp 4-word array that will contain nutations and rates, **
447 ** depending on the setting of list[10]. the order of **
448 ** quantities in nut is: **
449 ** **
450 ** d psi (nutation in longitude) **
451 ** d epsilon (nutation in obliquity) **
452 ** d psi dot **
453 ** d epsilon dot **
454 ** **
455 *****************************************************************************/
jpl_state(void * ephem,const double et,const int list[12],double pv[][6],double nut[4],const int bary)456 int DLL_FUNC jpl_state( void *ephem, const double et, const int list[12],
457 double pv[][6], double nut[4], const int bary)
458 {
459 struct jpl_eph_data *eph = (struct jpl_eph_data *)ephem;
460 int i,j, n_intervals;
461 long int nr;
462 double prev_midnight, time_of_day;
463 double *buf = eph->cache;
464 double s,t[2],aufac;
465 struct interpolation_info *iinfo = (struct interpolation_info *)eph->iinfo;
466
467
468 /* ********** main entry point ********** */
469
470 s=et - 0.5;
471 prev_midnight = floor( s);
472 time_of_day = s - prev_midnight;
473 prev_midnight += 0.5;
474
475 /* here prev_midnight contains last midnight before epoch desired (in JED: *.5)
476 and time_of_day contains the remaining, fractional part of the epoch */
477
478 /* error return for epoch out of range */
479
480 if( et < eph->ephem_start || et > eph->ephem_end)
481 return( -1);
482
483 /* calculate record # and relative time in interval */
484
485 nr = (long)((prev_midnight-eph->ephem_start)/eph->ephem_step)+2;
486 /* add 2 to adjus for the first two records containing header data */
487 if( prev_midnight == eph->ephem_end)
488 nr--;
489 t[0]=( prev_midnight-( (1.0*nr-2.0)*eph->ephem_step+eph->ephem_start) +
490 time_of_day )/eph->ephem_step;
491
492 /* read correct record if not in core (static vector buf[]) */
493
494 if( nr != eph->curr_cache_loc)
495 {
496 eph->curr_cache_loc = nr;
497 fseek( eph->ifile, nr * eph->recsize, SEEK_SET);
498 if( fread( buf, sizeof( double), (size_t)eph->ncoeff, eph->ifile) != (size_t)eph->ncoeff)
499 return( -2);
500 if( eph->swap_bytes)
501 swap_64_bit_val( buf, eph->ncoeff);
502 }
503 t[1] = eph->ephem_step;
504 aufac = 1.0 / eph->au;
505
506 for( n_intervals = 1; n_intervals <= 8; n_intervals *= 2)
507 for( i = 0; i < 11; i++)
508 if( n_intervals == eph->ipt[i][2] && (list[i] || i == 10))
509 {
510 int flag = ((i == 10) ? 2 : list[i]);
511 double *dest = ((i == 10) ? eph->pvsun : pv[i]);
512
513 interp( iinfo, &buf[eph->ipt[i][0]-1], t, (int)eph->ipt[i][1], 3,
514 n_intervals, flag, dest);
515 /* gotta convert units */
516 for( j = 0; j < flag * 3; j++)
517 dest[j] *= aufac;
518 }
519
520 if( !bary) /* gotta correct everybody for */
521 for( i = 0; i < 9; i++) /* the solar system barycenter */
522 for( j = 0; j < list[i] * 3; j++)
523 pv[i][j] -= eph->pvsun[j];
524
525 /* do nutations if requested (and if on file) */
526
527 if(list[10] > 0 && eph->ipt[11][1] > 0)
528 interp( iinfo, &buf[eph->ipt[11][0]-1], t, (int)eph->ipt[11][1], 2,
529 (int)eph->ipt[11][2], list[10], nut);
530
531 /* get librations if requested (and if on file) */
532
533 if(list[11] > 0 && eph->ipt[12][1] > 0)
534 {
535 double pefau[6];
536
537 interp( iinfo, &buf[eph->ipt[12][0]-1], t, (int)eph->ipt[12][1], 3,
538 (int)eph->ipt[12][2], list[11], pefau);
539 for( j = 0; j < 6; ++j) pv[10][j]=pefau[j];
540 }
541 return( 0);
542 }
543
544 /****************************************************************************
545 ** jpl_init_ephemeris( ephemeris_filename, nam, val, n_constants) **
546 *****************************************************************************
547 ** **
548 ** this function does the initial prep work for use of binary JPL **
549 ** ephemerides. **
550 ** const char *ephemeris_filename = full path/filename of the binary **
551 ** ephemeris (on the Willmann-Bell CDs, this is UNIX.200, 405, **
552 ** or 406)
553 ** char nam[][6] = array of constant names (max 6 characters each) **
554 ** You can pass nam=NULL if you don't care about the names **
555 ** double *val = array of values of constants **
556 ** You can pass val=NULL if you don't care about the constants **
557 ** Return value is a pointer to the jpl_eph_data structure **
558 ** NULL is returned if the file isn't opened or memory isn't alloced **
559 ****************************************************************************/
560
561 int jpl_init_err_code = 0;
562
jpl_init_ephemeris(const char * ephemeris_filename,char nam[][6],double * val)563 void * DLL_FUNC jpl_init_ephemeris( const char *ephemeris_filename,
564 char nam[][6], double *val)
565 {
566 int i, j;
567 long de_version;
568 char title[84];
569 FILE *ifile = fopen( ephemeris_filename, "rb");
570 struct jpl_eph_data *rval;
571 struct jpl_eph_data temp_data;
572 struct interpolation_info *iinfo;
573
574 jpl_init_err_code = 0;
575 if( !ifile)
576 {
577 jpl_init_err_code = -1;
578 return( NULL);
579 }
580 temp_data.ifile = ifile;
581 if( fread( title, 84, 1, ifile) != 1)
582 jpl_init_err_code = -3;
583 fseek( ifile, 2652L, SEEK_SET); /* skip title & constant name data */
584 if( fread( &temp_data, JPL_HEADER_SIZE, 1, ifile) != 1)
585 jpl_init_err_code = -4;
586
587 if( jpl_init_err_code)
588 {
589 fclose( ifile);
590 return( NULL);
591 }
592
593 de_version = atoi( title + 26);
594
595 /* A small piece of trickery: in the binary file, data is stored */
596 /* for ipt[0...11], then the ephemeris version, then the */
597 /* remaining ipt[12] data. A little switching is required to get */
598 /* the correct order. */
599 temp_data.ipt[12][0] = temp_data.ipt[12][1];
600 temp_data.ipt[12][1] = temp_data.ipt[12][2];
601 temp_data.ipt[12][2] = temp_data.ephemeris_version;
602 temp_data.ephemeris_version = de_version;
603
604 temp_data.swap_bytes = ( temp_data.ncon < 0 || temp_data.ncon > 65536L);
605 if( temp_data.swap_bytes) /* byte order is wrong for current platform */
606 {
607 swap_64_bit_val( &temp_data.ephem_start, 1);
608 swap_64_bit_val( &temp_data.ephem_end, 1);
609 swap_64_bit_val( &temp_data.ephem_step, 1);
610 swap_32_bit_val( &temp_data.ncon);
611 swap_64_bit_val( &temp_data.au, 1);
612 swap_64_bit_val( &temp_data.emrat, 1);
613 for( j = 0; j < 3; j++)
614 for( i = 0; i < 13; i++)
615 swap_32_bit_val( &temp_data.ipt[i][j]);
616 }
617 /* Once upon a time, the kernel size was determined from the */
618 /* DE version. This was not a terrible idea, except that it */
619 /* meant that when the code faced a new version, it broke. */
620 /* Now we use some logic to compute the kernel size. */
621 temp_data.kernel_size = 4;
622 for( i = 0; i < 13; i++)
623 temp_data.kernel_size +=
624 temp_data.ipt[i][1] * temp_data.ipt[i][2] * ((i == 11) ? 4 : 6);
625 temp_data.recsize = temp_data.kernel_size * 4L;
626 temp_data.ncoeff = temp_data.kernel_size / 2L;
627
628 /* Rather than do three separate allocations, everything */
629 /* we need is allocated in _one_ chunk, then parceled out. */
630 /* This looks a little weird, but it does simplify error */
631 /* handling and cleanup. */
632 rval = (struct jpl_eph_data *)calloc( sizeof( struct jpl_eph_data)
633 + sizeof( struct interpolation_info)
634 + temp_data.recsize, 1);
635 if( !rval)
636 {
637 jpl_init_err_code = -2;
638 fclose( ifile);
639 return( NULL);
640 }
641 memcpy( rval, &temp_data, sizeof( struct jpl_eph_data));
642 /* The 'iinfo' struct is right after the 'jpl_eph_data' struct: */
643 iinfo = (struct interpolation_info *)(rval + 1);
644 rval->iinfo = (void *)iinfo;
645 iinfo->np = 2;
646 iinfo->nv = 3;
647 iinfo->pc[0] = 1.0;
648 iinfo->pc[1] = 0.0;
649 iinfo->vc[1] = 1.0;
650 rval->curr_cache_loc = -1L;
651 /* The 'cache' data is right after the 'iinfo' struct: */
652 rval->cache = (double *)( iinfo + 1);
653
654 if( val)
655 {
656 fseek( ifile, rval->recsize, SEEK_SET);
657 if( fread( val, sizeof( double), (size_t)rval->ncon, ifile)
658 != (size_t)rval->ncon)
659 jpl_init_err_code = -5;
660 if( rval->swap_bytes) /* gotta swap the constants, too */
661 swap_64_bit_val( val, rval->ncon);
662 }
663
664 if( nam)
665 {
666 fseek( ifile, 84L * 3L, SEEK_SET); /* just after the 3 'title' lines */
667 for( i = 0; i < (int)rval->ncon; i++)
668 if( fread( nam[i], 6, 1, ifile) != 1)
669 jpl_init_err_code = -6;
670 }
671 return( rval);
672 }
673
674 /****************************************************************************
675 ** jpl_close_ephemeris( ephem) **
676 *****************************************************************************
677 ** **
678 ** this function closes files and frees up memory allocated by the **
679 ** jpl_init_ephemeris( ) function. **
680 ****************************************************************************/
jpl_close_ephemeris(void * ephem)681 void DLL_FUNC jpl_close_ephemeris( void *ephem)
682 {
683 struct jpl_eph_data *eph = (struct jpl_eph_data *)ephem;
684
685 fclose( eph->ifile);
686 free( ephem);
687 }
688 /*************************** THE END ***************************************/
689
690