ccount@ATHENA.MIT.EDU (03/26/90)
Posting-number: Volume 11, Issue 43
Submitted-by: ccount@ATHENA.MIT.EDU
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