home *** CD-ROM | disk | FTP | other *** search
/ Garbo / Garbo.cdr / pc / source / star.lzh / star.15 < prev    next >
Encoding:
Text File  |  1990-04-06  |  30.7 KB  |  1,124 lines

  1.  
  2. #! /bin/sh
  3. # This is a shell archive.  Remove anything before this line, then unpack
  4. # it by saving it into a file and typing "sh file".  To overwrite existing
  5. # files, type "sh file -c".  You can also feed this as standard input via
  6. # unshar, or by typing "sh <file", e.g..  If this archive is complete, you
  7. # will see the following message at the end:
  8. #        "End of archive 15 (of 32)."
  9. # Contents:  observe/planetcalc.c
  10. PATH=/bin:/usr/bin:/usr/ucb ; export PATH
  11. if test -f 'observe/planetcalc.c' -a "${1}" != "-c" ; then 
  12.   echo shar: Will not clobber existing file \"'observe/planetcalc.c'\"
  13. else
  14. echo shar: Extracting \"'observe/planetcalc.c'\" \(28915 characters\)
  15. sed "s/^X//" >'observe/planetcalc.c' <<'END_OF_FILE'
  16. X/*
  17. X * planetcalc.c
  18. X * planetary positions calculations
  19. X *
  20. X * Copyright (c) 1990 by Craig Counterman. All rights reserved.
  21. X *
  22. X * This software may be redistributed freely, not sold.
  23. X * This copyright notice and disclaimer of warranty must remain
  24. X *    unchanged. 
  25. X *
  26. X * No representation is made about the suitability of this
  27. X * software for any purpose.  It is provided "as is" without express or
  28. X * implied warranty, to the extent permitted by applicable law.
  29. X *
  30. X */
  31. X
  32. X#ifndef  lint
  33. Xstatic char rcsid[] =
  34. X  "$Header: planetcalc.c,v 1.7 90/02/19 17:21:21 ccount Exp $";
  35. X#endif
  36. X
  37. X
  38. X#include <stdio.h>
  39. X#include <math.h>
  40. X#include "observe.h"
  41. X#include "degree.h"
  42. X
  43. Xstatic double M_sun(), M_1(), M_2(), M_4(), M_5(), M_6();
  44. X
  45. X/* Approximate M for each planet 
  46. X   from Meeus, chapter 25 */
  47. X
  48. Xstatic double M_sun(jd)
  49. Xdouble jd;
  50. X{
  51. X  double T;
  52. X
  53. X  T = (jd - 2415020.0)/36525.0;
  54. X
  55. X  return into_range(358.47583 + 35999.04975*T - 0.000150*T*T -0.0000033*T*T*T);
  56. X}
  57. X
  58. Xstatic double M_1(jd)
  59. Xdouble jd;
  60. X{
  61. X  double T;
  62. X
  63. X  T = (jd - 2415020.0)/36525.0;
  64. X
  65. X  return into_range(102.27938 + 149472.51529*T + 0.000007*T*T);
  66. X}
  67. X
  68. Xstatic double M_2(jd)
  69. Xdouble jd;
  70. X{
  71. X  double T;
  72. X
  73. X  T = (jd - 2415020.0)/36525.0;
  74. X
  75. X  return into_range(212.60322 + 58517.80387*T +0.001286*T*T);
  76. X}
  77. X
  78. Xstatic double M_4(jd)
  79. Xdouble jd;
  80. X{
  81. X  double T;
  82. X
  83. X  T = (jd - 2415020.0)/36525.0;
  84. X
  85. X  return into_range(319.51913 + 19139.85475*T + 0.000181*T*T);
  86. X}
  87. X
  88. Xstatic double M_5(jd)
  89. Xdouble jd;
  90. X{
  91. X  double T;
  92. X
  93. X  T = (jd - 2415020.0)/36525.0;
  94. X
  95. X  return into_range(225.32833 + 3034.69202*T - 0.000722*T*T);
  96. X}
  97. X
  98. Xstatic double M_6(jd)
  99. Xdouble jd;
  100. X{
  101. X  double T;
  102. X
  103. X  T = (jd - 2415020.0)/36525.0;
  104. X
  105. X  return into_range(175.46622 +1221.55147*T - 0.000502*T*T);
  106. X}
  107. X
  108. Xstatic char *planet_name[] = {
  109. X  "Mercury",
  110. X  "Venus",
  111. X  "Mars",
  112. X  "Jupiter",
  113. X  "Saturn",
  114. X  "Uranus",
  115. X  "Neptune"
  116. X};
  117. X
  118. Xstatic char *planet_type[] = {
  119. X  "PM",
  120. X  "PV",
  121. X  "Pm",
  122. X  "PJ",
  123. X  "Ps",
  124. X  "PU",
  125. X  "PN"
  126. X};
  127. X
  128. Xstatic char *planet_colr[] = {
  129. X  "r9",
  130. X  "w9",
  131. X  "r9",
  132. X  "w9",
  133. X  "y9",
  134. X  "g9",
  135. X  "c9"
  136. X};
  137. X
  138. Xtypedef struct {
  139. X  double L[4];
  140. X  double a;
  141. X  double e[4];
  142. X  double i[4];
  143. X  double omega[4];
  144. X  double Omega[4];
  145. X  double size_1au;
  146. X  double mag0;
  147. X} pelements;
  148. X
  149. Xpelements peles[] = {
  150. X  {   /* Mercury */
  151. X    {178.179078, 149474.07078, 3.011e-4, 0.0},
  152. X    0.3870986,
  153. X    {.20561421, 2.046e-5, 3e-8,  0.0},
  154. X    {7.002881, 1.8608e-3, -1.83e-5, 0.0},
  155. X    {28.753753, 0.3702806, 0.0001208, 0.0},
  156. X    {47.145944, 1.1852083, 1.739e-4, 0.0},
  157. X    6.74,
  158. X    -0.42
  159. X    },
  160. X  {   /* Venus */
  161. X    {342.767053, 58519.21191, 0.0003097, 0.0},
  162. X    0.7233316,
  163. X    {0.00682069, -0.00004774, 0.000000091, 0.0},
  164. X    {3.393631, 0.0010058, -0.0000010, 0.0},
  165. X    {54.384186, 0.5081861, -0.0013864, 0.0},
  166. X    {75.779647, 0.8998500, 0.0004100, 0.0},
  167. X    16.92,
  168. X    -4.4
  169. X    },
  170. X  {   /* Mars */
  171. X    {293.737334, 19141.69551, 0.0003107, 0.0},
  172. X    1.5236883,
  173. X    {0.09331290, 0.000092064, -0.000000077, 0.0},
  174. X    {1.850333, -0.0006750, 0.0000126, 0.0},
  175. X    {285.431761, 1.0697667, 0.0001313, 0.00000414},
  176. X    {48.786442, 0.7709917, -0.0000014, -0.00000533},
  177. X    9.36,
  178. X    -1.52
  179. X    },
  180. X  {   /* Jupiter */
  181. X    {238.049257, 3036.301986, 0.0003347, -0.00000165},
  182. X    5.202561,
  183. X    {0.04833475, 0.000164180, -0.0000004676, -0.0000000017},
  184. X    {1.308736, -0.0056961, 0.0000039, 0.0},
  185. X    {273.277558, 0.5994317, 0.00070405, 0.00000508},
  186. X    {99.443414, 1.0105300, 0.00035222, -0.00000851},
  187. X    196.74,
  188. X    -9.4
  189. X    },
  190. X  {   /* Saturn */
  191. X    {266.564377, 1223.509884, 0.0003245, -0.0000058},
  192. X    9.554747,
  193. X    {0.05589232, -0.00034550, -0.000000728, 0.00000000074},
  194. X    {2.492519, -0.0039189, -0.00001549, 0.00000004},
  195. X    {338.307800, 1.0852207, 0.00097854, 0.00000992},
  196. X    {112.790414, 0.8731951, -0.00015218, -0.00000531},
  197. X    165.6,
  198. X    -8.88
  199. X    },
  200. X  {   /* Uranus */
  201. X    {244.197470, 429.863546, 0.0003160, -0.00000060},
  202. X    19.21814,
  203. X    {0.0463444, -0.00002658, 0.000000077, 0.0},
  204. X    {0.772464, 0.0006253, 0.0000395, 0.0},
  205. X    {98.071581, 0.9857650, -0.0010745, -0.00000061},
  206. X    {73.477111, 0.4986678, 0.0013117, 0.0},
  207. X    65.8,
  208. X    -7.19
  209. X    },
  210. X  {   /* Neptune */
  211. X    {84.457994, 219.885914, 0.0003205, -0.00000060},
  212. X    30.10957,
  213. X    {0.00899704, 0.000006330, -0.000000002, 0.0},
  214. X    {1.779242, -0.0095436, -0.0000091, 0.0},
  215. X    {276.045975, 0.3256394, 0.00014095, 0.000004113},
  216. X    {130.681389, 1.0989350, 0.00024987, -0.000004718},
  217. X    62.2,
  218. X    -6.87
  219. X    }
  220. X};
  221. X
  222. Xtypedef struct {
  223. X  double alpha_1, delta_1;
  224. X  double W_0, W_dot;
  225. X} rot_els_t;
  226. Xrot_els_t rot_els[] = {
  227. X  {                /* Mercury */
  228. X    280.98, 61.44,
  229. X    142.14, 6.13850
  230. X    },
  231. X  {                /* Venus */
  232. X    272.78, 67.21,
  233. X    353.00, -1.48142
  234. X    },
  235. X  {                /* Mars */
  236. X    317.61, 52.85,
  237. X    237.21, 350.89198
  238. X    },
  239. X  {                /* Jupiter III */
  240. X    268.04, 64.49,
  241. X    156.03, 870.53600
  242. X    },
  243. X  {                /* Saturn III */
  244. X    40.09, 83.49,
  245. X    223.60, 810.79390
  246. X    },
  247. X  {                /* Uranus */
  248. X    257.27, -15.09,
  249. X    287.17, -499.42197
  250. X    },
  251. X  {                /* Neptune */
  252. X    295.24, 40.62,
  253. X    315.33, 468.75000
  254. X    }
  255. X};
  256. X
  257. Xstatic double polynom(jd, a)
  258. X     double jd;
  259. X     double a[4];
  260. X{
  261. X  double T;
  262. X
  263. X  T = (jd - 2415020.0)/36525.0;
  264. X
  265. X  return (a[0] + a[1]*T + a[2]*T*T + a[3]*T*T*T);
  266. X}
  267. X
  268. Xstatic void mercury(), venus(), mars(),
  269. X  jupiter(), saturn(), uranus(), neptune();
  270. X
  271. X
  272. X/* Calculate alpha and delta
  273. X   from lambda and beta (and epsilon)
  274. X   which are from r, Delta, psi, b, l, and Theta (in sun_data)
  275. X   which are from u, i (given), and Omega (given)
  276. X   u is from L (given), nu, and M
  277. X   nu and M are calculated.
  278. X   r is from E, a (given) and e (given)
  279. X   E is calculated
  280. X
  281. X   calculate mag from Delta, size form Delta, phase (== beta).
  282. X */
  283. X
  284. Xvoid planet_pos(jd, sun_data, nplanet, data)
  285. X     double jd;            /* time, jd */
  286. X     sun_data_t sun_data;
  287. X     int nplanet;        /* Planet number, 0 = mercury */
  288. X     planet_data_t *data;
  289. X{
  290. X  double L_, a_, e_, i_, omega_, Omega_, M_;
  291. X  double r;            /* radius distance to sun */
  292. X  double l, b;            /* ecliptical longitude and latitude */
  293. X  double Delta;            /* Distance to earth */
  294. X  double lambda, beta;        /* geocentric longitude and latitude */
  295. X  double alpha, delta;        /* R.A. and dec. both degrees */
  296. X  double alpha2000, delta2000;    /* R.A. and dec. both degrees equin 2000.0 */
  297. X  double psi;            /* elongation */
  298. X  double N, D;            /* temporary variables */
  299. X  double Theta;            /* Theta of the sun */
  300. X  double epsilon;        /* obliquity */
  301. X  double Cen;            /* center */
  302. X  double A, B;            /* used in calculating p_n */
  303. X  double W_1, K;        /* Used in calculating lambda_e */
  304. X
  305. X  L_ = into_range(polynom(jd, peles[nplanet].L));
  306. X  a_ = peles[nplanet].a;
  307. X  e_ = polynom(jd, peles[nplanet].e);
  308. X  i_ = polynom(jd, peles[nplanet].i);
  309. X  omega_ = into_range(polynom(jd, peles[nplanet].omega));
  310. X  Omega_ = into_range(polynom(jd, peles[nplanet].Omega));
  311. X  M_ = into_range(L_ - omega_ - Omega_);
  312. X
  313. X  /* Perturb */
  314. X  switch (nplanet) {
  315. X  case 0:
  316. X    mercury(jd, L_, a_, e_, i_, omega_, Omega_, M_, &r, &l, &b, &Cen);
  317. X    break;
  318. X  case 1:
  319. X    venus(jd, L_, a_, e_, i_, omega_, Omega_, M_, &r, &l, &b, &Cen);
  320. X    break;
  321. X  case 2:
  322. X    mars(jd, L_, a_, e_, i_, omega_, Omega_, M_, &r, &l, &b, &Cen);
  323. X    break;
  324. X  case 3:
  325. X    jupiter(jd, L_, a_, e_, i_, omega_, Omega_, M_, &r, &l, &b, &Cen);
  326. X    break;
  327. X  case 4:
  328. X    saturn(jd, L_, a_, e_, i_, omega_, Omega_, M_, &r, &l, &b, &Cen);
  329. X    break;
  330. X  case 5:
  331. X    uranus(jd, L_, a_, e_, i_, omega_, Omega_, M_, &r, &l, &b, &Cen);
  332. X    break;
  333. X  case 6:
  334. X    neptune(jd, L_, a_, e_, i_, omega_, Omega_, M_, &r, &l, &b, &Cen);
  335. X    break;
  336. X  default:
  337. X    break;
  338. X  };
  339. X
  340. X  Theta = sun_data.Theta;
  341. X  N = r * DCOS(b) * DSIN(l - Theta);
  342. X  D = r * DCOS(b) * DCOS(l - Theta) + sun_data.R;
  343. X
  344. X  epsilon = obl_jd(jd);
  345. X
  346. X  lambda = into_range(RAD_TO_DEG * atan2(N, D)) + Theta;
  347. X  Delta = sqrt(N*N + D*D + (r * DSIN(b))*(r * DSIN(b)));
  348. X  beta = RAD_TO_DEG * asin(r * DSIN(b) / Delta);
  349. X  psi = RAD_TO_DEG * acos(DCOS(beta) * DCOS(lambda - Theta));
  350. X  if (into_range(lambda - Theta) > 180.0)
  351. X    psi = -psi;
  352. X  alpha = RAD_TO_DEG * atan2(DSIN(lambda)*DCOS(epsilon)
  353. X                 - DTAN(beta) * DSIN(epsilon),
  354. X                 DCOS(lambda));
  355. X  delta = RAD_TO_DEG * asin(DSIN(beta)*DCOS(epsilon)
  356. X                + DCOS(beta)*DSIN(epsilon)*DSIN(lambda));
  357. X  alpha = into_range(alpha);
  358. X
  359. X/* should correct for nutation and aberration */
  360. X  data->name = planet_name[nplanet];
  361. X  data->type = planet_type[nplanet];
  362. X  data->color = planet_colr[nplanet];
  363. X  data->alpha = alpha;
  364. X  data->delta = delta;
  365. X  precess(2000.0 - (2451545.0 - jd) / 365.25,
  366. X      2000.0, alpha, delta, &alpha2000, &delta2000);
  367. X  data->alpha2000 = alpha2000;
  368. X  data->delta2000 = delta2000;
  369. X  data->l = l;
  370. X  data->b = b;
  371. X  data->lambda = lambda;
  372. X  data->beta = beta;
  373. X  data->psi = psi;
  374. X  data->phase =
  375. X    DACOS((r*r + Delta*Delta - sun_data.R*sun_data.R) / (2*r*Delta));
  376. X  if (psi < 0) data->phase = -data->phase;
  377. X  data->r = r;
  378. X  data->Delta = Delta;
  379. X  data->illum_frac = ((r+Delta)*(r+Delta) - sun_data.R*sun_data.R)/(4*r*Delta);
  380. X  data->chi = /* position angle of bright limb */
  381. X    DATAN2(DCOS(sun_data.delta) * DSIN(sun_data.alpha - alpha),
  382. X      DCOS(delta) * DSIN(sun_data.delta)
  383. X      - DSIN(delta) * DCOS(sun_data.delta) * DCOS(sun_data.alpha - alpha));
  384. X  data->Cen = Cen;
  385. X
  386. X  data->mag = 5.0 * log10(r*Delta)
  387. X    - 2.5 * log10(data->illum_frac)
  388. X    + peles[nplanet].mag0;
  389. X
  390. X  data->size = peles[nplanet].size_1au / Delta;
  391. X
  392. X  data->rotation_elements.beta_e =
  393. X    DASIN(-DSIN(rot_els[nplanet].delta_1)*DSIN(delta)
  394. X      - DCOS(rot_els[nplanet].delta_1) * DCOS(delta)
  395. X          * DCOS(rot_els[nplanet].alpha_1 - alpha));
  396. X
  397. X  A = DCOS(rot_els[nplanet].delta_1) * DSIN(rot_els[nplanet].alpha_1 - alpha);
  398. X  A /= DCOS(data->rotation_elements.beta_e);
  399. X  B = DSIN(rot_els[nplanet].delta_1) * DCOS(delta)
  400. X    - DCOS(rot_els[nplanet].delta_1) * DSIN(delta)
  401. X      * DCOS(rot_els[nplanet].alpha_1 - alpha);
  402. X  B /= DCOS(data->rotation_elements.beta_e);
  403. X  data->rotation_elements.p_n = DATAN2(A, B);
  404. X
  405. X  W_1 = into_range(rot_els[nplanet].W_0
  406. X           + rot_els[nplanet].W_dot
  407. X             * (jd - 2447526.5 - 0.0057755 * Delta));
  408. X
  409. X  A = -DCOS(rot_els[nplanet].delta_1) * DSIN(delta)
  410. X    + DSIN(rot_els[nplanet].delta_1) * DCOS(delta)
  411. X      * DCOS(rot_els[nplanet].alpha_1 - alpha);
  412. X  A /= DCOS(data->rotation_elements.beta_e);
  413. X  B = DCOS(delta) * DSIN(rot_els[nplanet].alpha_1 - alpha);
  414. X  B /= DCOS(data->rotation_elements.beta_e);
  415. X  K = DATAN2(A, B);
  416. X  if (rot_els[nplanet].W_dot < 0.0)
  417. X    data->rotation_elements.lambda_e = into_range(K - W_1);
  418. X  else
  419. X    data->rotation_elements.lambda_e = into_range(W_1 - K);
  420. X}
  421. X
  422. X
  423. Xstatic void mercury(jd, L_, a_, e_, i_, omega_, Omega_, M_, r_p, l_p, b_p, C_p)
  424. Xdouble jd, L_, a_, e_, i_, omega_, Omega_, M_;
  425. Xdouble *r_p, *l_p, *b_p, *C_p;
  426. X{
  427. X  double E, nu;
  428. X  double M1, M2, M4, M5, M6;
  429. X  double r, l, b;
  430. X
  431. X  double u;            /* argument of latitude */
  432. X  double r_pert, l_pert;
  433. X
  434. X  M1 = M_1(jd);
  435. X  M2 = M_2(jd);
  436. X  M4 = M_4(jd);
  437. X  M5 = M_5(jd);
  438. X  M6 = M_6(jd);
  439. X
  440. X
  441. X  /* Calculate E and nu */
  442. X  anom_calc(M_, e_, &E, &nu);
  443. X  r = a_ * (1 - e_ * DCOS(E));
  444. X
  445. X  u = L_ + nu - M_ - Omega_;
  446. X  *C_p = nu - M_;
  447. X  l = into_range(RAD_TO_DEG * atan2(DCOS(i_) * DSIN(u), DCOS(u)) + Omega_);
  448. X  b = RAD_TO_DEG * asin(DSIN(u)*DSIN(i_));
  449. X
  450. X  /* Perturbations */
  451. X  l_pert = 0.00204*DCOS(5*M2-2*M1+12.220)
  452. X    +0.00103*DCOS(2*M2-M1-160.692)
  453. X      +0.00091*DCOS(2*M5-M1-37.003)
  454. X    +0.00078*DCOS(5*M2-3*M1+10.137);
  455. X  r_pert = 0.000007525*DCOS(2*M5-M1+53.013)
  456. X    +0.000006802*DCOS(5*M2-3*M1-259.918)
  457. X      +0.000005457*DCOS(2*M2-2*M1-71.188)
  458. X    +0.000003569*DCOS(5*M2-M1-77.75);
  459. X
  460. X  *r_p = r + r_pert;
  461. X  *l_p = l + l_pert;
  462. X  *b_p = b;
  463. X}
  464. X
  465. Xstatic void venus(jd, L_, a_, e_, i_, omega_, Omega_, M_, r_p, l_p, b_p, C_p)
  466. Xdouble jd, L_, a_, e_, i_, omega_, Omega_, M_;
  467. Xdouble *r_p, *l_p, *b_p, *C_p;
  468. X{
  469. X  double E, nu;
  470. X  double M, M1, M2, M4, M5, M6;
  471. X  double r, l, b;
  472. X
  473. X  double u;            /* argument of latitude */
  474. X  double T;
  475. X  double r_pert, l_pert;
  476. X
  477. X  M = M_sun(jd);
  478. X  M1 = M_1(jd);
  479. X  M2 = M_2(jd);
  480. X  M4 = M_4(jd);
  481. X  M5 = M_5(jd);
  482. X  M6 = M_6(jd);
  483. X
  484. X
  485. X  /* Long term perturbation */
  486. X  T = (jd - 2415020.0)/36525.0;
  487. X  l_pert = 0.00077 * DSIN(237.24 + 150.27*T);
  488. X  M_ += l_pert;
  489. X  L_ += l_pert;
  490. X
  491. X  /* Calculate E and nu */
  492. X  anom_calc(M_, e_, &E, &nu);
  493. X  r = a_ * (1 - e_ * DCOS(E));
  494. X
  495. X
  496. X  u = L_ + nu - M_ - Omega_;
  497. X  *C_p = nu - M_;
  498. X  l = into_range(RAD_TO_DEG * atan2(DCOS(i_) * DSIN(u), DCOS(u)) + Omega_);
  499. X  b = RAD_TO_DEG * asin(DSIN(u)*DSIN(i_));
  500. X
  501. X  /* Perturbations */
  502. X  l_pert = 0.00313*DCOS(2*M-2*M2 -148.225)
  503. X    +0.00198*DCOS(3*M-3*M2 +2.565)
  504. X      +0.00136*DCOS(M-M2-119.107)
  505. X    +0.00096*DCOS(3*M-2*M2-135.912)
  506. X      +0.00082*DCOS(M5-M2-208.087);
  507. X  r_pert = 0.000022501 * DCOS(2*M-2*M2-58.208)
  508. X    +0.000019045 * DCOS(3*M-3*M2+92.577)
  509. X      +0.000006887 * DCOS(M5-M2-118.090)
  510. X    +0.000005172 * DCOS(M-M2-29.110)
  511. X      +0.000003620 * DCOS(5*M-4*M2-104.208)
  512. X        +0.000003283 * DCOS(4*M-4*M2+63.513)
  513. X          +0.000003074 * DCOS(2*M5-2*M2-55.167);
  514. X
  515. X  *r_p = r + r_pert;
  516. X  *l_p = l + l_pert;
  517. X  *b_p = b;
  518. X}
  519. X
  520. Xstatic void mars(jd, L_, a_, e_, i_, omega_, Omega_, M_, r_p, l_p, b_p, C_p)
  521. Xdouble jd, L_, a_, e_, i_, omega_, Omega_, M_;
  522. Xdouble *r_p, *l_p, *b_p, *C_p;
  523. X{
  524. X  double E, nu;
  525. X  double M, M1, M2, M4, M5, M6;
  526. X  double r, l, b;
  527. X
  528. X  double u;            /* argument of latitude */
  529. X  double T;
  530. X  double r_pert, l_pert;
  531. X
  532. X  M = M_sun(jd);
  533. X  M1 = M_1(jd);
  534. X  M2 = M_2(jd);
  535. X  M4 = M_4(jd);
  536. X  M5 = M_5(jd);
  537. X  M6 = M_6(jd);
  538. X
  539. X
  540. X  /* Long term perturbation */
  541. X  T = (jd - 2415020.0)/36525.0;
  542. X  l_pert = -0.01133*DSIN(3*M5-8*M4 +4*M) - 0.00933*DCOS(3*M5-8*M4 +4*M);
  543. X
  544. X  M_ += l_pert;
  545. X  L_ += l_pert;
  546. X
  547. X  /* Calculate E and nu */
  548. X  anom_calc(M_, e_, &E, &nu);
  549. X  r = a_ * (1 - e_ * DCOS(E));
  550. X
  551. X
  552. X  u = L_ + nu - M_ - Omega_;
  553. X  *C_p = nu - M_;
  554. X  l = into_range(RAD_TO_DEG * atan2(DCOS(i_) * DSIN(u), DCOS(u)) + Omega_);
  555. X  b = RAD_TO_DEG * asin(DSIN(u)*DSIN(i_));
  556. X
  557. X  /* Perturbations */
  558. X  l_pert = 0.00705*DCOS(M5-M4-48.958)
  559. X    +0.00607*DCOS(2*M5-M4-188.350)
  560. X      +0.00445*DCOS(2*M5-2*M4-191.897)
  561. X    +0.00388*DCOS(M-2*M4+20.495)
  562. X      +0.00238*DCOS(M-M4+35.097)
  563. X        +0.00204*DCOS(2*M-3*M4+158.638)
  564. X          +0.00177*DCOS(3*M4-M2-57.602)
  565. X        +0.00136*DCOS(2*M-4*M4+154.093)
  566. X          +0.00104*DCOS(M5+17.618);
  567. X  r_pert = 0.000053227*DCOS(M5-M4+41.1306)
  568. X    +0.000050989*DCOS(2*M5-2*M4-101.9847)
  569. X      +0.000038278*DCOS(2*M5-M4-98.3292)
  570. X    +0.000015996*DCOS(M-M4-55.555)
  571. X      +0.000014764*DCOS(2*M-3*M4+68.622)
  572. X        +0.000008966*DCOS(M5-2*M4+43.615);
  573. X  r_pert += 0.000007914*DCOS(3*M5-2*M4-139.737)
  574. X    +0.000007004*DCOS(2*M5-3*M4-102.888)
  575. X      +0.000006620*DCOS(M-2*M4+113.202)
  576. X    +0.000004930*DCOS(3*M5-3*M4-76.243)
  577. X      +0.000004693*DCOS(3*M-5*M4+190.603)
  578. X        +0.000004571*DCOS(2*M-4*M4+244.702)
  579. X          +0.000004409*DCOS(3*M5-M4-115.828);
  580. X
  581. X  *r_p = r + r_pert;
  582. X  *l_p = l + l_pert;
  583. X  *b_p = b;
  584. X}
  585. X
  586. X
  587. Xstatic void jupiter(jd, L_, a_, e_, i_, omega_, Omega_, M_, r_p, l_p, b_p, C_p)
  588. Xdouble jd, L_, a_, e_, i_, omega_, Omega_, M_;
  589. Xdouble *r_p, *l_p, *b_p, *C_p;
  590. X{
  591. X  double E, nu;
  592. X  double M, M1, M2, M4, M5, M6;
  593. X  double r, l, b;
  594. X
  595. X  double u;            /* argument of latitude */
  596. X  double T;
  597. X  double A, B, e_pert, a_pert, v, zeta;
  598. X  double P, Q, S, V, W;
  599. X
  600. X  M = M_sun(jd);
  601. X  M1 = M_1(jd);
  602. X  M2 = M_2(jd);
  603. X  M4 = M_4(jd);
  604. X  M5 = M_5(jd);
  605. X  M6 = M_6(jd);
  606. X
  607. X  /* for perturbations */
  608. X  T = (jd - 2415020.0)/36525.0;
  609. X
  610. X  v = T/5.0 + 0.1;
  611. X  P = 237.47555 +3034.9061*T;
  612. X  Q = 265.91650 + 1222.1139*T;
  613. X  S = 243.51721 + 428.4677*T;
  614. X  V = 5.0*Q -2.0*P;
  615. X  W = 2.0*P - 6.0*Q +3.0*S;
  616. X  zeta = Q - P;
  617. X
  618. X
  619. X  A = (0.331364 - 0.010281*v - 0.004692*v*v)*DSIN(V)
  620. X    +(0.003228 - 0.064436*v + 0.002075*v*v)*DCOS(V)
  621. X      -(0.003083 + 0.000275*v - 0.000489*v*v)*DSIN(2*V)
  622. X    +0.002472*DSIN(W)
  623. X      +0.013619*DSIN(zeta)
  624. X        +0.018472*DSIN(2*zeta)
  625. X          +0.006717*DSIN(3*zeta)
  626. X        +0.002775*DSIN(4*zeta)
  627. X          +(0.007275 - 0.001253*v)*DSIN(zeta)*DSIN(Q)
  628. X            +0.006417*DSIN(2*zeta)*DSIN(Q)
  629. X              +0.002439*DSIN(3*zeta)*DSIN(Q);
  630. X  A += -(0.033839 + 0.001125*v)*DCOS(zeta)*DSIN(Q)
  631. X    -0.003767*DCOS(2*zeta)*DSIN(Q)
  632. X      -(0.035681 + 0.001208*v)*DSIN(zeta)*DCOS(Q)
  633. X    -0.004261*DSIN(2*zeta)*DCOS(Q)
  634. X      +0.002178*DCOS(Q)
  635. X        +(-0.006333 + 0.001161*v)*DCOS(zeta)*DCOS(Q)
  636. X          -0.006675*DCOS(2*zeta)*DCOS(Q)
  637. X        -0.002664*DCOS(3*zeta)*DCOS(Q)
  638. X          -0.002572*DSIN(zeta)*DSIN(2*Q)
  639. X            -0.003567*DSIN(2*zeta)*DSIN(2*Q)
  640. X              +0.002094*DCOS(zeta)*DCOS(2*Q)
  641. X            +0.003342*DCOS(2*zeta)*DCOS(2*Q);
  642. X
  643. X  e_pert = (.0003606 + .0000130*v - .0000043*v*v)*DSIN(V)
  644. X    +(.0001289 - .0000580*v)*DCOS(V)
  645. X      -.0006764*DSIN(zeta)*DSIN(Q)
  646. X    -.0001110*DSIN(2*zeta)*DSIN(Q)
  647. X      -.0000224*DSIN(3*zeta)*DSIN(Q)
  648. X        -.0000204*DSIN(Q)
  649. X          +(.0001284 + .0000116*v)*DCOS(zeta)*DSIN(Q)
  650. X        +.0000188*DCOS(2*zeta)*DSIN(Q)
  651. X          +(.0001460 + .0000130*v)*DSIN(zeta)*DCOS(Q)
  652. X            +.0000224*DSIN(2*zeta)*DCOS(Q)
  653. X              -.0000817*DCOS(Q);
  654. X
  655. X  e_pert += .0006074*DCOS(zeta)*DCOS(Q)
  656. X    +.0000992*DCOS(2*zeta)*DCOS(Q)
  657. X      +.0000508*DCOS(3*zeta)*DCOS(Q)
  658. X    +.0000230*DCOS(4*zeta)*DCOS(Q)
  659. X      +.0000108*DCOS(5*zeta)*DCOS(Q);
  660. X
  661. X  e_pert += -(.0000956 + .0000073*v)*DSIN(zeta)*DSIN(2*Q)
  662. X    +.0000448*DSIN(2*zeta)*DSIN(2*Q)
  663. X      +.0000137*DSIN(3*zeta)*DSIN(2*Q)
  664. X    +(-.0000997 + .0000108*v)*DCOS(zeta)*DSIN(2*Q)
  665. X      +.0000480*DCOS(2*zeta)*DSIN(2*Q);
  666. X
  667. X  e_pert += .0000148*DCOS(3*zeta)*DSIN(2*Q)
  668. X    +(-.0000956 +.0000099*v)*DSIN(zeta)*DCOS(2*Q)
  669. X      +.0000490*DSIN(2*zeta)*DCOS(2*Q)
  670. X    +.0000158*DSIN(3*zeta)*DCOS(2*Q);
  671. X
  672. X  e_pert += .0000179*DCOS(2*Q)
  673. X    +(.0001024 + .0000075*v)*DCOS(zeta)*DCOS(2*Q)
  674. X      -.0000437*DCOS(2*zeta)*DCOS(2*Q)
  675. X    -.0000132*DCOS(3*zeta)*DCOS(2*Q);
  676. X
  677. X  B = (0.007192 - 0.003147*v)*DSIN(V)
  678. X    +(-0.020428 - 0.000675*v + 0.000197*v*v)*DCOS(V)
  679. X      +(0.007269 + 0.000672*v)*DSIN(zeta)*DSIN(Q)
  680. X    -0.004344*DSIN(Q)
  681. X      +0.034036*DCOS(zeta)*DSIN(Q)
  682. X        +0.005614*DCOS(2*zeta)*DSIN(Q)
  683. X          +0.002964*DCOS(3*zeta)*DSIN(Q)
  684. X        +0.037761*DSIN(zeta)*DCOS(Q);
  685. X
  686. X  B += 0.006158*DSIN(2*zeta)*DCOS(Q)
  687. X      -0.006603*DCOS(zeta)*DCOS(Q)
  688. X    -0.005356*DSIN(zeta)*DSIN(2*Q)
  689. X      +0.002722*DSIN(2*zeta)*DSIN(2*Q)
  690. X        +0.004483*DCOS(zeta)*DSIN(2*Q);
  691. X
  692. X  B += -0.002642*DCOS(2*zeta)*DSIN(2*Q)
  693. X    +0.004403*DSIN(zeta)*DCOS(2*Q)
  694. X      -0.002536*DSIN(2*zeta)*DCOS(2*Q)
  695. X    +0.005547*DCOS(zeta)*DCOS(2*Q)
  696. X      -0.002689*DCOS(2*zeta)*DCOS(2*Q);
  697. X
  698. X  a_pert = -.000263*DCOS(V)
  699. X    +.000205*DCOS(zeta)
  700. X      +.000693*DCOS(2*zeta)
  701. X    +.000312*DCOS(3*zeta)
  702. X      +.000147*DCOS(4*zeta)
  703. X        +.000299*DSIN(zeta)*DSIN(Q)
  704. X          +.000181*DCOS(2*zeta)*DSIN(Q)
  705. X        +.000204*DSIN(2*zeta)*DCOS(Q)
  706. X          +.000111*DSIN(3*zeta)*DCOS(Q)
  707. X            -.000337*DCOS(zeta)*DCOS(Q)
  708. X              -.000111*DCOS(2*zeta)*DCOS(Q);
  709. X
  710. X  L_ += A;
  711. X  M_ += A - B / e_;
  712. X  e_ += e_pert;
  713. X  a_ += a_pert;
  714. X  omega_ += B;
  715. X
  716. X  /* Calculate E and nu */
  717. X  anom_calc(M_, e_, &E, &nu);
  718. X  r = a_ * (1 - e_ * DCOS(E));
  719. X
  720. X
  721. X  u = L_ + nu - M_ - Omega_;
  722. X  *C_p = nu - M_;
  723. X  l = into_range(RAD_TO_DEG * atan2(DCOS(i_) * DSIN(u), DCOS(u)) + Omega_);
  724. X  b = RAD_TO_DEG * asin(DSIN(u)*DSIN(i_));
  725. X
  726. X  *r_p = r;
  727. X  *l_p = l;
  728. X  *b_p = b;
  729. X}
  730. X
  731. X
  732. Xstatic void saturn(jd, L_, a_, e_, i_, omega_, Omega_, M_, r_p, l_p, b_p, C_p)
  733. Xdouble jd, L_, a_, e_, i_, omega_, Omega_, M_;
  734. Xdouble *r_p, *l_p, *b_p, *C_p;
  735. X{
  736. X  double E, nu;
  737. X  double M, M1, M2, M4, M5, M6;
  738. X  double r, l, b;
  739. X
  740. X  double u;            /* argument of latitude */
  741. X  double T;
  742. X  double A, B, e_pert, a_pert, b_pert, v, zeta, psi;
  743. X  double P, Q, S, V, W;
  744. X
  745. X  M = M_sun(jd);
  746. X  M1 = M_1(jd);
  747. X  M2 = M_2(jd);
  748. X  M4 = M_4(jd);
  749. X  M5 = M_5(jd);
  750. X  M6 = M_6(jd);
  751. X
  752. X  /* for perturbations */
  753. X  T = (jd - 2415020.0)/36525.0;
  754. X
  755. X  v = T/5.0 + 0.1;
  756. X  P = 237.47555 +3034.9061*T;
  757. X  Q = 265.91650 + 1222.1139*T;
  758. X  S = 243.51721 + 428.4677*T;
  759. X  V = 5.0*Q -2.0*P;
  760. X  W = 2.0*P - 6.0*Q +3.0*S;
  761. X  zeta = Q - P;
  762. X  psi = S - Q;
  763. X
  764. X  A = (-0.814181 + 0.018150*v + 0.016714*v*v)*DSIN(V)
  765. X    +(-0.010497 + 0.160906*v - 0.004100*v*v)*DCOS(V)
  766. X      +0.007581*DSIN(2*V)
  767. X    -0.007986*DSIN(W)
  768. X      -0.148811*DSIN(zeta)
  769. X        -0.040786*DSIN(2*zeta)
  770. X          -0.015208*DSIN(3*zeta)
  771. X        -0.006339*DSIN(4*zeta)
  772. X          -0.006244*DSIN(Q);
  773. X  A += (0.008931 + 0.002728*v)*DSIN(zeta)*DSIN(Q)
  774. X      -0.016500*DSIN(2*zeta)*DSIN(Q)
  775. X    -0.005775*DSIN(3*zeta)*DSIN(Q)
  776. X      +(0.081344 + 0.003206*v)*DCOS(zeta)*DSIN(Q)
  777. X        +0.015019*DCOS(2*zeta)*DSIN(Q)
  778. X          +(0.085581 + 0.002494*v)*DSIN(zeta)*DCOS(Q)
  779. X        +(0.025328 - 0.003117*v)*DCOS(zeta)*DCOS(Q);
  780. X  A += 0.014394*DCOS(2*zeta)*DCOS(Q)
  781. X      +0.006319*DCOS(3*zeta)*DCOS(Q)
  782. X    +0.006369*DSIN(zeta)*DSIN(2*Q)
  783. X      +0.009156*DSIN(2*zeta)*DSIN(2*Q)
  784. X        +0.007525*DSIN(3*psi)*DSIN(2*Q)
  785. X          -0.005236*DCOS(zeta)*DCOS(2*Q)
  786. X        -0.007736*DCOS(2*zeta)*DCOS(2*Q)
  787. X          -0.007528*DCOS(3*psi)*DCOS(2*Q);
  788. X
  789. X  e_pert = (-.0007927 + .0002548*v +.0000091*v*v)*DSIN(V)
  790. X    +(.0013381 + .0001226*v -.0000253*v*v)*DCOS(V)
  791. X      +(.0000248 - .0000121*v)*DSIN(2*V)
  792. X    -(.0000305 + .0000091*v)*DCOS(2*V)
  793. X      +.0000412*DSIN(2*zeta)
  794. X        +.0012415*DSIN(Q)
  795. X          +(.0000390 -.0000617*v)*DSIN(zeta)*DSIN(Q)
  796. X        +(.0000165 - .0000204*v)*DSIN(2*zeta)*DSIN(Q)
  797. X          +.0026599*DCOS(zeta)*DSIN(Q)
  798. X            -.0004687*DCOS(2*zeta)*DSIN(Q);
  799. X  e_pert += -.0001870*DCOS(3*zeta)*DSIN(Q)
  800. X      -.0000821*DCOS(4*zeta)*DSIN(Q)
  801. X    -.0000377*DCOS(5*zeta)*DSIN(Q)
  802. X      +.0000497*DCOS(2*psi)*DSIN(Q)
  803. X        +(.0000163 - .0000611*v)*DCOS(Q)
  804. X          -.0012696*DSIN(zeta)*DCOS(Q)
  805. X        -.0004200*DSIN(2*zeta)*DCOS(Q)
  806. X          -.0001503*DSIN(3*zeta)*DCOS(Q)
  807. X            -.0000619*DSIN(4*zeta)*DCOS(Q)
  808. X              -.0000268*DSIN(5*zeta)*DCOS(Q);
  809. X  e_pert += -(.0000282 + .0001306*v)*DCOS(zeta)*DCOS(Q)
  810. X      +(-.0000086 + .0000230*v)*DCOS(2*zeta)*DCOS(Q)
  811. X    +.0000461*DSIN(2*psi)*DCOS(Q)
  812. X      -.0000350*DSIN(2*Q)
  813. X        +(.0002211 - .0000286*v)*DSIN(zeta)*DSIN(2*Q)
  814. X          -.0002208*DSIN(2*zeta)*DSIN(2*Q)
  815. X        -.0000568*DSIN(3*zeta)*DSIN(2*Q)
  816. X          -.0000346*DSIN(4*zeta)*DSIN(2*Q)
  817. X            -(.0002780 + .0000222*v)*DCOS(zeta)*DSIN(2*Q)
  818. X              +(.0002022 + .0000263*v)*DCOS(2*zeta)*DSIN(2*Q);
  819. X  e_pert += .0000248*DCOS(3*zeta)*DSIN(2*Q)
  820. X      +.0000242*DSIN(3*psi)*DSIN(2*Q)
  821. X    +.0000467*DCOS(3*psi)*DSIN(2*Q)
  822. X      -.0000490*DCOS(2*Q)
  823. X        -(.0002842 + .0000279*v)*DSIN(zeta)*DCOS(2*Q)
  824. X          +(.0000128 + .0000226*v)*DSIN(2*zeta)*DCOS(2*Q)
  825. X        +.0000224*DSIN(3*zeta)*DCOS(2*Q)
  826. X          +(-.0001594 + .0000282*v)*DCOS(zeta)*DCOS(2*Q)
  827. X            +(.0002162 - .0000207*v)*DCOS(2*zeta)*DCOS(2*Q)
  828. X              +.0000561*DCOS(3*zeta)*DCOS(2*Q);
  829. X  e_pert += .0000343*DCOS(4*zeta)*DCOS(2*Q)
  830. X      +.0000469*DSIN(3*psi)*DCOS(2*Q)
  831. X    -.0000242*DCOS(3*psi)*DCOS(2*Q)
  832. X      -.0000205*DSIN(zeta)*DSIN(3*Q)
  833. X        +.0000262*DSIN(3*zeta)*DSIN(3*Q)
  834. X          +.0000208*DCOS(zeta)*DCOS(3*Q)
  835. X        -.0000271*DCOS(3*zeta)*DCOS(3*Q)
  836. X          -.0000382*DCOS(3*zeta)*DSIN(4*Q)
  837. X            -.0000376*DSIN(3*zeta)*DCOS(4*Q);
  838. X  B = (0.077108 + 0.007186*v - 0.001533*v*v)*DSIN(V)
  839. X    +(0.045803 - 0.014766*v - 0.000536*v*v)*DCOS(V)
  840. X      -0.007075*DSIN(zeta)
  841. X    -0.075825*DSIN(zeta)*DSIN(Q)
  842. X      -0.024839*DSIN(2*zeta)*DSIN(Q)
  843. X        -0.008631*DSIN(3*zeta)*DSIN(Q)
  844. X          -0.072586*DCOS(Q)
  845. X        -0.150383*DCOS(zeta)*DCOS(Q)
  846. X          +0.026897*DCOS(2*zeta)*DCOS(Q)
  847. X            +0.010053*DCOS(3*zeta)*DCOS(Q);
  848. X  B += -(0.013597 +0.001719*v)*DSIN(zeta)*DSIN(2*Q)
  849. X      +(-0.007742 + 0.001517*v)*DCOS(zeta)*DSIN(2*Q)
  850. X    +(0.013586 - 0.001375*v)*DCOS(2*zeta)*DSIN(2*Q)
  851. X      +(-0.013667 + 0.001239*v)*DSIN(zeta)*DCOS(2*Q)
  852. X        +0.011981*DSIN(2*zeta)*DCOS(2*Q)
  853. X          +(0.014861 + 0.001136*v)*DCOS(zeta)*DCOS(2*Q)
  854. X        -(0.013064 + 0.001628*v)*DCOS(2*zeta)*DCOS(2*Q);
  855. X
  856. X  a_pert = .000572*DSIN(V) -.001590*DSIN(2*zeta)*DCOS(Q)
  857. X    +.002933*DCOS(V) -.000647*DSIN(3*zeta)*DCOS(Q)
  858. X      +.033629*DCOS(zeta) -.000344*DSIN(4*zeta)*DCOS(Q)
  859. X    -.003081*DCOS(2*zeta) +.002885*DCOS(zeta)*DCOS(Q)
  860. X      -.001423*DCOS(3*zeta) +(.002172 + .000102*v)*DCOS(2*zeta)*DCOS(Q)
  861. X        -.000671*DCOS(4*zeta) +.000296*DCOS(3*zeta)*DCOS(Q)
  862. X          -.000320*DCOS(5*zeta) -.000267*DSIN(2*zeta)*DSIN(2*Q);
  863. X  a_pert += .001098*DSIN(Q) -.000778*DCOS(zeta)*DSIN(2*Q)
  864. X      -.002812*DSIN(zeta)*DSIN(Q) +.000495*DCOS(2*zeta)*DSIN(2*Q)
  865. X    +.000688*DSIN(2*zeta)*DSIN(Q) +.000250*DCOS(3*zeta)*DSIN(2*Q);
  866. X  a_pert += -.000393*DSIN(3*zeta)*DSIN(Q)
  867. X      -.000228*DSIN(4*zeta)*DSIN(Q)
  868. X    +.002138*DCOS(zeta)*DSIN(Q)
  869. X      -.000999*DCOS(2*zeta)*DSIN(Q)
  870. X        -.000642*DCOS(3*zeta)*DSIN(Q)
  871. X          -.000325*DCOS(4*zeta)*DSIN(Q)
  872. X        -.000890*DCOS(Q)
  873. X          +.002206*DSIN(zeta)*DCOS(Q);
  874. X  a_pert += -.000856*DSIN(zeta)*DCOS(2*Q)
  875. X      +.000441*DSIN(2*zeta)*DCOS(2*Q)
  876. X    +.000296*DCOS(2*zeta)*DCOS(2*Q)
  877. X      +.000211*DCOS(3*zeta)*DCOS(2*Q)
  878. X        -.000427*DSIN(zeta)*DSIN(3*Q)
  879. X          +.000398*DSIN(3*zeta)*DSIN(3*Q)
  880. X        +.000344*DCOS(zeta)*DCOS(3*Q)
  881. X          -.000427*DCOS(3*zeta)*DCOS(3*Q);
  882. X
  883. X  L_ += A;
  884. X  M_ += A - B / e_;
  885. X  e_ += e_pert;
  886. X  a_ += a_pert;
  887. X  omega_ += B;
  888. X
  889. X  /* Calculate E and nu */
  890. X  anom_calc(M_, e_, &E, &nu);
  891. X  r = a_ * (1 - e_ * DCOS(E));
  892. X
  893. X
  894. X  u = L_ + nu - M_ - Omega_;
  895. X  *C_p = nu - M_;
  896. X  l = into_range(RAD_TO_DEG * atan2(DCOS(i_) * DSIN(u), DCOS(u)) + Omega_);
  897. X  b = RAD_TO_DEG * asin(DSIN(u)*DSIN(i_));
  898. X
  899. X  b_pert = 0.000747*DCOS(zeta)*DSIN(Q)
  900. X      +0.001069*DCOS(zeta)*DCOS(Q)
  901. X    +0.002108*DSIN(2*zeta)*DSIN(2*Q)
  902. X      +0.001261*DCOS(2*zeta)*DSIN(2*Q)
  903. X        +0.001236*DSIN(2*zeta)*DCOS(2*Q)
  904. X          -0.002075*DCOS(2*zeta)*DCOS(2*Q);
  905. X
  906. X
  907. X  *r_p = r;
  908. X  *l_p = l;
  909. X  *b_p = b + b_pert;
  910. X}
  911. X
  912. Xstatic void uranus(jd, L_, a_, e_, i_, omega_, Omega_, M_, r_p, l_p, b_p, C_p)
  913. Xdouble jd, L_, a_, e_, i_, omega_, Omega_, M_;
  914. Xdouble *r_p, *l_p, *b_p, *C_p;
  915. X{
  916. X  double E, nu;
  917. X  double M, M1, M2, M4, M5, M6;
  918. X  double r, l, b;
  919. X
  920. X  double u;            /* argument of latitude */
  921. X  double T;
  922. X  double A, B, e_pert, a_pert, b_pert, v, zeta, eta, theta;
  923. X  double P, Q, S, V, W, G, H;
  924. X  double r_pert, l_pert;
  925. X
  926. X  M = M_sun(jd);
  927. X  M1 = M_1(jd);
  928. X  M2 = M_2(jd);
  929. X  M4 = M_4(jd);
  930. X  M5 = M_5(jd);
  931. X  M6 = M_6(jd);
  932. X
  933. X  /* for perturbations */
  934. X  T = (jd - 2415020.0)/36525.0;
  935. X
  936. X  v = T/5.0 + 0.1;
  937. X  P = 237.47555 +3034.9061*T;
  938. X  Q = 265.91650 + 1222.1139*T;
  939. X  S = 243.51721 + 428.4677*T;
  940. X  V = 5.0*Q -2.0*P;
  941. X  W = 2.0*P - 6.0*Q +3.0*S;
  942. X  G = 83.76922 + 218.4901*T;
  943. X  H = 2.0*G - S;
  944. X  zeta = S - P;
  945. X  eta = S - Q;
  946. X  theta = G - S;
  947. X
  948. X  A = (0.864319 - 0.001583*v)*DSIN(H)
  949. X    +(0.082222 - 0.006833*v)*DCOS(H)
  950. X      +0.036017*DSIN(2*H)
  951. X    -0.003019*DCOS(2*H)
  952. X      +0.008122*DSIN(W);
  953. X
  954. X  e_pert = (-.0003349 + .0000163*v)*DSIN(H)
  955. X    +.0020981*DCOS(H)
  956. X      +.0001311*DCOS(H);
  957. X
  958. X  B = 0.120303*DSIN(H)
  959. X    +(0.019472 - 0.000947*v)*DCOS(H)
  960. X      +0.006197*DSIN(2*H);
  961. X
  962. X  a_pert = - 0.003825*DCOS(H);
  963. X
  964. X  L_ += A;
  965. X  M_ += A - B / e_;
  966. X  e_ += e_pert;
  967. X  a_ += a_pert;
  968. X  omega_ += B;
  969. X
  970. X  /* Calculate E and nu */
  971. X  anom_calc(M_, e_, &E, &nu);
  972. X  r = a_ * (1 - e_ * DCOS(E));
  973. X
  974. X
  975. X  u = L_ + nu - M_ - Omega_;
  976. X  *C_p = nu - M_;
  977. X  l = into_range(RAD_TO_DEG * atan2(DCOS(i_) * DSIN(u), DCOS(u)) + Omega_);
  978. X  b = RAD_TO_DEG * asin(DSIN(u)*DSIN(i_));
  979. X
  980. X  l_pert = (0.010122 - 0.000988*v)*DSIN(S+eta)
  981. X      +(-0.038581 + 0.002031*v - 0.001910*v*v)*DCOS(S+eta)
  982. X    +(0.034964 - 0.001038*v + 0.000868*v*v)*DCOS(2*S+eta)
  983. X      +0.005594*DSIN(S +3*theta);
  984. X  l_pert += -0.014808*DSIN(zeta)
  985. X      -0.005794*DSIN(eta)
  986. X    +0.002347*DCOS(eta)
  987. X      +0.009872*DSIN(theta)
  988. X        +0.008803*DSIN(2*theta)
  989. X          -0.004308*DSIN(3*theta);
  990. X  b_pert = (0.000458*DSIN(eta) - 0.000642*DCOS(eta) - 0.000517*DCOS(4*theta))
  991. X      *DSIN(S)
  992. X    -(0.000347*DSIN(eta) + 0.000853*DCOS(eta) + 0.000517*DSIN(4*eta))
  993. X      *DCOS(S)
  994. X        +0.000403*(DCOS(2*theta)*DSIN(2*S) + DSIN(2*theta)*DCOS(2*S));
  995. X  r_pert = -.025948
  996. X      +.004985*DCOS(zeta)
  997. X    -.001230*DCOS(S)
  998. X      +.003354*DCOS(eta)
  999. X        +(.005795*DCOS(S) - .001165*DSIN(S) + .001388*DCOS(2*S))*DSIN(eta)
  1000. X          +(.001351*DCOS(S) + .005702*DSIN(S) + .001388*DSIN(2*S))*DCOS(eta)
  1001. X        +.000904*DCOS(2*theta)
  1002. X          +.000894*(DCOS(theta) - DCOS(3*theta));
  1003. X
  1004. X  *r_p = r + r_pert;
  1005. X  *l_p = l + l_pert;
  1006. X  *b_p = b + b_pert;
  1007. X}
  1008. X
  1009. Xstatic void neptune(jd, L_, a_, e_, i_, omega_, Omega_, M_, r_p, l_p, b_p, C_p)
  1010. Xdouble jd, L_, a_, e_, i_, omega_, Omega_, M_;
  1011. Xdouble *r_p, *l_p, *b_p, *C_p;
  1012. X{
  1013. X  double E, nu;
  1014. X  double M, M1, M2, M4, M5, M6;
  1015. X  double r, l, b;
  1016. X
  1017. X  double u;            /* argument of latitude */
  1018. X  double T;
  1019. X  double A, B, e_pert, a_pert, b_pert, v, zeta, eta, theta;
  1020. X  double P, Q, S, V, W, G, H;
  1021. X  double r_pert, l_pert;
  1022. X
  1023. X  M = M_sun(jd);
  1024. X  M1 = M_1(jd);
  1025. X  M2 = M_2(jd);
  1026. X  M4 = M_4(jd);
  1027. X  M5 = M_5(jd);
  1028. X  M6 = M_6(jd);
  1029. X
  1030. X  /* for perturbations */
  1031. X  T = (jd - 2415020.0)/36525.0;
  1032. X
  1033. X  v = T/5.0 + 0.1;
  1034. X  P = 237.47555 +3034.9061*T;
  1035. X  Q = 265.91650 + 1222.1139*T;
  1036. X  S = 243.51721 + 428.4677*T;
  1037. X  V = 5.0*Q -2.0*P;
  1038. X  W = 2.0*P - 6.0*Q +3.0*S;
  1039. X  G = 83.76922 + 218.4901*T;
  1040. X  H = 2.0*G - S;
  1041. X  zeta = S - P;
  1042. X  eta = S - Q;
  1043. X  theta = G - S;
  1044. X
  1045. X  A = (-0.589833 + 0.001089*v)*DSIN(H)
  1046. X    +(-0.056094 + 0.004658*v)*DCOS(H)
  1047. X      -0.024286*DSIN(2*H);
  1048. X
  1049. X  e_pert = .0004389*DSIN(H)
  1050. X    +.0004262*DCOS(H)
  1051. X      +.0001129*DSIN(2*H)
  1052. X    +.0001089*DCOS(2*H);
  1053. X
  1054. X  B = 0.024039*DSIN(H)
  1055. X    -0.025303*DCOS(H)
  1056. X      +0.006206*DSIN(2*H)
  1057. X    -0.005992*DCOS(2*H);
  1058. X
  1059. X  a_pert = -0.000817*DSIN(H)
  1060. X    +0.008189*DCOS(H)
  1061. X      +0.000781*DCOS(2*H);
  1062. X
  1063. X  L_ += A;
  1064. X  M_ += A - B / e_;
  1065. X  e_ += e_pert;
  1066. X  a_ += a_pert;
  1067. X  omega_ += B;
  1068. X
  1069. X  /* Calculate E and nu */
  1070. X  anom_calc(M_, e_, &E, &nu);
  1071. X  r = a_ * (1 - e_ * DCOS(E));
  1072. X
  1073. X
  1074. X  u = L_ + nu - M_ - Omega_;
  1075. X  *C_p = nu - M_;
  1076. X  l = into_range(RAD_TO_DEG * atan2(DCOS(i_) * DSIN(u), DCOS(u)) + Omega_);
  1077. X  b = RAD_TO_DEG * asin(DSIN(u)*DSIN(i_));
  1078. X
  1079. X  l_pert = -0.009556*DSIN(zeta)
  1080. X      -0.005178*DSIN(eta)
  1081. X    +0.002572*DSIN(2*theta)
  1082. X      -0.002972*DCOS(2*theta)*DSIN(G)
  1083. X        -0.002833*DSIN(2*theta)*DCOS(G);
  1084. X
  1085. X  b_pert = 0.000336*DCOS(2*theta)*DSIN(G)
  1086. X      +0.000364*DSIN(2*theta)*DCOS(G);
  1087. X
  1088. X  r_pert = -.040596
  1089. X      +.004992*DCOS(zeta)
  1090. X    +.002744*DCOS(eta)
  1091. X      +.002044*DCOS(theta)
  1092. X        +.001051*DCOS(2*theta);
  1093. X
  1094. X  *r_p = r + r_pert;
  1095. X  *l_p = l + l_pert;
  1096. X  *b_p = b + b_pert;
  1097. X}
  1098. X
  1099. END_OF_FILE
  1100. if test 28915 -ne `wc -c <'observe/planetcalc.c'`; then
  1101.     echo shar: \"'observe/planetcalc.c'\" unpacked with wrong size!
  1102. fi
  1103. # end of 'observe/planetcalc.c'
  1104. fi
  1105. echo shar: End of archive 15 \(of 32\).
  1106. cp /dev/null ark15isdone
  1107. MISSING=""
  1108. for I in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 ; do
  1109.     if test ! -f ark${I}isdone ; then
  1110.     MISSING="${MISSING} ${I}"
  1111.     fi
  1112. done
  1113. if test "${MISSING}" = "" ; then
  1114.     echo You have unpacked all 32 archives.
  1115.     rm -f ark[1-9]isdone ark[1-9][0-9]isdone
  1116. else
  1117.     echo You still need to unpack the following archives:
  1118.     echo "        " ${MISSING}
  1119. fi
  1120. ##  End of shell archive.
  1121. exit 0
  1122.  
  1123.  
  1124.