00001 /***************************************************************************** 00002 * ***** jpl planetary and lunar ephemerides ***** ver.1.0 * 00003 * Last modified: 4 March 2004 by RSB * 00004 *****************************************************************************/ 00005 /* 00006 * jpleph.c - function pleph() for reading from JPL binary ephemeris table 00007 * and returning relative position-velocity coordinates for selected pair 00008 * of solar system objects 00009 * 00010 * Author: Richard S. Bogart (rbogart@spd.aas.org), Stanford University 00011 * 00012 * This code is based on the JPL routines included in the file testeph.f 00013 * available at ftp://ssd.jpl.nasa.gov/pub/eph/export/fortran/ 00014 * (last modified 23 July 2001) and translated into C by 00015 * Piotr A. Dybczynski (dybol@phys.amu.edu.pl), version 1.2 (23 July 1997) 00016 * available at ftp://ftp.astro.amu.edu.pl/pub/jpleph (see also 00017 * http://www.projectpluto.com/jpl_eph.htm 00018 * 00019 * The following substantive changes have been made to the original code: 00020 * 00021 * An internal byte-reordering function has been added in order to 00022 * support reading of the big-endian JPL-distributed tables on 00023 * little-endian platforms that support the IEEE-754 numeric data 00024 * representations (notably Intel processors). This allows the 00025 * sharing of the same binary tables by machines with multiple 00026 * architectures and obviates the necessity of doing private ASCII 00027 * to binary table conversions. The type of platform is automatically 00028 * detected by comparing redundant binary and ASCII information in the 00029 * JPL tables. 00030 * The type of the JPL table is automatically detected and the read buffer 00031 * sizes and record lengths allocated appropriately. 00032 * An additional function pleph_ext() has been provided with added arguments 00033 * for supporting options that were previously fixed on compilation: 00034 * - the option for returning the results in km and km/sec rather than 00035 * AU and AU/day 00036 * - specification of the filename of the binary ephemeris table 00037 * The basic function pleph() uses the (as-distributed) pre-compiled 00038 * values for these options. 00039 * The compiler option for non-barycentric values was removed, as it was 00040 * in fact unconditionally overridden. 00041 * Removed option for positions only from interp() 00042 * Moved checks for nutations and librations in pleph_ext() to follow 00043 * state() call, where they must logically be in order to guarantee 00044 * valid values 00045 * Modified et to split integer and fractional parts for allegedly greater 00046 * interpolation precision, though it does not seem to make much difference 00047 * Added static variables to (a) ensure that the earth-moon mass ratio 00048 * is properly set when there are repeated calls to pleph_ext within 00049 * a program; and (b) allow for reading from different tables in calls 00050 * from the same calling program 00051 * Replaced several initialization loops with memset and memcpy calls 00052 * 00053 */ 00054 00055 00056 #define DEFAULT_EPHEM_FILE "/home/rick/src/ephem/tables/unxp2000.405" 00057 /* 00058 #define DEFAULT_EPHEM_FILE "/home/soi/CM/tables/jpleph.403.1950-2050" 00059 */ 00060 #include <stdio.h> 00061 #include <math.h> 00062 00063 /* 00064 * Byte reordering function added to support reading of big-endian binary 00065 * tables on little-endian architectures that still use the IEEE-754 00066 * floating point representation 00067 */ 00068 static void bytereorder (void *data, int size, int num) { 00069 int n, s; 00070 unsigned char *c = (unsigned char *)data; 00071 unsigned char t; 00072 00073 for (n=0; n<num; n++, c+=size) 00074 for (s=0; s<(size/2); s++) { 00075 t = *(c + s); 00076 *(c + s) = *(c + size - s - 1); 00077 *(c + size - s - 1) = t; 00078 } 00079 } 00080 00081 /**************************************************************************** 00082 **** split(tt,fr) **** 00083 ***************************************************************************** 00084 **** this subroutine breaks a d.p. number into a d.p. integer **** 00085 **** and a d.p. fractional part. **** 00086 **** **** 00087 **** calling sequence parameters: **** 00088 **** **** 00089 **** tt = d.p. input number **** 00090 **** **** 00091 **** fr = d.p. 2-word output array. **** 00092 **** fr(1) contains integer part **** 00093 **** fr(2) contains fractional part **** 00094 **** **** 00095 **** for negative input numbers, fr(1) contains the next **** 00096 **** more negative integer; fr(2) contains a positive fraction. **** 00097 ****************************************************************************/ 00098 static void split (double tt, double fr[2]) { 00099 fr[1] = modf (tt, &fr[0]); 00100 if (tt >= 0.0 || fr[1] == 0.0) return; 00101 /* make adjustments for negative input number */ 00102 fr[0] = fr[0] - 1.0; 00103 fr[1] = fr[1] + 1.0; 00104 return; 00105 } 00106 00107 /***************************************************************************** 00108 ** interp(buf,t,ncf,ncm,na,pv) ** 00109 ****************************************************************************** 00110 ** ** 00111 ** this subroutine differentiates and interpolates a ** 00112 ** set of chebyshev coefficients to give position and velocity ** 00113 ** ** 00114 ** calling sequence parameters: ** 00115 ** ** 00116 ** input: ** 00117 ** ** 00118 ** buf 1st location of array of d.p. chebyshev coefficients ** 00119 ** of position ** 00120 ** ** 00121 ** tt fractional time in interval covered by coefficients at ** 00122 ** which interpolation is wanted (0 <= tt <= 1) ** 00123 ** ** 00124 ** intrvl length of whole interval in input time units ** 00125 ** ** 00126 ** ncf # of coefficients per component ** 00127 ** ** 00128 ** ncm # of components per set of coefficients ** 00129 ** ** 00130 ** na # of sets of coefficients in full array ** 00131 ** (i.e., # of sub-intervals in full interval) ** 00132 ** ** 00133 ** output: ** 00134 ** ** 00135 ** pv interpolated quantities requested. dimension ** 00136 ** expected is pv(ncm,ifl), dp. ** 00137 ** ** 00138 *****************************************************************************/ 00139 00140 static void interp (double *coef, double tt, double intrvl, int ncf, int ncm, 00141 int na, double posvel[6]) { 00142 static double pc[18], vc[18]; 00143 static double twot = 0.0; 00144 static int first_call = 1, np = 2, nv = 3; 00145 double dna, dt1, tc, temp, temp1, vfac; 00146 int i, j, l, lcfcm, offset; 00147 00148 if (first_call ) { 00149 /* initialize static vectors */ 00150 pc[0] = 1.0; 00151 pc[1] = 0.0; 00152 vc[1] = 1.0; 00153 first_call = 0; 00154 } 00155 /* 00156 entry point. get correct sub-interval number for this set 00157 of coefficients and then get normalized chebyshev time 00158 within that subinterval. 00159 */ 00160 dna = (double)na; 00161 modf (tt, &dt1); 00162 temp = dna * tt; 00163 l = (int)(temp - dt1); 00164 lcfcm = l * ncf * ncm; 00165 /* tc is the normalized chebyshev time (-1 <= tc <= 1) */ 00166 tc = 2.0 * (modf (temp, &temp1) + dt1) - 1.0; 00167 /* 00168 check to see whether chebyshev time has changed, 00169 and compute new polynomial values if it has. 00170 (the element pc[1] is the value of t1[tc] and hence 00171 contains the value of tc on the previous call.) 00172 */ 00173 if (tc != pc[1]) { 00174 np = 2; 00175 nv = 3; 00176 pc[1] = tc; 00177 twot = tc + tc; 00178 } 00179 /* 00180 be sure that at least 'ncf' polynomials have been evaluated 00181 and are stored in the array 'pc'. 00182 */ 00183 if (np < ncf) { 00184 for (i=np; i<ncf; i++) pc[i] = twot*pc[i-1] - pc[i-2]; 00185 np = ncf; 00186 } 00187 /* interpolate to get position for each component */ 00188 for (i=0; i<ncm; i++){ /* ncm is a number of coordinates */ 00189 posvel[i] = 0.0; 00190 offset = i*ncf + lcfcm; 00191 for (j=ncf-1; j>=0; j--) 00192 posvel[i] += pc[j] * coef[j + offset]; 00193 } 00194 /* be sure enough derivative polynomials have been generated and stored */ 00195 /* for velocity interpolation */ 00196 vfac = (dna + dna) / intrvl; 00197 vc[2] = twot + twot; 00198 if (nv < ncf) { 00199 for (i=nv; i<ncf; i++) vc[i] = twot*vc[i-1] + pc[i-1] + pc[i-1] - vc[i-2]; 00200 nv = ncf; 00201 } 00202 /* interpolate to get velocity for each component */ 00203 for (i=0; i<ncm; i++) { 00204 posvel[i+ncm] = 0.0; 00205 offset = i*ncf + lcfcm; 00206 for (j=ncf-1; j>0; j--) 00207 posvel[i+ncm] += vc[j] * coef[j + offset]; 00208 posvel[i+ncm] *= vfac; 00209 } 00210 return; 00211 } 00212 00213 /***************************************************************************** 00214 ** static int state (double *et2, int *list, char *table, int in_kms, 00215 double pv[][6], double *nut, double *emrat, double pvsun[6], 00216 int *nut_flag, int *lib_flag) 00217 ****************************************************************************** 00218 ** This subroutine reads and interpolates the jpl planetary ephemeris file ** 00219 ** ** 00220 ** Calling sequence parameters: ** 00221 ** ** 00222 ** Input: ** 00223 ** ** 00224 ** et2[] double, 2-element JED epoch at which interpolation ** 00225 ** is wanted. Any combination of et2[0]+et2[1] which falls ** 00226 ** within the time span on the file is a permissible epoch. ** 00227 ** ** 00228 ** a. for ease in programming, the user may put the ** 00229 ** entire epoch in et2[0] and set et2[1]=0.0 ** 00230 ** ** 00231 ** b. for maximum interpolation accuracy, set et2[0] = ** 00232 ** the most recent midnight at or before interpolation ** 00233 ** epoch and set et2[1] = fractional part of a day ** 00234 ** elapsed between et2[0] and epoch. ** 00235 ** ** 00236 ** c. as an alternative, it may prove convenient to set ** 00237 ** et2[0] = some fixed epoch, such as start of integration,** 00238 ** and et2[1] = elapsed interval between then and epoch. ** 00239 ** ** 00240 ** list 12-element integer array specifying what interpolation ** 00241 ** is wanted for each of the "bodies" on the file. ** 00242 ** ** 00243 ** list[i]=0, no interpolation for body i ** 00244 ** !=0, position and velocity ** 00245 ** ** 00246 ** The designation of the astronomical bodies by i is: ** 00247 ** ** 00248 ** i = 0: Mercury ** 00249 ** = 1: Venus ** 00250 ** = 2: Earth-Moon barycenter ** 00251 ** = 3: Mars ** 00252 ** = 4: Jupiter ** 00253 ** = 5: Saturn ** 00254 ** = 6: Uranus ** 00255 ** = 7: Neptune ** 00256 ** = 8: Pluto ** 00257 ** = 9: geocentric Moon ** 00258 ** =10: nutations in longitude and obliquity ** 00259 ** =11: lunar librations (if on file) ** 00260 ** ** 00261 * table character string containing name of file containing binary 00262 * ephemeris data 00263 * 00264 * in_kms integer flag: if non-zero, values returned in km and km/s 00265 * otherwise AU and AU/day 00266 * 00267 ** output: ** 00268 ** ** 00269 ** pv[][6] double array that will contain requested interpolated ** 00270 ** quantities. The body specified by list[i] will have its ** 00271 ** state in the array starting at pv[i][0] (on any given ** 00272 ** call, only those words in 'pv' which are affected by the ** 00273 ** first 10 'list' entries (and by list[11] if librations are ** 00274 ** on the file) are set. The rest of the 'pv' array ** 00275 ** is untouched.) The order of components in pv[][] is: ** 00276 ** pv[][0]=x,pv[][1]=y,....pv[][5]=dz. ** 00277 ** ** 00278 ** All output vectors are referenced to the Earth mean ** 00279 ** equator and equinox of epoch. The Moon state is always ** 00280 ** geocentric; the other nine states are either heliocentric ** 00281 ** or solar-system barycentric, depending on the setting of ** 00282 ** global variables (see below). ** 00283 ** ** 00284 ** Lunar librations, if on file, are put into pv[10][k] if ** 00285 ** list[11] is 1 or 2. ** 00286 ** ** 00287 ** nut dp 4-element array that will contain nutations and rates, ** 00288 ** depending on the setting of list[10]. the order of ** 00289 ** quantities in nut is: ** 00290 ** ** 00291 ** d psi (nutation in longitude) ** 00292 ** d epsilon (nutation in obliquity) ** 00293 ** d psi dot ** 00294 ** d epsilon dot ** 00295 * 00296 * emrat pointer to variable containing earth-moon mass ratio 00297 * 00298 * pvsun[6] double 6-element array containing the barycentric position 00299 * and velocity of the Sun 00300 * 00301 * nut_flag flag value telling whether table contains nutations 00302 * 00303 * lib_flag flag value telling whether table contains librations 00304 * 00305 ****************************************************************************/ 00306 static int state (double *et2, int *list, char *table, int in_kms, 00307 double pv[][6], double *nut, double *emrat, double pvsun[6], 00308 int *nut_flag, int *lib_flag) { 00309 static FILE *binfile; 00310 static double *buf; 00311 static double val[400]; 00312 static double au, tstart, tstop, tstep; 00313 static double earth_moon_mass_ratio; 00314 double pefau[6], pjd[4], ss[3]; 00315 double aufac, jda0, jda1, s, t, tstp; 00316 static int ipt[13][3]; 00317 static int first_call = 1, reorder = 0, recl = 0; 00318 static int coeffs, recsize; 00319 int i, j, n, ncon, numde, rec, series_num; 00320 static char *last_table = NULL; 00321 char title[3][84], name[400][6]; 00322 00323 if (last_table) { 00324 if (strcmp (table, last_table)) { 00325 free (last_table); 00326 fclose (binfile); 00327 first_call = 1; 00328 } 00329 } 00330 if (first_call) { 00331 if (sizeof (double) != 8) { 00332 fprintf (stderr, "** Error: Local architecture does not use 8-byte DP floats\n"); 00333 return (1); 00334 } 00335 if (!(binfile = fopen (table, "r"))) { 00336 fprintf (stderr, "** Error: Unable to open file \"%s\" for read\n", table); 00337 return (1); 00338 } 00339 for (n = 0; n < 3; n++) fread (title[n], 84, 1, binfile); 00340 if (sscanf (title[0], "JPL Planetary Ephemeris DE%d/LE", &series_num) != 1) { 00341 fprintf (stderr, "** Error: Unexpected title format in file\n %s\n", 00342 table); 00343 return (1); 00344 } 00345 if (sscanf (title[1], "Start Epoch: JED= %lf", &jda0) != 1) { 00346 fprintf (stderr, "** Error: Unexpected title format in file\n %s\n", 00347 table); 00348 return (1); 00349 } 00350 if (sscanf (title[2], "Final Epoch: JED= %lf", &jda1) != 1) { 00351 fprintf (stderr, "** Error: Unexpected title format in file\n %s\n", 00352 table); 00353 return (1); 00354 } 00355 if (series_num == 200 || series_num == 202) coeffs = 826; 00356 else if (series_num == 403 || series_num == 405) coeffs = 1018; 00357 else if (series_num == 404 || series_num == 406) coeffs = 728; 00358 else { 00359 fprintf (stderr, "** Error: Unknown ephemeris series %d in file\n %s\\n", 00360 series_num, table); 00361 return (1); 00362 } 00363 recsize = coeffs * sizeof (double); 00364 buf = (double *)malloc (recsize); 00365 /* May not need these */ 00366 for (n = 0; n < 400; n++) fread (name[n], 6, 1, binfile); 00367 00368 fread (ss, sizeof (ss), 1, binfile); 00369 if ((fabs (ss[0] - jda0) > 0.5) || (fabs (ss[1] - jda1) > 0.5)) { 00370 bytereorder (ss, sizeof (double), 3); 00371 if ((fabs (ss[0] - jda0) > 0.5) || (fabs (ss[1] - jda1) > 0.5)) { 00372 fprintf (stderr, "Error: Local architecture does not support IEEE-754\n"); 00373 fprintf (stderr, " or binary epoch limits inconsistent with header\n"); 00374 return (1); 00375 } 00376 reorder = 1; 00377 } 00378 tstart = ss[0]; 00379 tstop = ss[1]; 00380 tstep = ss[2]; 00381 /* ncon is 144 in DE403, 156 in DE405 */ 00382 fread (&ncon, sizeof (int), 1, binfile); 00383 fread (&au, sizeof (double), 1, binfile); 00384 if (reorder) bytereorder (&au, sizeof (double), 1); 00385 fread (emrat, sizeof (double), 1, binfile); 00386 if (reorder) bytereorder (emrat, sizeof (double), 1); 00387 earth_moon_mass_ratio = *emrat; 00388 for (n = 0; n < 12; n++) fread (ipt[n], sizeof (int), 3, binfile); 00389 /* Ignore? */ 00390 fread (&numde, sizeof (int), 1, binfile); 00391 fread (ipt[12], sizeof (int), 3, binfile); 00392 if (reorder) bytereorder (ipt, sizeof (int), 39); 00393 *nut_flag = ipt[11][1]; 00394 *lib_flag = ipt[12][1]; 00395 fseek (binfile, recsize, SEEK_SET); 00396 /* May not need these */ 00397 fread (val, sizeof (val), 1, binfile); 00398 if (reorder) bytereorder (val, sizeof (double), 400); 00399 last_table = malloc (strlen (table) + 1); 00400 strcpy (last_table, table); 00401 first_call = 0; 00402 } 00403 00404 /****************************** main entry point *****************************/ 00405 00406 if (et2[0] == 0.0) return 0; 00407 s = et2[0] - 0.5; 00408 split (s, &pjd[0]); 00409 split (et2[1], &pjd[2]); 00410 pjd[0] = pjd[0] + pjd[2] + 0.5; 00411 pjd[1] = pjd[1] + pjd[3]; 00412 split (pjd[1], &pjd[2]); 00413 pjd[0] = pjd[0] + pjd[2]; 00414 00415 /* here pjd[0] contains last midnight before epoch desired (in JED: *.5) 00416 and pjd[3] contains the remaining, fractional part of the epoch */ 00417 /* error return for epoch out of range */ 00418 if ((pjd[0] + pjd[3]) < tstart || (pjd[0] + pjd[3]) > tstop) { 00419 fprintf (stderr, "Requested JED not within ephemeris limits:\n"); 00420 fprintf (stderr, "%f < %f || > %f\n", pjd[0] + pjd[3], tstart, tstop); 00421 return 1; 00422 } 00423 *emrat = earth_moon_mass_ratio; 00424 /* Calculate record # and relative time in interval */ 00425 /* add 2 to adjust for the first two records containing header data */ 00426 rec = (long)((pjd[0] - tstart) / tstep) + 2; 00427 if (pjd[0] == tstop) rec--; 00428 t = (pjd[0] - ((rec - 2) * tstep + tstart) + pjd[3]) / tstep; 00429 /* read correct record if not in core (static vector buf[]) */ 00430 if (rec != recl) { 00431 if (fseek (binfile, rec*recsize, SEEK_SET)) { 00432 fprintf (stderr, "Error seeking to new record within file \"%s\"\n", 00433 table); 00434 return 1; 00435 } 00436 fread (buf, sizeof (*buf), coeffs, binfile); 00437 if (reorder) bytereorder (buf, sizeof (double), coeffs); 00438 recl = rec; 00439 } 00440 tstp = tstep; 00441 aufac = 1.0; 00442 if (in_kms) tstp *= 86400.0; 00443 else aufac /= au; 00444 /* interpolate Solar System barycentric Sun */ 00445 interp (&buf[ipt[10][0]-1], t, tstp, ipt[10][1], 3, ipt[10][2], pefau); 00446 for (n=0; n<6; n++) pvsun[n] = pefau[n] * aufac; 00447 /* check and interpolate whichever bodies are requested */ 00448 for (i=0; i<10; i++) { 00449 if (list[i] == 0) continue; 00450 interp (&buf[ipt[i][0]-1], t, tstp, ipt[i][1], 3, ipt[i][2], pefau); 00451 for (j=0; j<6; j++) pv[i][j] = pefau[j]*aufac; 00452 } 00453 /* do nutations if requested (and if on file) */ 00454 if (list[10] && ipt[11][1] > 0) 00455 interp (&buf[ipt[11][0]-1], t, tstp, ipt[11][1], 2, ipt[11][2], nut); 00456 /* get librations if requested (and if on file) */ 00457 if( list[11] && ipt[12][1] > 0) { 00458 interp (&buf[ipt[12][0]-1], t, tstp, ipt[12][1], 3, ipt[12][2], pefau); 00459 for (j=0; j<6; j++) pv[10][j] = pefau[j]; 00460 } 00461 return 0; 00462 } 00463 00464 /* 00465 * pleph_ext (double jed, int target, int center, int in_kms, double *rrd 00466 * char *table) 00467 * 00468 * Added function to provide extended option capability via additional 00469 * arguments; the traditional function pleph merely calls this function 00470 * with the extra arguments set to their (as-distributed) pre-compiled 00471 * values. 00472 * 00473 * Parameters: 00474 * 00475 * jed = julian ephemeris date at which interpolation is wanted 00476 * 00477 * target = code number of 'target' point 00478 * 00479 * center = code number of center point 00480 * 00481 ** The numbering convention for 'target' and 'center' is: ** 00482 ** ** 00483 ** 1 = mercury 8 = neptune ** 00484 ** 2 = venus 9 = pluto ** 00485 ** 3 = earth 10 = moon ** 00486 ** 4 = mars 11 = sun ** 00487 ** 5 = jupiter 12 = solar-system barycenter ** 00488 ** 6 = saturn 13 = earth-moon barycenter ** 00489 ** 7 = uranus 14 = nutations (longitude and obliq) ** 00490 ** 15 = librations, if on eph. file ** 00491 ** ** 00492 ** (If nutations are wanted, set target = 14. ** 00493 ** For librations, set target = 15. set center= 0) ** 00494 ** 00495 * in_kms units for position velocity array rrd 00496 * 0 -> AU, AU/day; else km, km/sec 00497 * 00498 * rrd output 6-element, double array of position and velocity 00499 * of point 'target' relative to 'center', in units requested. 00500 * For librations the units are radians and radians per day. 00501 * In the case of nutations the first four words of 00502 * rrd will be set to nutations and rates, having units of 00503 * radians and radians/day. 00504 * 00505 */ 00506 void pleph_ext (double jed, int target, int center, int in_kms, double *rrd, 00507 char *table) { 00508 double et2[2], pv[13][6]; /* pv is the position/velocity array 00509 NUMBERED FROM ZERO: 0:Mercury,1:Venus,... 00510 8=Pluto,9=Moon,10=Sun,11=SSBary,12=EMBary 00511 First 10 elements (0-9) are affected by state(), 00512 all are adjusted here. */ 00513 double pvsun[6]; 00514 00515 double emrat; 00516 int i, k, nflag, lflag; 00517 int list[12]; /* list is a vector denoting, for which "body" 00518 ephemeris values should be calculated by state(): 00519 0=Mercury,1=Venus,2=EMBary,...,8=Pluto, 00520 9=geocentric Moon, 10=nutations in long. & obliq. 00521 11= lunar librations */ 00522 00523 /* initialize et2 for 'state' and set up component count */ 00524 et2[0] = trunc (jed); 00525 et2[1] = fmod (jed, 1.0); 00526 memset (rrd, 0, 6 * sizeof (double)); 00527 /* 00528 for (i=0; i<6; i++) rrd[i] = 0.0; 00529 */ 00530 00531 if (target == center) return; 00532 00533 memset (list, 0, 12 * sizeof (int)); 00534 /* 00535 for (i=0; i<12; ++i) list[i] = 0; 00536 */ 00537 /* set up proper entries in 'list' array for state call */ 00538 for (i=0; i<2; i++) { 00539 k = (i) ? center - 1 : target - 1; 00540 if (k <= 9) list[k] = 2; /* Major planets */ 00541 if (k == 9) list[2] = 2; /* for moon state earth state is necessary */ 00542 if (k == 2) list[9] = 2; /* for earth state moon state is necessary */ 00543 if (k == 12) list[2] = 2; /* EMBary state additional */ 00544 } 00545 /* make call to state */ 00546 state (et2, list, table, in_kms, pv, rrd, &emrat, pvsun, &nflag, &lflag); 00547 /* check for nutations */ 00548 if (target == 14) { 00549 if (nflag > 0) { 00550 list[10] = 2; 00551 state (et2, list, table, in_kms, pv, rrd, &emrat, pvsun, &nflag, &lflag); 00552 } 00553 else fprintf (stderr, "***** no nutations on the ephemeris file ******\n"); 00554 return; 00555 } 00556 /* check for librations */ 00557 if (target == 15) { 00558 if (lflag > 0) { 00559 /* there are librations on ephemeris file */ 00560 list[11] = 2; 00561 state (et2, list, table, in_kms, pv, rrd, &emrat, pvsun, &nflag, &lflag); 00562 memcpy (rrd, pv[10], 6 * sizeof (double)); 00563 /* 00564 for (i=0; i<6; i++) rrd[i] = pv[10][i]; 00565 */ 00566 } 00567 else fprintf (stderr, "***** no librations on the ephemeris file ******\n"); 00568 return; 00569 } 00570 /* Solar System barycentric Sun state goes to pv[10][] */ 00571 if (target == 11 || center == 11) 00572 memcpy (pv[10], pvsun, 6 * sizeof (double)); 00573 /* 00574 for (i=0; i<6; i++) pv[10][i] = pvsun[i]; 00575 */ 00576 /* Solar System Barycenter coordinates & velocities equal to zero */ 00577 if (target == 12 || center == 12) 00578 memset (pv[11], 0, 6 * sizeof (double)); 00579 /* 00580 for (i=0; i<6; i++) pv[11][i] = 0.0; 00581 */ 00582 /* Solar System barycentric EMBary state: */ 00583 if (target == 13 || center == 13) 00584 memcpy (pv[12], pv[2], 6 * sizeof (double)); 00585 /* 00586 for (i=0; i<6; i++) pv[12][i] = pv[2][i]; 00587 */ 00588 /* if Moon from Earth or Earth from Moon ..... */ 00589 if ((target * center) == 30 && (target + center) == 13) 00590 memset (pv[2], 0, 6 * sizeof (double)); 00591 /* 00592 for(i=0; i<6; i++) pv[2][i] = 0.0; 00593 */ 00594 else { 00595 if (list[2] == 2) /* calculate earth state from EMBary */ 00596 for (i=0; i<6; i++) pv[2][i] -= pv[9][i] / (1.0 + emrat); 00597 if (list[9] == 2) /* calculate Solar System barycentric Moon state */ 00598 for (i=0; i<6; i++) pv[9][i] += pv[2][i]; 00599 } 00600 00601 for (i=0; i<6; i++) rrd[i] = pv[target-1][i] - pv[center-1][i]; 00602 return; 00603 } 00604 00605 /***************************************************************************** 00606 ** pleph (jed, target, center, rrd) ** 00607 ** ** 00608 ** This subroutine reads the jpl planetary ephemeris ** 00609 ** and gives the position and velocity of the point 'target' ** 00610 ** with respect to 'center'. ** 00611 ** ** 00612 ** Calling sequence parameters: ** 00613 ** ** 00614 ** jed = (double) julian ephemeris date at which interpolation ** 00615 ** is wanted. ** 00616 ** ** 00617 ** target = integer number of 'target' point. ** 00618 ** ** 00619 ** center = integer number of center point. ** 00620 ** ** 00621 ** The numbering convention for 'target' and 'center' is: ** 00622 ** ** 00623 ** 1 = mercury 8 = neptune ** 00624 ** 2 = venus 9 = pluto ** 00625 ** 3 = earth 10 = moon ** 00626 ** 4 = mars 11 = sun ** 00627 ** 5 = jupiter 12 = solar-system barycenter ** 00628 ** 6 = saturn 13 = earth-moon barycenter ** 00629 ** 7 = uranus 14 = nutations (longitude and obliq) ** 00630 ** 15 = librations, if on eph. file ** 00631 ** ** 00632 ** (If nutations are wanted, set target = 14. ** 00633 ** For librations, set target = 15. set center= 0) ** 00634 ** ** 00635 ** rrd = output 6-element, double array of position and velocity ** 00636 ** of point 'target' relative to 'center'. The units are au and ** 00637 ** au/day. For librations the units are radians and radians ** 00638 ** per day. In the case of nutations the first four words of ** 00639 ** rrd will be set to nutations and rates, having units of ** 00640 ** radians and radians/day. ** 00641 ** ** 00642 *****************************************************************************/ 00643 void pleph (double jed, int target, int center, double *rrd) { 00644 pleph_ext (jed, target, center, 0, rrd, DEFAULT_EPHEM_FILE); 00645 }
1.5.7.1