Usenet 1994 October/usenetsourcesnewsgroupsinfomagicoctober1994disk2.iso/misc/volume11/starchart/part15

Text File | 1990-03-25 | 30.9 KB | 1,131 lines

Newsgroups: comp.sources.misc
subject: v11i043: starchart 3.2 Part 15/32
from: [e-mail]
Sender: [e-mail] (Brandon S. Allbery - comp.sources.misc)

Posting-number: Volume 11, Issue 43
Submitted-by: [e-mail]
Archive-name: starchart/part15

#! /bin/sh
# This is a shell archive.  Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file".  To overwrite existing
# files, type "sh file -c".  You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g..  If this archive is complete, you
# will see the following message at the end:
#        "End of archive 15 (of 32)."
# Contents:  observe/planetcalc.c
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'observe/planetcalc.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'observe/planetcalc.c'\"
else
echo shar: Extracting \"'observe/planetcalc.c'\" \(28915 characters\)
sed "s/^X//" >'observe/planetcalc.c' <<'END_OF_FILE'
X/*
X * planetcalc.c
X * planetary positions calculations
X *
X * Copyright (c) 1990 by Craig Counterman. All rights reserved.
X *
X * This software may be redistributed freely, not sold.
X * This copyright notice and disclaimer of warranty must remain
X *    unchanged.
X *
X * No representation is made about the suitability of this
X * software for any purpose.  It is provided "as is" without express or
X * implied warranty, to the extent permitted by applicable law.
X *
X */
X
X#ifndef  lint
Xstatic char rcsid[] =
X  "$Header: planetcalc.c,v 1.7 90/02/19 17:21:21 ccount Exp $";
X#endif
X
X
X#include <stdio.h>
X#include <math.h>
X#include "observe.h"
X#include "degree.h"
X
Xstatic double M_sun(), M_1(), M_2(), M_4(), M_5(), M_6();
X
X/* Approximate M for each planet
X   from Meeus, chapter 25 */
X
Xstatic double M_sun(jd)
Xdouble jd;
X{
X  double T;
X
X  T = (jd - 2415020.0)/36525.0;
X
X  return into_range(358.47583 + 35999.04975*T - 0.000150*T*T -0.0000033*T*T*T);
X}
X
Xstatic double M_1(jd)
Xdouble jd;
X{
X  double T;
X
X  T = (jd - 2415020.0)/36525.0;
X
X  return into_range(102.27938 + 149472.51529*T + 0.000007*T*T);
X}
X
Xstatic double M_2(jd)
Xdouble jd;
X{
X  double T;
X
X  T = (jd - 2415020.0)/36525.0;
X
X  return into_range(212.60322 + 58517.80387*T +0.001286*T*T);
X}
X
Xstatic double M_4(jd)
Xdouble jd;
X{
X  double T;
X
X  T = (jd - 2415020.0)/36525.0;
X
X  return into_range(319.51913 + 19139.85475*T + 0.000181*T*T);
X}
X
Xstatic double M_5(jd)
Xdouble jd;
X{
X  double T;
X
X  T = (jd - 2415020.0)/36525.0;
X
X  return into_range(225.32833 + 3034.69202*T - 0.000722*T*T);
X}
X
Xstatic double M_6(jd)
Xdouble jd;
X{
X  double T;
X
X  T = (jd - 2415020.0)/36525.0;
X
X  return into_range(175.46622 +1221.55147*T - 0.000502*T*T);
X}
X
Xstatic char *planet_name[] = {
X  "Mercury",
X  "Venus",
X  "Mars",
X  "Jupiter",
X  "Saturn",
X  "Uranus",
X  "Neptune"
X};
X
Xstatic char *planet_type[] = {
X  "PM",
X  "PV",
X  "Pm",
X  "PJ",
X  "Ps",
X  "PU",
X  "PN"
X};
X
Xstatic char *planet_colr[] = {
X  "r9",
X  "w9",
X  "r9",
X  "w9",
X  "y9",
X  "g9",
X  "c9"
X};
X
Xtypedef struct {
X  double L[4];
X  double a;
X  double e[4];
X  double i[4];
X  double omega[4];
X  double Omega[4];
X  double size_1au;
X  double mag0;
X} pelements;
X
Xpelements peles[] = {
X  {   /* Mercury */
X    {178.179078, 149474.07078, 3.011e-4, 0.0},
X    0.3870986,
X    {.20561421, 2.046e-5, 3e-8,  0.0},
X    {7.002881, 1.8608e-3, -1.83e-5, 0.0},
X    {28.753753, 0.3702806, 0.0001208, 0.0},
X    {47.145944, 1.1852083, 1.739e-4, 0.0},
X    6.74,
X    -0.42
X    },
X  {   /* Venus */
X    {342.767053, 58519.21191, 0.0003097, 0.0},
X    0.7233316,
X    {0.00682069, -0.00004774, 0.000000091, 0.0},
X    {3.393631, 0.0010058, -0.0000010, 0.0},
X    {54.384186, 0.5081861, -0.0013864, 0.0},
X    {75.779647, 0.8998500, 0.0004100, 0.0},
X    16.92,
X    -4.4
X    },
X  {   /* Mars */
X    {293.737334, 19141.69551, 0.0003107, 0.0},
X    1.5236883,
X    {0.09331290, 0.000092064, -0.000000077, 0.0},
X    {1.850333, -0.0006750, 0.0000126, 0.0},
X    {285.431761, 1.0697667, 0.0001313, 0.00000414},
X    {48.786442, 0.7709917, -0.0000014, -0.00000533},
X    9.36,
X    -1.52
X    },
X  {   /* Jupiter */
X    {238.049257, 3036.301986, 0.0003347, -0.00000165},
X    5.202561,
X    {0.04833475, 0.000164180, -0.0000004676, -0.0000000017},
X    {1.308736, -0.0056961, 0.0000039, 0.0},
X    {273.277558, 0.5994317, 0.00070405, 0.00000508},
X    {99.443414, 1.0105300, 0.00035222, -0.00000851},
X    196.74,
X    -9.4
X    },
X  {   /* Saturn */
X    {266.564377, 1223.509884, 0.0003245, -0.0000058},
X    9.554747,
X    {0.05589232, -0.00034550, -0.000000728, 0.00000000074},
X    {2.492519, -0.0039189, -0.00001549, 0.00000004},
X    {338.307800, 1.0852207, 0.00097854, 0.00000992},
X    {112.790414, 0.8731951, -0.00015218, -0.00000531},
X    165.6,
X    -8.88
X    },
X  {   /* Uranus */
X    {244.197470, 429.863546, 0.0003160, -0.00000060},
X    19.21814,
X    {0.0463444, -0.00002658, 0.000000077, 0.0},
X    {0.772464, 0.0006253, 0.0000395, 0.0},
X    {98.071581, 0.9857650, -0.0010745, -0.00000061},
X    {73.477111, 0.4986678, 0.0013117, 0.0},
X    65.8,
X    -7.19
X    },
X  {   /* Neptune */
X    {84.457994, 219.885914, 0.0003205, -0.00000060},
X    30.10957,
X    {0.00899704, 0.000006330, -0.000000002, 0.0},
X    {1.779242, -0.0095436, -0.0000091, 0.0},
X    {276.045975, 0.3256394, 0.00014095, 0.000004113},
X    {130.681389, 1.0989350, 0.00024987, -0.000004718},
X    62.2,
X    -6.87
X    }
X};
X
Xtypedef struct {
X  double alpha_1, delta_1;
X  double W_0, W_dot;
X} rot_els_t;
Xrot_els_t rot_els[] = {
X  {                /* Mercury */
X    280.98, 61.44,
X    142.14, 6.13850
X    },
X  {                /* Venus */
X    272.78, 67.21,
X    353.00, -1.48142
X    },
X  {                /* Mars */
X    317.61, 52.85,
X    237.21, 350.89198
X    },
X  {                /* Jupiter III */
X    268.04, 64.49,
X    156.03, 870.53600
X    },
X  {                /* Saturn III */
X    40.09, 83.49,
X    223.60, 810.79390
X    },
X  {                /* Uranus */
X    257.27, -15.09,
X    287.17, -499.42197
X    },
X  {                /* Neptune */
X    295.24, 40.62,
X    315.33, 468.75000
X    }
X};
X
Xstatic double polynom(jd, a)
X     double jd;
X     double a[4];
X{
X  double T;
X
X  T = (jd - 2415020.0)/36525.0;
X
X  return (a[0] + a[1]*T + a[2]*T*T + a[3]*T*T*T);
X}
X
Xstatic void mercury(), venus(), mars(),
X  jupiter(), saturn(), uranus(), neptune();
X
X
X/* Calculate alpha and delta
X   from lambda and beta (and epsilon)
X   which are from r, Delta, psi, b, l, and Theta (in sun_data)
X   which are from u, i (given), and Omega (given)
X   u is from L (given), nu, and M
X   nu and M are calculated.
X   r is from E, a (given) and e (given)
X   E is calculated
X
X   calculate mag from Delta, size form Delta, phase (== beta).
X */
X
Xvoid planet_pos(jd, sun_data, nplanet, data)
X     double jd;            /* time, jd */
X     sun_data_t sun_data;
X     int nplanet;        /* Planet number, 0 = mercury */
X     planet_data_t *data;
X{
X  double L_, a_, e_, i_, omega_, Omega_, M_;
X  double r;            /* radius distance to sun */
X  double l, b;            /* ecliptical longitude and latitude */
X  double Delta;            /* Distance to earth */
X  double lambda, beta;        /* geocentric longitude and latitude */
X  double alpha, delta;        /* R.A. and dec. both degrees */
X  double alpha2000, delta2000;    /* R.A. and dec. both degrees equin 2000.0 */
X  double psi;            /* elongation */
X  double N, D;            /* temporary variables */
X  double Theta;            /* Theta of the sun */
X  double epsilon;        /* obliquity */
X  double Cen;            /* center */
X  double A, B;            /* used in calculating p_n */
X  double W_1, K;        /* Used in calculating lambda_e */
X
X  L_ = into_range(polynom(jd, peles[nplanet].L));
X  a_ = peles[nplanet].a;
X  e_ = polynom(jd, peles[nplanet].e);
X  i_ = polynom(jd, peles[nplanet].i);
X  omega_ = into_range(polynom(jd, peles[nplanet].omega));
X  Omega_ = into_range(polynom(jd, peles[nplanet].Omega));
X  M_ = into_range(L_ - omega_ - Omega_);
X
X  /* Perturb */
X  switch (nplanet) {
X  case 0:
X    mercury(jd, L_, a_, e_, i_, omega_, Omega_, M_, &r, &l, &b, &Cen);
X    break;
X  case 1:
X    venus(jd, L_, a_, e_, i_, omega_, Omega_, M_, &r, &l, &b, &Cen);
X    break;
X  case 2:
X    mars(jd, L_, a_, e_, i_, omega_, Omega_, M_, &r, &l, &b, &Cen);
X    break;
X  case 3:
X    jupiter(jd, L_, a_, e_, i_, omega_, Omega_, M_, &r, &l, &b, &Cen);
X    break;
X  case 4:
X    saturn(jd, L_, a_, e_, i_, omega_, Omega_, M_, &r, &l, &b, &Cen);
X    break;
X  case 5:
X    uranus(jd, L_, a_, e_, i_, omega_, Omega_, M_, &r, &l, &b, &Cen);
X    break;
X  case 6:
X    neptune(jd, L_, a_, e_, i_, omega_, Omega_, M_, &r, &l, &b, &Cen);
X    break;
X  default:
X    break;
X  };
X
X  Theta = sun_data.Theta;
X  N = r * DCOS(b) * DSIN(l - Theta);
X  D = r * DCOS(b) * DCOS(l - Theta) + sun_data.R;
X
X  epsilon = obl_jd(jd);
X
X  lambda = into_range(RAD_TO_DEG * atan2(N, D)) + Theta;
X  Delta = sqrt(N*N + D*D + (r * DSIN(b))*(r * DSIN(b)));
X  beta = RAD_TO_DEG * asin(r * DSIN(b) / Delta);
X  psi = RAD_TO_DEG * acos(DCOS(beta) * DCOS(lambda - Theta));
X  if (into_range(lambda - Theta) > 180.0)
X    psi = -psi;
X  alpha = RAD_TO_DEG * atan2(DSIN(lambda)*DCOS(epsilon)
X                 - DTAN(beta) * DSIN(epsilon),
X                 DCOS(lambda));
X  delta = RAD_TO_DEG * asin(DSIN(beta)*DCOS(epsilon)
X                + DCOS(beta)*DSIN(epsilon)*DSIN(lambda));
X  alpha = into_range(alpha);
X
X/* should correct for nutation and aberration */
X  data->name = planet_name[nplanet];
X  data->type = planet_type[nplanet];
X  data->color = planet_colr[nplanet];
X  data->alpha = alpha;
X  data->delta = delta;
X  precess(2000.0 - (2451545.0 - jd) / 365.25,
X      2000.0, alpha, delta, &alpha2000, &delta2000);
X  data->alpha2000 = alpha2000;
X  data->delta2000 = delta2000;
X  data->l = l;
X  data->b = b;
X  data->lambda = lambda;
X  data->beta = beta;
X  data->psi = psi;
X  data->phase =
X    DACOS((r*r + Delta*Delta - sun_data.R*sun_data.R) / (2*r*Delta));
X  if (psi < 0) data->phase = -data->phase;
X  data->r = r;
X  data->Delta = Delta;
X  data->illum_frac = ((r+Delta)*(r+Delta) - sun_data.R*sun_data.R)/(4*r*Delta);
X  data->chi = /* position angle of bright limb */
X    DATAN2(DCOS(sun_data.delta) * DSIN(sun_data.alpha - alpha),
X      DCOS(delta) * DSIN(sun_data.delta)
X      - DSIN(delta) * DCOS(sun_data.delta) * DCOS(sun_data.alpha - alpha));
X  data->Cen = Cen;
X
X  data->mag = 5.0 * log10(r*Delta)
X    - 2.5 * log10(data->illum_frac)
X    + peles[nplanet].mag0;
X
X  data->size = peles[nplanet].size_1au / Delta;
X
X  data->rotation_elements.beta_e =
X    DASIN(-DSIN(rot_els[nplanet].delta_1)*DSIN(delta)
X      - DCOS(rot_els[nplanet].delta_1) * DCOS(delta)
X          * DCOS(rot_els[nplanet].alpha_1 - alpha));
X
X  A = DCOS(rot_els[nplanet].delta_1) * DSIN(rot_els[nplanet].alpha_1 - alpha);
X  A /= DCOS(data->rotation_elements.beta_e);
X  B = DSIN(rot_els[nplanet].delta_1) * DCOS(delta)
X    - DCOS(rot_els[nplanet].delta_1) * DSIN(delta)
X      * DCOS(rot_els[nplanet].alpha_1 - alpha);
X  B /= DCOS(data->rotation_elements.beta_e);
X  data->rotation_elements.p_n = DATAN2(A, B);
X
X  W_1 = into_range(rot_els[nplanet].W_0
X           + rot_els[nplanet].W_dot
X             * (jd - 2447526.5 - 0.0057755 * Delta));
X
X  A = -DCOS(rot_els[nplanet].delta_1) * DSIN(delta)
X    + DSIN(rot_els[nplanet].delta_1) * DCOS(delta)
X      * DCOS(rot_els[nplanet].alpha_1 - alpha);
X  A /= DCOS(data->rotation_elements.beta_e);
X  B = DCOS(delta) * DSIN(rot_els[nplanet].alpha_1 - alpha);
X  B /= DCOS(data->rotation_elements.beta_e);
X  K = DATAN2(A, B);
X  if (rot_els[nplanet].W_dot < 0.0)
X    data->rotation_elements.lambda_e = into_range(K - W_1);
X  else
X    data->rotation_elements.lambda_e = into_range(W_1 - K);
X}
X
X
Xstatic void mercury(jd, L_, a_, e_, i_, omega_, Omega_, M_, r_p, l_p, b_p, C_p)
Xdouble jd, L_, a_, e_, i_, omega_, Omega_, M_;
Xdouble *r_p, *l_p, *b_p, *C_p;
X{
X  double E, nu;
X  double M1, M2, M4, M5, M6;
X  double r, l, b;
X
X  double u;            /* argument of latitude */
X  double r_pert, l_pert;
X
X  M1 = M_1(jd);
X  M2 = M_2(jd);
X  M4 = M_4(jd);
X  M5 = M_5(jd);
X  M6 = M_6(jd);
X
X
X  /* Calculate E and nu */
X  anom_calc(M_, e_, &E, &nu);
X  r = a_ * (1 - e_ * DCOS(E));
X
X  u = L_ + nu - M_ - Omega_;
X  *C_p = nu - M_;
X  l = into_range(RAD_TO_DEG * atan2(DCOS(i_) * DSIN(u), DCOS(u)) + Omega_);
X  b = RAD_TO_DEG * asin(DSIN(u)*DSIN(i_));
X
X  /* Perturbations */
X  l_pert = 0.00204*DCOS(5*M2-2*M1+12.220)
X    +0.00103*DCOS(2*M2-M1-160.692)
X      +0.00091*DCOS(2*M5-M1-37.003)
X    +0.00078*DCOS(5*M2-3*M1+10.137);
X  r_pert = 0.000007525*DCOS(2*M5-M1+53.013)
X    +0.000006802*DCOS(5*M2-3*M1-259.918)
X      +0.000005457*DCOS(2*M2-2*M1-71.188)
X    +0.000003569*DCOS(5*M2-M1-77.75);
X
X  *r_p = r + r_pert;
X  *l_p = l + l_pert;
X  *b_p = b;
X}
X
Xstatic void venus(jd, L_, a_, e_, i_, omega_, Omega_, M_, r_p, l_p, b_p, C_p)
Xdouble jd, L_, a_, e_, i_, omega_, Omega_, M_;
Xdouble *r_p, *l_p, *b_p, *C_p;
X{
X  double E, nu;
X  double M, M1, M2, M4, M5, M6;
X  double r, l, b;
X
X  double u;            /* argument of latitude */
X  double T;
X  double r_pert, l_pert;
X
X  M = M_sun(jd);
X  M1 = M_1(jd);
X  M2 = M_2(jd);
X  M4 = M_4(jd);
X  M5 = M_5(jd);
X  M6 = M_6(jd);
X
X
X  /* Long term perturbation */
X  T = (jd - 2415020.0)/36525.0;
X  l_pert = 0.00077 * DSIN(237.24 + 150.27*T);
X  M_ += l_pert;
X  L_ += l_pert;
X
X  /* Calculate E and nu */
X  anom_calc(M_, e_, &E, &nu);
X  r = a_ * (1 - e_ * DCOS(E));
X
X
X  u = L_ + nu - M_ - Omega_;
X  *C_p = nu - M_;
X  l = into_range(RAD_TO_DEG * atan2(DCOS(i_) * DSIN(u), DCOS(u)) + Omega_);
X  b = RAD_TO_DEG * asin(DSIN(u)*DSIN(i_));
X
X  /* Perturbations */
X  l_pert = 0.00313*DCOS(2*M-2*M2 -148.225)
X    +0.00198*DCOS(3*M-3*M2 +2.565)
X      +0.00136*DCOS(M-M2-119.107)
X    +0.00096*DCOS(3*M-2*M2-135.912)
X      +0.00082*DCOS(M5-M2-208.087);
X  r_pert = 0.000022501 * DCOS(2*M-2*M2-58.208)
X    +0.000019045 * DCOS(3*M-3*M2+92.577)
X      +0.000006887 * DCOS(M5-M2-118.090)
X    +0.000005172 * DCOS(M-M2-29.110)
X      +0.000003620 * DCOS(5*M-4*M2-104.208)
X        +0.000003283 * DCOS(4*M-4*M2+63.513)
X          +0.000003074 * DCOS(2*M5-2*M2-55.167);
X
X  *r_p = r + r_pert;
X  *l_p = l + l_pert;
X  *b_p = b;
X}
X
Xstatic void mars(jd, L_, a_, e_, i_, omega_, Omega_, M_, r_p, l_p, b_p, C_p)
Xdouble jd, L_, a_, e_, i_, omega_, Omega_, M_;
Xdouble *r_p, *l_p, *b_p, *C_p;
X{
X  double E, nu;
X  double M, M1, M2, M4, M5, M6;
X  double r, l, b;
X
X  double u;            /* argument of latitude */
X  double T;
X  double r_pert, l_pert;
X
X  M = M_sun(jd);
X  M1 = M_1(jd);
X  M2 = M_2(jd);
X  M4 = M_4(jd);
X  M5 = M_5(jd);
X  M6 = M_6(jd);
X
X
X  /* Long term perturbation */
X  T = (jd - 2415020.0)/36525.0;
X  l_pert = -0.01133*DSIN(3*M5-8*M4 +4*M) - 0.00933*DCOS(3*M5-8*M4 +4*M);
X
X  M_ += l_pert;
X  L_ += l_pert;
X
X  /* Calculate E and nu */
X  anom_calc(M_, e_, &E, &nu);
X  r = a_ * (1 - e_ * DCOS(E));
X
X
X  u = L_ + nu - M_ - Omega_;
X  *C_p = nu - M_;
X  l = into_range(RAD_TO_DEG * atan2(DCOS(i_) * DSIN(u), DCOS(u)) + Omega_);
X  b = RAD_TO_DEG * asin(DSIN(u)*DSIN(i_));
X
X  /* Perturbations */
X  l_pert = 0.00705*DCOS(M5-M4-48.958)
X    +0.00607*DCOS(2*M5-M4-188.350)
X      +0.00445*DCOS(2*M5-2*M4-191.897)
X    +0.00388*DCOS(M-2*M4+20.495)
X      +0.00238*DCOS(M-M4+35.097)
X        +0.00204*DCOS(2*M-3*M4+158.638)
X          +0.00177*DCOS(3*M4-M2-57.602)
X        +0.00136*DCOS(2*M-4*M4+154.093)
X          +0.00104*DCOS(M5+17.618);
X  r_pert = 0.000053227*DCOS(M5-M4+41.1306)
X    +0.000050989*DCOS(2*M5-2*M4-101.9847)
X      +0.000038278*DCOS(2*M5-M4-98.3292)
X    +0.000015996*DCOS(M-M4-55.555)
X      +0.000014764*DCOS(2*M-3*M4+68.622)
X        +0.000008966*DCOS(M5-2*M4+43.615);
X  r_pert += 0.000007914*DCOS(3*M5-2*M4-139.737)
X    +0.000007004*DCOS(2*M5-3*M4-102.888)
X      +0.000006620*DCOS(M-2*M4+113.202)
X    +0.000004930*DCOS(3*M5-3*M4-76.243)
X      +0.000004693*DCOS(3*M-5*M4+190.603)
X        +0.000004571*DCOS(2*M-4*M4+244.702)
X          +0.000004409*DCOS(3*M5-M4-115.828);
X
X  *r_p = r + r_pert;
X  *l_p = l + l_pert;
X  *b_p = b;
X}
X
X
Xstatic void jupiter(jd, L_, a_, e_, i_, omega_, Omega_, M_, r_p, l_p, b_p, C_p)
Xdouble jd, L_, a_, e_, i_, omega_, Omega_, M_;
Xdouble *r_p, *l_p, *b_p, *C_p;
X{
X  double E, nu;
X  double M, M1, M2, M4, M5, M6;
X  double r, l, b;
X
X  double u;            /* argument of latitude */
X  double T;
X  double A, B, e_pert, a_pert, v, zeta;
X  double P, Q, S, V, W;
X
X  M = M_sun(jd);
X  M1 = M_1(jd);
X  M2 = M_2(jd);
X  M4 = M_4(jd);
X  M5 = M_5(jd);
X  M6 = M_6(jd);
X
X  /* for perturbations */
X  T = (jd - 2415020.0)/36525.0;
X
X  v = T/5.0 + 0.1;
X  P = 237.47555 +3034.9061*T;
X  Q = 265.91650 + 1222.1139*T;
X  S = 243.51721 + 428.4677*T;
X  V = 5.0*Q -2.0*P;
X  W = 2.0*P - 6.0*Q +3.0*S;
X  zeta = Q - P;
X
X
X  A = (0.331364 - 0.010281*v - 0.004692*v*v)*DSIN(V)
X    +(0.003228 - 0.064436*v + 0.002075*v*v)*DCOS(V)
X      -(0.003083 + 0.000275*v - 0.000489*v*v)*DSIN(2*V)
X    +0.002472*DSIN(W)
X      +0.013619*DSIN(zeta)
X        +0.018472*DSIN(2*zeta)
X          +0.006717*DSIN(3*zeta)
X        +0.002775*DSIN(4*zeta)
X          +(0.007275 - 0.001253*v)*DSIN(zeta)*DSIN(Q)
X            +0.006417*DSIN(2*zeta)*DSIN(Q)
X              +0.002439*DSIN(3*zeta)*DSIN(Q);
X  A += -(0.033839 + 0.001125*v)*DCOS(zeta)*DSIN(Q)
X    -0.003767*DCOS(2*zeta)*DSIN(Q)
X      -(0.035681 + 0.001208*v)*DSIN(zeta)*DCOS(Q)
X    -0.004261*DSIN(2*zeta)*DCOS(Q)
X      +0.002178*DCOS(Q)
X        +(-0.006333 + 0.001161*v)*DCOS(zeta)*DCOS(Q)
X          -0.006675*DCOS(2*zeta)*DCOS(Q)
X        -0.002664*DCOS(3*zeta)*DCOS(Q)
X          -0.002572*DSIN(zeta)*DSIN(2*Q)
X            -0.003567*DSIN(2*zeta)*DSIN(2*Q)
X              +0.002094*DCOS(zeta)*DCOS(2*Q)
X            +0.003342*DCOS(2*zeta)*DCOS(2*Q);
X
X  e_pert = (.0003606 + .0000130*v - .0000043*v*v)*DSIN(V)
X    +(.0001289 - .0000580*v)*DCOS(V)
X      -.0006764*DSIN(zeta)*DSIN(Q)
X    -.0001110*DSIN(2*zeta)*DSIN(Q)
X      -.0000224*DSIN(3*zeta)*DSIN(Q)
X        -.0000204*DSIN(Q)
X          +(.0001284 + .0000116*v)*DCOS(zeta)*DSIN(Q)
X        +.0000188*DCOS(2*zeta)*DSIN(Q)
X          +(.0001460 + .0000130*v)*DSIN(zeta)*DCOS(Q)
X            +.0000224*DSIN(2*zeta)*DCOS(Q)
X              -.0000817*DCOS(Q);
X
X  e_pert += .0006074*DCOS(zeta)*DCOS(Q)
X    +.0000992*DCOS(2*zeta)*DCOS(Q)
X      +.0000508*DCOS(3*zeta)*DCOS(Q)
X    +.0000230*DCOS(4*zeta)*DCOS(Q)
X      +.0000108*DCOS(5*zeta)*DCOS(Q);
X
X  e_pert += -(.0000956 + .0000073*v)*DSIN(zeta)*DSIN(2*Q)
X    +.0000448*DSIN(2*zeta)*DSIN(2*Q)
X      +.0000137*DSIN(3*zeta)*DSIN(2*Q)
X    +(-.0000997 + .0000108*v)*DCOS(zeta)*DSIN(2*Q)
X      +.0000480*DCOS(2*zeta)*DSIN(2*Q);
X
X  e_pert += .0000148*DCOS(3*zeta)*DSIN(2*Q)
X    +(-.0000956 +.0000099*v)*DSIN(zeta)*DCOS(2*Q)
X      +.0000490*DSIN(2*zeta)*DCOS(2*Q)
X    +.0000158*DSIN(3*zeta)*DCOS(2*Q);
X
X  e_pert += .0000179*DCOS(2*Q)
X    +(.0001024 + .0000075*v)*DCOS(zeta)*DCOS(2*Q)
X      -.0000437*DCOS(2*zeta)*DCOS(2*Q)
X    -.0000132*DCOS(3*zeta)*DCOS(2*Q);
X
X  B = (0.007192 - 0.003147*v)*DSIN(V)
X    +(-0.020428 - 0.000675*v + 0.000197*v*v)*DCOS(V)
X      +(0.007269 + 0.000672*v)*DSIN(zeta)*DSIN(Q)
X    -0.004344*DSIN(Q)
X      +0.034036*DCOS(zeta)*DSIN(Q)
X        +0.005614*DCOS(2*zeta)*DSIN(Q)
X          +0.002964*DCOS(3*zeta)*DSIN(Q)
X        +0.037761*DSIN(zeta)*DCOS(Q);
X
X  B += 0.006158*DSIN(2*zeta)*DCOS(Q)
X      -0.006603*DCOS(zeta)*DCOS(Q)
X    -0.005356*DSIN(zeta)*DSIN(2*Q)
X      +0.002722*DSIN(2*zeta)*DSIN(2*Q)
X        +0.004483*DCOS(zeta)*DSIN(2*Q);
X
X  B += -0.002642*DCOS(2*zeta)*DSIN(2*Q)
X    +0.004403*DSIN(zeta)*DCOS(2*Q)
X      -0.002536*DSIN(2*zeta)*DCOS(2*Q)
X    +0.005547*DCOS(zeta)*DCOS(2*Q)
X      -0.002689*DCOS(2*zeta)*DCOS(2*Q);
X
X  a_pert = -.000263*DCOS(V)
X    +.000205*DCOS(zeta)
X      +.000693*DCOS(2*zeta)
X    +.000312*DCOS(3*zeta)
X      +.000147*DCOS(4*zeta)
X        +.000299*DSIN(zeta)*DSIN(Q)
X          +.000181*DCOS(2*zeta)*DSIN(Q)
X        +.000204*DSIN(2*zeta)*DCOS(Q)
X          +.000111*DSIN(3*zeta)*DCOS(Q)
X            -.000337*DCOS(zeta)*DCOS(Q)
X              -.000111*DCOS(2*zeta)*DCOS(Q);
X
X  L_ += A;
X  M_ += A - B / e_;
X  e_ += e_pert;
X  a_ += a_pert;
X  omega_ += B;
X
X  /* Calculate E and nu */
X  anom_calc(M_, e_, &E, &nu);
X  r = a_ * (1 - e_ * DCOS(E));
X
X
X  u = L_ + nu - M_ - Omega_;
X  *C_p = nu - M_;
X  l = into_range(RAD_TO_DEG * atan2(DCOS(i_) * DSIN(u), DCOS(u)) + Omega_);
X  b = RAD_TO_DEG * asin(DSIN(u)*DSIN(i_));
X
X  *r_p = r;
X  *l_p = l;
X  *b_p = b;
X}
X
X
Xstatic void saturn(jd, L_, a_, e_, i_, omega_, Omega_, M_, r_p, l_p, b_p, C_p)
Xdouble jd, L_, a_, e_, i_, omega_, Omega_, M_;
Xdouble *r_p, *l_p, *b_p, *C_p;
X{
X  double E, nu;
X  double M, M1, M2, M4, M5, M6;
X  double r, l, b;
X
X  double u;            /* argument of latitude */
X  double T;
X  double A, B, e_pert, a_pert, b_pert, v, zeta, psi;
X  double P, Q, S, V, W;
X
X  M = M_sun(jd);
X  M1 = M_1(jd);
X  M2 = M_2(jd);
X  M4 = M_4(jd);
X  M5 = M_5(jd);
X  M6 = M_6(jd);
X
X  /* for perturbations */
X  T = (jd - 2415020.0)/36525.0;
X
X  v = T/5.0 + 0.1;
X  P = 237.47555 +3034.9061*T;
X  Q = 265.91650 + 1222.1139*T;
X  S = 243.51721 + 428.4677*T;
X  V = 5.0*Q -2.0*P;
X  W = 2.0*P - 6.0*Q +3.0*S;
X  zeta = Q - P;
X  psi = S - Q;
X
X  A = (-0.814181 + 0.018150*v + 0.016714*v*v)*DSIN(V)
X    +(-0.010497 + 0.160906*v - 0.004100*v*v)*DCOS(V)
X      +0.007581*DSIN(2*V)
X    -0.007986*DSIN(W)
X      -0.148811*DSIN(zeta)
X        -0.040786*DSIN(2*zeta)
X          -0.015208*DSIN(3*zeta)
X        -0.006339*DSIN(4*zeta)
X          -0.006244*DSIN(Q);
X  A += (0.008931 + 0.002728*v)*DSIN(zeta)*DSIN(Q)
X      -0.016500*DSIN(2*zeta)*DSIN(Q)
X    -0.005775*DSIN(3*zeta)*DSIN(Q)
X      +(0.081344 + 0.003206*v)*DCOS(zeta)*DSIN(Q)
X        +0.015019*DCOS(2*zeta)*DSIN(Q)
X          +(0.085581 + 0.002494*v)*DSIN(zeta)*DCOS(Q)
X        +(0.025328 - 0.003117*v)*DCOS(zeta)*DCOS(Q);
X  A += 0.014394*DCOS(2*zeta)*DCOS(Q)
X      +0.006319*DCOS(3*zeta)*DCOS(Q)
X    +0.006369*DSIN(zeta)*DSIN(2*Q)
X      +0.009156*DSIN(2*zeta)*DSIN(2*Q)
X        +0.007525*DSIN(3*psi)*DSIN(2*Q)
X          -0.005236*DCOS(zeta)*DCOS(2*Q)
X        -0.007736*DCOS(2*zeta)*DCOS(2*Q)
X          -0.007528*DCOS(3*psi)*DCOS(2*Q);
X
X  e_pert = (-.0007927 + .0002548*v +.0000091*v*v)*DSIN(V)
X    +(.0013381 + .0001226*v -.0000253*v*v)*DCOS(V)
X      +(.0000248 - .0000121*v)*DSIN(2*V)
X    -(.0000305 + .0000091*v)*DCOS(2*V)
X      +.0000412*DSIN(2*zeta)
X        +.0012415*DSIN(Q)
X          +(.0000390 -.0000617*v)*DSIN(zeta)*DSIN(Q)
X        +(.0000165 - .0000204*v)*DSIN(2*zeta)*DSIN(Q)
X          +.0026599*DCOS(zeta)*DSIN(Q)
X            -.0004687*DCOS(2*zeta)*DSIN(Q);
X  e_pert += -.0001870*DCOS(3*zeta)*DSIN(Q)
X      -.0000821*DCOS(4*zeta)*DSIN(Q)
X    -.0000377*DCOS(5*zeta)*DSIN(Q)
X      +.0000497*DCOS(2*psi)*DSIN(Q)
X        +(.0000163 - .0000611*v)*DCOS(Q)
X          -.0012696*DSIN(zeta)*DCOS(Q)
X        -.0004200*DSIN(2*zeta)*DCOS(Q)
X          -.0001503*DSIN(3*zeta)*DCOS(Q)
X            -.0000619*DSIN(4*zeta)*DCOS(Q)
X              -.0000268*DSIN(5*zeta)*DCOS(Q);
X  e_pert += -(.0000282 + .0001306*v)*DCOS(zeta)*DCOS(Q)
X      +(-.0000086 + .0000230*v)*DCOS(2*zeta)*DCOS(Q)
X    +.0000461*DSIN(2*psi)*DCOS(Q)
X      -.0000350*DSIN(2*Q)
X        +(.0002211 - .0000286*v)*DSIN(zeta)*DSIN(2*Q)
X          -.0002208*DSIN(2*zeta)*DSIN(2*Q)
X        -.0000568*DSIN(3*zeta)*DSIN(2*Q)
X          -.0000346*DSIN(4*zeta)*DSIN(2*Q)
X            -(.0002780 + .0000222*v)*DCOS(zeta)*DSIN(2*Q)
X              +(.0002022 + .0000263*v)*DCOS(2*zeta)*DSIN(2*Q);
X  e_pert += .0000248*DCOS(3*zeta)*DSIN(2*Q)
X      +.0000242*DSIN(3*psi)*DSIN(2*Q)
X    +.0000467*DCOS(3*psi)*DSIN(2*Q)
X      -.0000490*DCOS(2*Q)
X        -(.0002842 + .0000279*v)*DSIN(zeta)*DCOS(2*Q)
X          +(.0000128 + .0000226*v)*DSIN(2*zeta)*DCOS(2*Q)
X        +.0000224*DSIN(3*zeta)*DCOS(2*Q)
X          +(-.0001594 + .0000282*v)*DCOS(zeta)*DCOS(2*Q)
X            +(.0002162 - .0000207*v)*DCOS(2*zeta)*DCOS(2*Q)
X              +.0000561*DCOS(3*zeta)*DCOS(2*Q);
X  e_pert += .0000343*DCOS(4*zeta)*DCOS(2*Q)
X      +.0000469*DSIN(3*psi)*DCOS(2*Q)
X    -.0000242*DCOS(3*psi)*DCOS(2*Q)
X      -.0000205*DSIN(zeta)*DSIN(3*Q)
X        +.0000262*DSIN(3*zeta)*DSIN(3*Q)
X          +.0000208*DCOS(zeta)*DCOS(3*Q)
X        -.0000271*DCOS(3*zeta)*DCOS(3*Q)
X          -.0000382*DCOS(3*zeta)*DSIN(4*Q)
X            -.0000376*DSIN(3*zeta)*DCOS(4*Q);
X  B = (0.077108 + 0.007186*v - 0.001533*v*v)*DSIN(V)
X    +(0.045803 - 0.014766*v - 0.000536*v*v)*DCOS(V)
X      -0.007075*DSIN(zeta)
X    -0.075825*DSIN(zeta)*DSIN(Q)
X      -0.024839*DSIN(2*zeta)*DSIN(Q)
X        -0.008631*DSIN(3*zeta)*DSIN(Q)
X          -0.072586*DCOS(Q)
X        -0.150383*DCOS(zeta)*DCOS(Q)
X          +0.026897*DCOS(2*zeta)*DCOS(Q)
X            +0.010053*DCOS(3*zeta)*DCOS(Q);
X  B += -(0.013597 +0.001719*v)*DSIN(zeta)*DSIN(2*Q)
X      +(-0.007742 + 0.001517*v)*DCOS(zeta)*DSIN(2*Q)
X    +(0.013586 - 0.001375*v)*DCOS(2*zeta)*DSIN(2*Q)
X      +(-0.013667 + 0.001239*v)*DSIN(zeta)*DCOS(2*Q)
X        +0.011981*DSIN(2*zeta)*DCOS(2*Q)
X          +(0.014861 + 0.001136*v)*DCOS(zeta)*DCOS(2*Q)
X        -(0.013064 + 0.001628*v)*DCOS(2*zeta)*DCOS(2*Q);
X
X  a_pert = .000572*DSIN(V) -.001590*DSIN(2*zeta)*DCOS(Q)
X    +.002933*DCOS(V) -.000647*DSIN(3*zeta)*DCOS(Q)
X      +.033629*DCOS(zeta) -.000344*DSIN(4*zeta)*DCOS(Q)
X    -.003081*DCOS(2*zeta) +.002885*DCOS(zeta)*DCOS(Q)
X      -.001423*DCOS(3*zeta) +(.002172 + .000102*v)*DCOS(2*zeta)*DCOS(Q)
X        -.000671*DCOS(4*zeta) +.000296*DCOS(3*zeta)*DCOS(Q)
X          -.000320*DCOS(5*zeta) -.000267*DSIN(2*zeta)*DSIN(2*Q);
X  a_pert += .001098*DSIN(Q) -.000778*DCOS(zeta)*DSIN(2*Q)
X      -.002812*DSIN(zeta)*DSIN(Q) +.000495*DCOS(2*zeta)*DSIN(2*Q)
X    +.000688*DSIN(2*zeta)*DSIN(Q) +.000250*DCOS(3*zeta)*DSIN(2*Q);
X  a_pert += -.000393*DSIN(3*zeta)*DSIN(Q)
X      -.000228*DSIN(4*zeta)*DSIN(Q)
X    +.002138*DCOS(zeta)*DSIN(Q)
X      -.000999*DCOS(2*zeta)*DSIN(Q)
X        -.000642*DCOS(3*zeta)*DSIN(Q)
X          -.000325*DCOS(4*zeta)*DSIN(Q)
X        -.000890*DCOS(Q)
X          +.002206*DSIN(zeta)*DCOS(Q);
X  a_pert += -.000856*DSIN(zeta)*DCOS(2*Q)
X      +.000441*DSIN(2*zeta)*DCOS(2*Q)
X    +.000296*DCOS(2*zeta)*DCOS(2*Q)
X      +.000211*DCOS(3*zeta)*DCOS(2*Q)
X        -.000427*DSIN(zeta)*DSIN(3*Q)
X          +.000398*DSIN(3*zeta)*DSIN(3*Q)
X        +.000344*DCOS(zeta)*DCOS(3*Q)
X          -.000427*DCOS(3*zeta)*DCOS(3*Q);
X
X  L_ += A;
X  M_ += A - B / e_;
X  e_ += e_pert;
X  a_ += a_pert;
X  omega_ += B;
X
X  /* Calculate E and nu */
X  anom_calc(M_, e_, &E, &nu);
X  r = a_ * (1 - e_ * DCOS(E));
X
X
X  u = L_ + nu - M_ - Omega_;
X  *C_p = nu - M_;
X  l = into_range(RAD_TO_DEG * atan2(DCOS(i_) * DSIN(u), DCOS(u)) + Omega_);
X  b = RAD_TO_DEG * asin(DSIN(u)*DSIN(i_));
X
X  b_pert = 0.000747*DCOS(zeta)*DSIN(Q)
X      +0.001069*DCOS(zeta)*DCOS(Q)
X    +0.002108*DSIN(2*zeta)*DSIN(2*Q)
X      +0.001261*DCOS(2*zeta)*DSIN(2*Q)
X        +0.001236*DSIN(2*zeta)*DCOS(2*Q)
X          -0.002075*DCOS(2*zeta)*DCOS(2*Q);
X
X
X  *r_p = r;
X  *l_p = l;
X  *b_p = b + b_pert;
X}
X
Xstatic void uranus(jd, L_, a_, e_, i_, omega_, Omega_, M_, r_p, l_p, b_p, C_p)
Xdouble jd, L_, a_, e_, i_, omega_, Omega_, M_;
Xdouble *r_p, *l_p, *b_p, *C_p;
X{
X  double E, nu;
X  double M, M1, M2, M4, M5, M6;
X  double r, l, b;
X
X  double u;            /* argument of latitude */
X  double T;
X  double A, B, e_pert, a_pert, b_pert, v, zeta, eta, theta;
X  double P, Q, S, V, W, G, H;
X  double r_pert, l_pert;
X
X  M = M_sun(jd);
X  M1 = M_1(jd);
X  M2 = M_2(jd);
X  M4 = M_4(jd);
X  M5 = M_5(jd);
X  M6 = M_6(jd);
X
X  /* for perturbations */
X  T = (jd - 2415020.0)/36525.0;
X
X  v = T/5.0 + 0.1;
X  P = 237.47555 +3034.9061*T;
X  Q = 265.91650 + 1222.1139*T;
X  S = 243.51721 + 428.4677*T;
X  V = 5.0*Q -2.0*P;
X  W = 2.0*P - 6.0*Q +3.0*S;
X  G = 83.76922 + 218.4901*T;
X  H = 2.0*G - S;
X  zeta = S - P;
X  eta = S - Q;
X  theta = G - S;
X
X  A = (0.864319 - 0.001583*v)*DSIN(H)
X    +(0.082222 - 0.006833*v)*DCOS(H)
X      +0.036017*DSIN(2*H)
X    -0.003019*DCOS(2*H)
X      +0.008122*DSIN(W);
X
X  e_pert = (-.0003349 + .0000163*v)*DSIN(H)
X    +.0020981*DCOS(H)
X      +.0001311*DCOS(H);
X
X  B = 0.120303*DSIN(H)
X    +(0.019472 - 0.000947*v)*DCOS(H)
X      +0.006197*DSIN(2*H);
X
X  a_pert = - 0.003825*DCOS(H);
X
X  L_ += A;
X  M_ += A - B / e_;
X  e_ += e_pert;
X  a_ += a_pert;
X  omega_ += B;
X
X  /* Calculate E and nu */
X  anom_calc(M_, e_, &E, &nu);
X  r = a_ * (1 - e_ * DCOS(E));
X
X
X  u = L_ + nu - M_ - Omega_;
X  *C_p = nu - M_;
X  l = into_range(RAD_TO_DEG * atan2(DCOS(i_) * DSIN(u), DCOS(u)) + Omega_);
X  b = RAD_TO_DEG * asin(DSIN(u)*DSIN(i_));
X
X  l_pert = (0.010122 - 0.000988*v)*DSIN(S+eta)
X      +(-0.038581 + 0.002031*v - 0.001910*v*v)*DCOS(S+eta)
X    +(0.034964 - 0.001038*v + 0.000868*v*v)*DCOS(2*S+eta)
X      +0.005594*DSIN(S +3*theta);
X  l_pert += -0.014808*DSIN(zeta)
X      -0.005794*DSIN(eta)
X    +0.002347*DCOS(eta)
X      +0.009872*DSIN(theta)
X        +0.008803*DSIN(2*theta)
X          -0.004308*DSIN(3*theta);
X  b_pert = (0.000458*DSIN(eta) - 0.000642*DCOS(eta) - 0.000517*DCOS(4*theta))
X      *DSIN(S)
X    -(0.000347*DSIN(eta) + 0.000853*DCOS(eta) + 0.000517*DSIN(4*eta))
X      *DCOS(S)
X        +0.000403*(DCOS(2*theta)*DSIN(2*S) + DSIN(2*theta)*DCOS(2*S));
X  r_pert = -.025948
X      +.004985*DCOS(zeta)
X    -.001230*DCOS(S)
X      +.003354*DCOS(eta)
X        +(.005795*DCOS(S) - .001165*DSIN(S) + .001388*DCOS(2*S))*DSIN(eta)
X          +(.001351*DCOS(S) + .005702*DSIN(S) + .001388*DSIN(2*S))*DCOS(eta)
X        +.000904*DCOS(2*theta)
X          +.000894*(DCOS(theta) - DCOS(3*theta));
X
X  *r_p = r + r_pert;
X  *l_p = l + l_pert;
X  *b_p = b + b_pert;
X}
X
Xstatic void neptune(jd, L_, a_, e_, i_, omega_, Omega_, M_, r_p, l_p, b_p, C_p)
Xdouble jd, L_, a_, e_, i_, omega_, Omega_, M_;
Xdouble *r_p, *l_p, *b_p, *C_p;
X{
X  double E, nu;
X  double M, M1, M2, M4, M5, M6;
X  double r, l, b;
X
X  double u;            /* argument of latitude */
X  double T;
X  double A, B, e_pert, a_pert, b_pert, v, zeta, eta, theta;
X  double P, Q, S, V, W, G, H;
X  double r_pert, l_pert;
X
X  M = M_sun(jd);
X  M1 = M_1(jd);
X  M2 = M_2(jd);
X  M4 = M_4(jd);
X  M5 = M_5(jd);
X  M6 = M_6(jd);
X
X  /* for perturbations */
X  T = (jd - 2415020.0)/36525.0;
X
X  v = T/5.0 + 0.1;
X  P = 237.47555 +3034.9061*T;
X  Q = 265.91650 + 1222.1139*T;
X  S = 243.51721 + 428.4677*T;
X  V = 5.0*Q -2.0*P;
X  W = 2.0*P - 6.0*Q +3.0*S;
X  G = 83.76922 + 218.4901*T;
X  H = 2.0*G - S;
X  zeta = S - P;
X  eta = S - Q;
X  theta = G - S;
X
X  A = (-0.589833 + 0.001089*v)*DSIN(H)
X    +(-0.056094 + 0.004658*v)*DCOS(H)
X      -0.024286*DSIN(2*H);
X
X  e_pert = .0004389*DSIN(H)
X    +.0004262*DCOS(H)
X      +.0001129*DSIN(2*H)
X    +.0001089*DCOS(2*H);
X
X  B = 0.024039*DSIN(H)
X    -0.025303*DCOS(H)
X      +0.006206*DSIN(2*H)
X    -0.005992*DCOS(2*H);
X
X  a_pert = -0.000817*DSIN(H)
X    +0.008189*DCOS(H)
X      +0.000781*DCOS(2*H);
X
X  L_ += A;
X  M_ += A - B / e_;
X  e_ += e_pert;
X  a_ += a_pert;
X  omega_ += B;
X
X  /* Calculate E and nu */
X  anom_calc(M_, e_, &E, &nu);
X  r = a_ * (1 - e_ * DCOS(E));
X
X
X  u = L_ + nu - M_ - Omega_;
X  *C_p = nu - M_;
X  l = into_range(RAD_TO_DEG * atan2(DCOS(i_) * DSIN(u), DCOS(u)) + Omega_);
X  b = RAD_TO_DEG * asin(DSIN(u)*DSIN(i_));
X
X  l_pert = -0.009556*DSIN(zeta)
X      -0.005178*DSIN(eta)
X    +0.002572*DSIN(2*theta)
X      -0.002972*DCOS(2*theta)*DSIN(G)
X        -0.002833*DSIN(2*theta)*DCOS(G);
X
X  b_pert = 0.000336*DCOS(2*theta)*DSIN(G)
X      +0.000364*DSIN(2*theta)*DCOS(G);
X
X  r_pert = -.040596
X      +.004992*DCOS(zeta)
X    +.002744*DCOS(eta)
X      +.002044*DCOS(theta)
X        +.001051*DCOS(2*theta);
X
X  *r_p = r + r_pert;
X  *l_p = l + l_pert;
X  *b_p = b + b_pert;
X}
X
END_OF_FILE
if test 28915 -ne `wc -c <'observe/planetcalc.c'`; then
    echo shar: \"'observe/planetcalc.c'\" unpacked with wrong size!
fi
# end of 'observe/planetcalc.c'
fi
echo shar: End of archive 15 \(of 32\).
cp /dev/null ark15isdone
MISSING=""
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
    if test ! -f ark${I}isdone ; then
    MISSING="${MISSING} ${I}"
    fi
done
if test "${MISSING}" = "" ; then
    echo You have unpacked all 32 archives.
    rm -f ark[1-9]isdone ark[1-9][0-9]isdone
else
    echo You still need to unpack the following archives:
    echo "        " ${MISSING}
fi
##  End of shell archive.
exit 0