FACT++  1.0
int eraPlan94 ( double  date1,
double  date2,
int  np,
double  pv[2][3] 
)

Definition at line 3 of file plan94.c.

References eraAnpm(), ERFA_D2PI, ERFA_DAS2R, ERFA_DJ00, ERFA_DJM, i, and t.

Referenced by palPlanet(), and t_plan94().

163 {
164 /* Gaussian constant */
165  static const double GK = 0.017202098950;
166 
167 /* Sin and cos of J2000.0 mean obliquity (IAU 1976) */
168  static const double SINEPS = 0.3977771559319137;
169  static const double COSEPS = 0.9174820620691818;
170 
171 /* Maximum number of iterations allowed to solve Kepler's equation */
172  static const int KMAX = 10;
173 
174  int jstat, i, k;
175  double t, da, dl, de, dp, di, dom, dmu, arga, argl, am,
176  ae, dae, ae2, at, r, v, si2, xq, xp, tl, xsw,
177  xcw, xm2, xf, ci2, xms, xmc, xpxq2, x, y, z;
178 
179 /* Planetary inverse masses */
180  static const double amas[] = { 6023600.0, /* Mercury */
181  408523.5, /* Venus */
182  328900.5, /* EMB */
183  3098710.0, /* Mars */
184  1047.355, /* Jupiter */
185  3498.5, /* Saturn */
186  22869.0, /* Uranus */
187  19314.0 }; /* Neptune */
188 
189 /*
190 ** Tables giving the mean Keplerian elements, limited to t^2 terms:
191 **
192 ** a semi-major axis (AU)
193 ** dlm mean longitude (degree and arcsecond)
194 ** e eccentricity
195 ** pi longitude of the perihelion (degree and arcsecond)
196 ** dinc inclination (degree and arcsecond)
197 ** omega longitude of the ascending node (degree and arcsecond)
198 */
199 
200  static const double a[][3] = {
201  { 0.3870983098, 0.0, 0.0 }, /* Mercury */
202  { 0.7233298200, 0.0, 0.0 }, /* Venus */
203  { 1.0000010178, 0.0, 0.0 }, /* EMB */
204  { 1.5236793419, 3e-10, 0.0 }, /* Mars */
205  { 5.2026032092, 19132e-10, -39e-10 }, /* Jupiter */
206  { 9.5549091915, -0.0000213896, 444e-10 }, /* Saturn */
207  { 19.2184460618, -3716e-10, 979e-10 }, /* Uranus */
208  { 30.1103868694, -16635e-10, 686e-10 } /* Neptune */
209  };
210 
211  static const double dlm[][3] = {
212  { 252.25090552, 5381016286.88982, -1.92789 },
213  { 181.97980085, 2106641364.33548, 0.59381 },
214  { 100.46645683, 1295977422.83429, -2.04411 },
215  { 355.43299958, 689050774.93988, 0.94264 },
216  { 34.35151874, 109256603.77991, -30.60378 },
217  { 50.07744430, 43996098.55732, 75.61614 },
218  { 314.05500511, 15424811.93933, -1.75083 },
219  { 304.34866548, 7865503.20744, 0.21103 }
220  };
221 
222  static const double e[][3] = {
223  { 0.2056317526, 0.0002040653, -28349e-10 },
224  { 0.0067719164, -0.0004776521, 98127e-10 },
225  { 0.0167086342, -0.0004203654, -0.0000126734 },
226  { 0.0934006477, 0.0009048438, -80641e-10 },
227  { 0.0484979255, 0.0016322542, -0.0000471366 },
228  { 0.0555481426, -0.0034664062, -0.0000643639 },
229  { 0.0463812221, -0.0002729293, 0.0000078913 },
230  { 0.0094557470, 0.0000603263, 0.0 }
231  };
232 
233  static const double pi[][3] = {
234  { 77.45611904, 5719.11590, -4.83016 },
235  { 131.56370300, 175.48640, -498.48184 },
236  { 102.93734808, 11612.35290, 53.27577 },
237  { 336.06023395, 15980.45908, -62.32800 },
238  { 14.33120687, 7758.75163, 259.95938 },
239  { 93.05723748, 20395.49439, 190.25952 },
240  { 173.00529106, 3215.56238, -34.09288 },
241  { 48.12027554, 1050.71912, 27.39717 }
242  };
243 
244  static const double dinc[][3] = {
245  { 7.00498625, -214.25629, 0.28977 },
246  { 3.39466189, -30.84437, -11.67836 },
247  { 0.0, 469.97289, -3.35053 },
248  { 1.84972648, -293.31722, -8.11830 },
249  { 1.30326698, -71.55890, 11.95297 },
250  { 2.48887878, 91.85195, -17.66225 },
251  { 0.77319689, -60.72723, 1.25759 },
252  { 1.76995259, 8.12333, 0.08135 }
253  };
254 
255  static const double omega[][3] = {
256  { 48.33089304, -4515.21727, -31.79892 },
257  { 76.67992019, -10008.48154, -51.32614 },
258  { 174.87317577, -8679.27034, 15.34191 },
259  { 49.55809321, -10620.90088, -230.57416 },
260  { 100.46440702, 6362.03561, 326.52178 },
261  { 113.66550252, -9240.19942, -66.23743 },
262  { 74.00595701, 2669.15033, 145.93964 },
263  { 131.78405702, -221.94322, -0.78728 }
264  };
265 
266 /* Tables for trigonometric terms to be added to the mean elements of */
267 /* the semi-major axes */
268 
269  static const double kp[][9] = {
270  { 69613, 75645, 88306, 59899, 15746, 71087, 142173, 3086, 0 },
271  { 21863, 32794, 26934, 10931, 26250, 43725, 53867, 28939, 0 },
272  { 16002, 21863, 32004, 10931, 14529, 16368, 15318, 32794, 0 },
273  { 6345, 7818, 15636, 7077, 8184, 14163, 1107, 4872, 0 },
274  { 1760, 1454, 1167, 880, 287, 2640, 19, 2047, 1454 },
275  { 574, 0, 880, 287, 19, 1760, 1167, 306, 574 },
276  { 204, 0, 177, 1265, 4, 385, 200, 208, 204 },
277  { 0, 102, 106, 4, 98, 1367, 487, 204, 0 }
278  };
279 
280  static const double ca[][9] = {
281  { 4, -13, 11, -9, -9, -3, -1, 4, 0 },
282  { -156, 59, -42, 6, 19, -20, -10, -12, 0 },
283  { 64, -152, 62, -8, 32, -41, 19, -11, 0 },
284  { 124, 621, -145, 208, 54, -57, 30, 15, 0 },
285  { -23437, -2634, 6601, 6259, -1507,-1821, 2620, -2115, -1489 },
286  { 62911,-119919, 79336,17814,-24241,12068, 8306, -4893, 8902 },
287  { 389061,-262125,-44088, 8387,-22976,-2093, -615, -9720, 6633 },
288  { -412235,-157046,-31430,37817, -9740, -13, -7449, 9644, 0 }
289  };
290 
291  static const double sa[][9] = {
292  { -29, -1, 9, 6, -6, 5, 4, 0, 0 },
293  { -48, -125, -26, -37, 18, -13, -20, -2, 0 },
294  { -150, -46, 68, 54, 14, 24, -28, 22, 0 },
295  { -621, 532, -694, -20, 192, -94, 71, -73, 0 },
296  { -14614,-19828, -5869, 1881, -4372, -2255, 782, 930, 913 },
297  { 139737, 0, 24667, 51123, -5102, 7429, -4095, -1976, -9566 },
298  { -138081, 0, 37205,-49039,-41901,-33872,-27037,-12474, 18797 },
299  { 0, 28492,133236, 69654, 52322,-49577,-26430, -3593, 0 }
300  };
301 
302 /* Tables giving the trigonometric terms to be added to the mean */
303 /* elements of the mean longitudes */
304 
305  static const double kq[][10] = {
306  { 3086,15746,69613,59899,75645,88306, 12661, 2658, 0, 0 },
307  { 21863,32794,10931, 73, 4387,26934, 1473, 2157, 0, 0 },
308  { 10,16002,21863,10931, 1473,32004, 4387, 73, 0, 0 },
309  { 10, 6345, 7818, 1107,15636, 7077, 8184, 532, 10, 0 },
310  { 19, 1760, 1454, 287, 1167, 880, 574, 2640, 19, 1454 },
311  { 19, 574, 287, 306, 1760, 12, 31, 38, 19, 574 },
312  { 4, 204, 177, 8, 31, 200, 1265, 102, 4, 204 },
313  { 4, 102, 106, 8, 98, 1367, 487, 204, 4, 102 }
314  };
315 
316  static const double cl[][10] = {
317  { 21, -95, -157, 41, -5, 42, 23, 30, 0, 0 },
318  { -160, -313, -235, 60, -74, -76, -27, 34, 0, 0 },
319  { -325, -322, -79, 232, -52, 97, 55, -41, 0, 0 },
320  { 2268, -979, 802, 602, -668, -33, 345, 201, -55, 0 },
321  { 7610, -4997,-7689,-5841,-2617, 1115,-748,-607, 6074, 354 },
322  { -18549, 30125,20012, -730, 824, 23,1289,-352, -14767, -2062 },
323  { -135245,-14594, 4197,-4030,-5630,-2898,2540,-306, 2939, 1986 },
324  { 89948, 2103, 8963, 2695, 3682, 1648, 866,-154, -1963, -283 }
325  };
326 
327  static const double sl[][10] = {
328  { -342, 136, -23, 62, 66, -52, -33, 17, 0, 0 },
329  { 524, -149, -35, 117, 151, 122, -71, -62, 0, 0 },
330  { -105, -137, 258, 35, -116, -88,-112, -80, 0, 0 },
331  { 854, -205, -936, -240, 140, -341, -97, -232, 536, 0 },
332  { -56980, 8016, 1012, 1448,-3024,-3710, 318, 503, 3767, 577 },
333  { 138606,-13478,-4964, 1441,-1319,-1482, 427, 1236, -9167, -1918 },
334  { 71234,-41116, 5334,-4935,-1848, 66, 434, -1748, 3780, -701 },
335  { -47645, 11647, 2166, 3194, 679, 0,-244, -419, -2531, 48 }
336  };
337 
338 /*--------------------------------------------------------------------*/
339 
340 /* Validate the planet number. */
341  if ((np < 1) || (np > 8)) {
342  jstat = -1;
343 
344  /* Reset the result in case of failure. */
345  for (k = 0; k < 2; k++) {
346  for (i = 0; i < 3; i++) {
347  pv[k][i] = 0.0;
348  }
349  }
350 
351  } else {
352 
353  /* Decrement the planet number to start at zero. */
354  np--;
355 
356  /* Time: Julian millennia since J2000.0. */
357  t = ((date1 - ERFA_DJ00) + date2) / ERFA_DJM;
358 
359  /* OK status unless remote date. */
360  jstat = fabs(t) <= 1.0 ? 0 : 1;
361 
362  /* Compute the mean elements. */
363  da = a[np][0] +
364  (a[np][1] +
365  a[np][2] * t) * t;
366  dl = (3600.0 * dlm[np][0] +
367  (dlm[np][1] +
368  dlm[np][2] * t) * t) * ERFA_DAS2R;
369  de = e[np][0] +
370  ( e[np][1] +
371  e[np][2] * t) * t;
372  dp = eraAnpm((3600.0 * pi[np][0] +
373  (pi[np][1] +
374  pi[np][2] * t) * t) * ERFA_DAS2R);
375  di = (3600.0 * dinc[np][0] +
376  (dinc[np][1] +
377  dinc[np][2] * t) * t) * ERFA_DAS2R;
378  dom = eraAnpm((3600.0 * omega[np][0] +
379  (omega[np][1] +
380  omega[np][2] * t) * t) * ERFA_DAS2R);
381 
382  /* Apply the trigonometric terms. */
383  dmu = 0.35953620 * t;
384  for (k = 0; k < 8; k++) {
385  arga = kp[np][k] * dmu;
386  argl = kq[np][k] * dmu;
387  da += (ca[np][k] * cos(arga) +
388  sa[np][k] * sin(arga)) * 1e-7;
389  dl += (cl[np][k] * cos(argl) +
390  sl[np][k] * sin(argl)) * 1e-7;
391  }
392  arga = kp[np][8] * dmu;
393  da += t * (ca[np][8] * cos(arga) +
394  sa[np][8] * sin(arga)) * 1e-7;
395  for (k = 8; k < 10; k++) {
396  argl = kq[np][k] * dmu;
397  dl += t * (cl[np][k] * cos(argl) +
398  sl[np][k] * sin(argl)) * 1e-7;
399  }
400  dl = fmod(dl, ERFA_D2PI);
401 
402  /* Iterative soln. of Kepler's equation to get eccentric anomaly. */
403  am = dl - dp;
404  ae = am + de * sin(am);
405  k = 0;
406  dae = 1.0;
407  while (k < KMAX && fabs(dae) > 1e-12) {
408  dae = (am - ae + de * sin(ae)) / (1.0 - de * cos(ae));
409  ae += dae;
410  k++;
411  if (k == KMAX-1) jstat = 2;
412  }
413 
414  /* True anomaly. */
415  ae2 = ae / 2.0;
416  at = 2.0 * atan2(sqrt((1.0 + de) / (1.0 - de)) * sin(ae2),
417  cos(ae2));
418 
419  /* Distance (AU) and speed (radians per day). */
420  r = da * (1.0 - de * cos(ae));
421  v = GK * sqrt((1.0 + 1.0 / amas[np]) / (da * da * da));
422 
423  si2 = sin(di / 2.0);
424  xq = si2 * cos(dom);
425  xp = si2 * sin(dom);
426  tl = at + dp;
427  xsw = sin(tl);
428  xcw = cos(tl);
429  xm2 = 2.0 * (xp * xcw - xq * xsw);
430  xf = da / sqrt(1 - de * de);
431  ci2 = cos(di / 2.0);
432  xms = (de * sin(dp) + xsw) * xf;
433  xmc = (de * cos(dp) + xcw) * xf;
434  xpxq2 = 2 * xp * xq;
435 
436  /* Position (J2000.0 ecliptic x,y,z in AU). */
437  x = r * (xcw - xm2 * xp);
438  y = r * (xsw + xm2 * xq);
439  z = r * (-xm2 * ci2);
440 
441  /* Rotate to equatorial. */
442  pv[0][0] = x;
443  pv[0][1] = y * COSEPS - z * SINEPS;
444  pv[0][2] = y * SINEPS + z * COSEPS;
445 
446  /* Velocity (J2000.0 ecliptic xdot,ydot,zdot in AU/d). */
447  x = v * (( -1.0 + 2.0 * xp * xp) * xms + xpxq2 * xmc);
448  y = v * (( 1.0 - 2.0 * xq * xq) * xmc - xpxq2 * xms);
449  z = v * (2.0 * ci2 * (xp * xms + xq * xmc));
450 
451  /* Rotate to equatorial. */
452  pv[1][0] = x;
453  pv[1][1] = y * COSEPS - z * SINEPS;
454  pv[1][2] = y * SINEPS + z * COSEPS;
455 
456  }
457 
458 /* Return the status. */
459  return jstat;
460 
461 }
ERFA_DJ00
#define ERFA_DJ00
Definition: erfam.h:87
i
int i
Definition: db_dim_client.c:21
ERFA_DAS2R
#define ERFA_DAS2R
Definition: erfam.h:60
ERFA_DJM
#define ERFA_DJM
Definition: erfam.h:84
ERFA_D2PI
#define ERFA_D2PI
Definition: erfam.h:48
eraAnpm
double eraAnpm(double a)
Definition: anpm.c:3
t
TT t
Definition: test_client.c:26

+ Here is the call graph for this function:

+ Here is the caller graph for this function: