allbery@uunet.UU.NET (Brandon S. Allbery - comp.sources.misc) (11/28/89)
Posting-number: Volume 9, Issue 41
Submitted-by: ecd@umn-cs.cs.umn.edu@ncs-med.UUCP (Elwood C. Downey)
Archive-name: ephem2/part04
# This is a "shell archive" file; run it with sh.
# This is file 4.
echo x plans.c
cat > plans.c << 'xXx'
#include <stdio.h>
#include <math.h>
#include "astro.h"
#define TWOPI (2*PI)
#define mod2PI(x) ((x) - (int)((x)/TWOPI)*TWOPI)
/* given a modified Julian date, mjd, and a planet, p, find:
* lpd0: heliocentric longitude,
* psi0: heliocentric latitude,
* rp0: distance from the sun to the planet,
* rho0: distance from the Earth to the planet,
* none corrected for light time, ie, they are the true values for the
* given instant.
* lam: geocentric ecliptic longitude,
* bet: geocentric ecliptic latitude,
* each corrected for light time, ie, they are the apparent values as
* seen from the center of the Earth for the given instant.
* dia: angular diameter in arcsec at 1 AU,
* mag: visual magnitude when 1 AU from sun and earth at 0 phase angle.
*
* all angles are in radians, all distances in AU.
* the mean orbital elements are found by calling pelement(), then mutual
* perturbation corrections are applied as necessary.
*
* corrections for nutation and abberation must be made by the caller. The RA
* and DEC calculated from the fully-corrected ecliptic coordinates are then
* the apparent geocentric coordinates. Further corrections can be made, if
* required, for atmospheric refraction and geocentric parallax although the
* intrinsic error herein of about 10 arcseconds is usually the dominant
* error at this stage.
* TODO: combine the several intermediate expressions when get a good compiler.
*/
plans (mjd, p, lpd0, psi0, rp0, rho0, lam, bet, dia, mag)
double mjd;
int p;
double *lpd0, *psi0, *rp0, *rho0, *lam, *bet, *dia, *mag;
{
static double plan[8][9];
static lastmjd = -10000;
double dl; /* perturbation correction for longitude */
double dr; /* " orbital radius */
double dml; /* " mean longitude */
double ds; /* " eccentricity */
double dm; /* " mean anomaly */
double da; /* " semi-major axis */
double dhl; /* " heliocentric longitude */
double lsn, rsn; /* true geocentric longitude of sun and sun-earth rad */
double mas; /* mean anomaly of the sun */
double re; /* radius of earth's orbit */
double lg; /* longitude of earth */
double map[8]; /* array of mean anomalies for each planet */
double lpd, psi, rp, rho;
double ll, sll, cll;
double t;
double dt;
int pass;
int j;
double s, ma;
double nu, ea;
double lp, om;
double lo, slo, clo;
double inc, y;
double spsi, cpsi;
double rpd;
/* only need to fill in plan[] once for a given mjd */
if (mjd != lastmjd) {
pelement (mjd, plan);
lastmjd = mjd;
}
dt = 0;
t = mjd/36525.;
sunpos (mjd, &lsn, &rsn);
masun (mjd, &mas);
re = rsn;
lg = lsn+PI;
/* first find the true position of the planet at mjd.
* then repeat a second time for a slightly different time based
* on the position found in the first pass to account for light-travel
* time.
*/
for (pass = 0; pass < 2; pass++) {
for (j = 0; j < 8; j++)
map[j] = degrad(plan[j][0]-plan[j][2]-dt*plan[j][1]);
/* set initial corrections to 0.
* then modify as necessary for the planet of interest.
*/
dl = 0;
dr = 0;
dml = 0;
ds = 0;
dm = 0;
da = 0;
dhl = 0;
switch (p) {
case MERCURY:
p_mercury (map, &dl, &dr);
break;
case VENUS:
p_venus (t, mas, map, &dl, &dr, &dml, &dm);
break;
case MARS:
p_mars (mas, map, &dl, &dr, &dml, &dm);
break;
case JUPITER:
p_jupiter (t, plan[p][3], &dml, &ds, &dm, &da);
break;
case SATURN:
p_saturn (t, plan[p][3], &dml, &ds, &dm, &da, &dhl);
break;
case URANUS:
p_uranus (t, plan[p][3], &dl, &dr, &dml, &ds, &dm, &da, &dhl);
break;
case NEPTUNE:
p_neptune (t, plan[p][3], &dl, &dr, &dml, &ds, &dm, &da, &dhl);
break;
case PLUTO:
/* no perturbation theory for pluto */
break;
}
s = plan[p][3]+ds;
ma = map[p]+dm;
anomaly (ma, s, &nu, &ea);
rp = (plan[p][6]+da)*(1-s*s)/(1+s*cos(nu));
lp = raddeg(nu)+plan[p][2]+raddeg(dml-dm);
lp = degrad(lp);
om = degrad(plan[p][5]);
lo = lp-om;
slo = sin(lo);
clo = cos(lo);
inc = degrad(plan[p][4]);
rp = rp+dr;
spsi = slo*sin(inc);
y = slo*cos(inc);
psi = asin(spsi)+dhl;
spsi = sin(psi);
lpd = atan(y/clo)+om+degrad(dl);
if (clo<0) lpd += PI;
if (lpd>TWOPI) lpd -= TWOPI;
cpsi = cos(psi);
rpd = rp*cpsi;
ll = lpd-lg;
rho = sqrt(re*re+rp*rp-2*re*rp*cpsi*cos(ll));
/* when we view a planet we see it in the position it occupied
* dt hours ago, where rho is the distance between it and earth,
* in AU. use this as the new time for the next pass.
*/
dt = rho*5.775518e-3;
if (pass == 0) {
/* save heliocentric coordinates after first pass since, being
* true, they are NOT to be corrected for light-travel time.
*/
*lpd0 = lpd;
*psi0 = psi;
*rp0 = rp;
*rho0 = rho;
}
}
sll = sin(ll);
cll = cos(ll);
if (p < MARS)
*lam = atan(-1*rpd*sll/(re-rpd*cll))+lg+PI;
else
*lam = atan(re*sll/(rpd-re*cll))+lpd;
range (lam, TWOPI);
*bet = atan(rpd*spsi*sin(*lam-lpd)/(cpsi*re*sll));
*dia = plan[p][7];
*mag = plan[p][8];
}
/* set auxilliary variables used for jupiter, saturn, uranus, and neptune */
static
aux_jsun (t, j1, j2, j3, j4, j5, j6)
double t;
double *j1, *j2, *j3, *j4, *j5, *j6;
{
*j1 = t/5+0.1;
*j2 = mod2PI(4.14473+5.29691e1*t);
*j3 = mod2PI(4.641118+2.132991e1*t);
*j4 = mod2PI(4.250177+7.478172*t);
*j5 = 5 * *j3 - 2 * *j2;
*j6 = 2 * *j2 - 6 * *j3 + 3 * *j4;
}
/* find the mean anomaly of the sun at mjd.
* this is the same as that used in sun() but when it was converted to C it
* was not known it would be required outside that routine.
* TODO: add an argument to sun() to return mas and eliminate this routine.
*/
static
masun (mjd, mas)
double mjd;
double *mas;
{
double t, t2;
double a, b;
t = mjd/36525;
t2 = t*t;
a = 9.999736042e1*t;
b = 360.*(a-(int)a);
*mas = degrad (3.5847583e2-(1.5e-4+3.3e-6*t)*t2+b);
}
/* perturbations for mercury */
static
p_mercury (map, dl, dr)
double map[];
double *dl, *dr;
{
*dl = 2.04e-3*cos(5*map[2-1]-2*map[1-1]+2.1328e-1)+
1.03e-3*cos(2*map[2-1]-map[1-1]-2.8046)+
9.1e-4*cos(2*map[3]-map[1-1]-6.4582e-1)+
7.8e-4*cos(5*map[2-1]-3*map[1-1]+1.7692e-1);
*dr = 7.525e-6*cos(2*map[3]-map[1-1]+9.25251e-1)+
6.802e-6*cos(5*map[2-1]-3*map[1-1]-4.53642)+
5.457e-6*cos(2*map[2-1]-2*map[1-1]-1.24246)+
3.569e-6*cos(5*map[2-1]-map[1-1]-1.35699);
}
/* ....venus */
static
p_venus (t, mas, map, dl, dr, dml, dm)
double t, mas, map[];
double *dl, *dr, *dml, *dm;
{
*dml = degrad (7.7e-4*sin(4.1406+t*2.6227));
*dm = *dml;
*dl = 3.13e-3*cos(2*mas-2*map[2-1]-2.587)+
1.98e-3*cos(3*mas-3*map[2-1]+4.4768e-2)+
1.36e-3*cos(mas-map[2-1]-2.0788)+
9.6e-4*cos(3*mas-2*map[2-1]-2.3721)+
8.2e-4*cos(map[3]-map[2-1]-3.6318);
*dr = 2.2501e-5*cos(2*mas-2*map[2-1]-1.01592)+
1.9045e-5*cos(3*mas-3*map[2-1]+1.61577)+
6.887e-6*cos(map[3]-map[2-1]-2.06106)+
5.172e-6*cos(mas-map[2-1]-5.08065e-1)+
3.62e-6*cos(5*mas-4*map[2-1]-1.81877)+
3.283e-6*cos(4*mas-4*map[2-1]+1.10851)+
3.074e-6*cos(2*map[3]-2*map[2-1]-9.62846e-1);
}
/* ....mars */
static
p_mars (mas, map, dl, dr, dml, dm)
double mas, map[];
double *dl, *dr, *dml, *dm;
{
double a;
a = 3*map[3]-8*map[2]+4*mas;
*dml = degrad (-1*(1.133e-2*sin(a)+9.33e-3*cos(a)));
*dm = *dml;
*dl = 7.05e-3*cos(map[3]-map[2]-8.5448e-1)+
6.07e-3*cos(2*map[3]-map[2]-3.2873)+
4.45e-3*cos(2*map[3]-2*map[2]-3.3492)+
3.88e-3*cos(mas-2*map[2]+3.5771e-1)+
2.38e-3*cos(mas-map[2]+6.1256e-1)+
2.04e-3*cos(2*mas-3*map[2]+2.7688)+
1.77e-3*cos(3*map[2]-map[2-1]-1.0053)+
1.36e-3*cos(2*mas-4*map[2]+2.6894)+
1.04e-3*cos(map[3]+3.0749e-1);
*dr = 5.3227e-5*cos(map[3]-map[2]+7.17864e-1)+
5.0989e-5*cos(2*map[3]-2*map[2]-1.77997)+
3.8278e-5*cos(2*map[3]-map[2]-1.71617)+
1.5996e-5*cos(mas-map[2]-9.69618e-1)+
1.4764e-5*cos(2*mas-3*map[2]+1.19768)+
8.966e-6*cos(map[3]-2*map[2]+7.61225e-1);
*dr += 7.914e-6*cos(3*map[3]-2*map[2]-2.43887)+
7.004e-6*cos(2*map[3]-3*map[2]-1.79573)+
6.62e-6*cos(mas-2*map[2]+1.97575)+
4.93e-6*cos(3*map[3]-3*map[2]-1.33069)+
4.693e-6*cos(3*mas-5*map[2]+3.32665)+
4.571e-6*cos(2*mas-4*map[2]+4.27086)+
4.409e-6*cos(3*map[3]-map[2]-2.02158);
}
/* ....jupiter */
static
p_jupiter (t, s, dml, ds, dm, da)
double t, s;
double *dml, *ds, *dm, *da;
{
double dp;
double j1, j2, j3, j4, j5, j6, j7;
double sj3, cj3, s2j3, c2j3;
double sj5, cj5, s2j5;
double sj6;
double sj7, cj7, s2j7, c2j7, s3j7, c3j7, s4j7, c4j7, c5j7;
aux_jsun (t, &j1, &j2, &j3, &j4, &j5, &j6);
j7 = j3-j2;
sj3 = sin(j3);
cj3 = cos(j3);
s2j3 = sin(2*j3);
c2j3 = cos(2*j3);
sj5 = sin(j5);
cj5 = cos(j5);
s2j5 = sin(2*j5);
sj6 = sin(j6);
sj7 = sin(j7);
cj7 = cos(j7);
s2j7 = sin(2*j7);
c2j7 = cos(2*j7);
s3j7 = sin(3*j7);
c3j7 = cos(3*j7);
s4j7 = sin(4*j7);
c4j7 = cos(4*j7);
c5j7 = cos(5*j7);
*dml = (3.31364e-1-(1.0281e-2+4.692e-3*j1)*j1)*sj5+
(3.228e-3-(6.4436e-2-2.075e-3*j1)*j1)*cj5-
(3.083e-3+(2.75e-4-4.89e-4*j1)*j1)*s2j5+
2.472e-3*sj6+1.3619e-2*sj7+1.8472e-2*s2j7+6.717e-3*s3j7+
2.775e-3*s4j7+6.417e-3*s2j7*sj3+
(7.275e-3-1.253e-3*j1)*sj7*sj3+
2.439e-3*s3j7*sj3-(3.5681e-2+1.208e-3*j1)*sj7*cj3;
*dml += -3.767e-3*c2j7*sj3-(3.3839e-2+1.125e-3*j1)*cj7*sj3-
4.261e-3*s2j7*cj3+
(1.161e-3*j1-6.333e-3)*cj7*cj3+
2.178e-3*cj3-6.675e-3*c2j7*cj3-2.664e-3*c3j7*cj3-
2.572e-3*sj7*s2j3-3.567e-3*s2j7*s2j3+2.094e-3*cj7*c2j3+
3.342e-3*c2j7*c2j3;
*dml = degrad(*dml);
*ds = (3606+(130-43*j1)*j1)*sj5+(1289-580*j1)*cj5-6764*sj7*sj3-
1110*s2j7*sj3-224*s3j7*sj3-204*sj3+(1284+116*j1)*cj7*sj3+
188*c2j7*sj3+(1460+130*j1)*sj7*cj3+224*s2j7*cj3-817*cj3+
6074*cj3*cj7+992*c2j7*cj3+
508*c3j7*cj3+230*c4j7*cj3+108*c5j7*cj3;
*ds += -(956+73*j1)*sj7*s2j3+448*s2j7*s2j3+137*s3j7*s2j3+
(108*j1-997)*cj7*s2j3+480*c2j7*s2j3+148*c3j7*s2j3+
(99*j1-956)*sj7*c2j3+490*s2j7*c2j3+
158*s3j7*c2j3+179*c2j3+(1024+75*j1)*cj7*c2j3-
437*c2j7*c2j3-132*c3j7*c2j3;
*ds *= 1e-7;
dp = (7.192e-3-3.147e-3*j1)*sj5-4.344e-3*sj3+
(j1*(1.97e-4*j1-6.75e-4)-2.0428e-2)*cj5+
3.4036e-2*cj7*sj3+(7.269e-3+6.72e-4*j1)*sj7*sj3+
5.614e-3*c2j7*sj3+2.964e-3*c3j7*sj3+3.7761e-2*sj7*cj3+
6.158e-3*s2j7*cj3-
6.603e-3*cj7*cj3-5.356e-3*sj7*s2j3+2.722e-3*s2j7*s2j3+
4.483e-3*cj7*s2j3-2.642e-3*c2j7*s2j3+4.403e-3*sj7*c2j3-
2.536e-3*s2j7*c2j3+5.547e-3*cj7*c2j3-2.689e-3*c2j7*c2j3;
*dm = *dml-(degrad(dp)/s);
*da = 205*cj7-263*cj5+693*c2j7+312*c3j7+147*c4j7+299*sj7*sj3+
181*c2j7*sj3+204*s2j7*cj3+111*s3j7*cj3-337*cj7*cj3-
111*c2j7*cj3;
*da *= 1e-6;
}
/* ....saturn */
static
p_saturn (t, s, dml, ds, dm, da, dhl)
double t, s;
double *dml, *ds, *dm, *da, *dhl;
{
double dp;
double j1, j2, j3, j4, j5, j6, j7, j8;
double sj3, cj3, s2j3, c2j3, s3j3, c3j3, s4j3, c4j3;
double sj5, cj5, s2j5, c2j5;
double sj6;
double sj7, cj7, s2j7, c2j7, s3j7, c3j7, s4j7, c4j7, c5j7, s5j7;
double s2j8, c2j8, s3j8, c3j8;
aux_jsun (t, &j1, &j2, &j3, &j4, &j5, &j6);
j7 = j3-j2;
sj3 = sin(j3);
cj3 = cos(j3);
s2j3 = sin(2*j3);
c2j3 = cos(2*j3);
sj5 = sin(j5);
cj5 = cos(j5);
s2j5 = sin(2*j5);
sj6 = sin(j6);
sj7 = sin(j7);
cj7 = cos(j7);
s2j7 = sin(2*j7);
c2j7 = cos(2*j7);
s3j7 = sin(3*j7);
c3j7 = cos(3*j7);
s4j7 = sin(4*j7);
c4j7 = cos(4*j7);
c5j7 = cos(5*j7);
s3j3 = sin(3*j3);
c3j3 = cos(3*j3);
s4j3 = sin(4*j3);
c4j3 = cos(4*j3);
c2j5 = cos(2*j5);
s5j7 = sin(5*j7);
j8 = j4-j3;
s2j8 = sin(2*j8);
c2j8 = cos(2*j8);
s3j8 = sin(3*j8);
c3j8 = cos(3*j8);
*dml = 7.581e-3*s2j5-7.986e-3*sj6-1.48811e-1*sj7-4.0786e-2*s2j7-
(8.14181e-1-(1.815e-2-1.6714e-2*j1)*j1)*sj5-
(1.0497e-2-(1.60906e-1-4.1e-3*j1)*j1)*cj5-1.5208e-2*s3j7-
6.339e-3*s4j7-6.244e-3*sj3-1.65e-2*s2j7*sj3+
(8.931e-3+2.728e-3*j1)*sj7*sj3-5.775e-3*s3j7*sj3+
(8.1344e-2+3.206e-3*j1)*cj7*sj3+1.5019e-2*c2j7*sj3;
*dml += (8.5581e-2+2.494e-3*j1)*sj7*cj3+1.4394e-2*c2j7*cj3+
(2.5328e-2-3.117e-3*j1)*cj7*cj3+
6.319e-3*c3j7*cj3+6.369e-3*sj7*s2j3+9.156e-3*s2j7*s2j3+
7.525e-3*s3j8*s2j3-5.236e-3*cj7*c2j3-7.736e-3*c2j7*c2j3-
7.528e-3*c3j8*c2j3;
*dml = degrad(*dml);
*ds = (-7927+(2548+91*j1)*j1)*sj5+(13381+(1226-253*j1)*j1)*cj5+
(248-121*j1)*s2j5-(305+91*j1)*c2j5+412*s2j7+12415*sj3+
(390-617*j1)*sj7*sj3+(165-204*j1)*s2j7*sj3+26599*cj7*sj3-
4687*c2j7*sj3-1870*c3j7*sj3-821*c4j7*sj3-
377*c5j7*sj3+497*c2j8*sj3+(163-611*j1)*cj3;
*ds += -12696*sj7*cj3-4200*s2j7*cj3-1503*s3j7*cj3-619*s4j7*cj3-
268*s5j7*cj3-(282+1306*j1)*cj7*cj3+(-86+230*j1)*c2j7*cj3+
461*s2j8*cj3-350*s2j3+(2211-286*j1)*sj7*s2j3-
2208*s2j7*s2j3-568*s3j7*s2j3-346*s4j7*s2j3-
(2780+222*j1)*cj7*s2j3+(2022+263*j1)*c2j7*s2j3+248*c3j7*s2j3+
242*s3j8*s2j3+467*c3j8*s2j3-490*c2j3-(2842+279*j1)*sj7*c2j3;
*ds += (128+226*j1)*s2j7*c2j3+224*s3j7*c2j3+
(-1594+282*j1)*cj7*c2j3+(2162-207*j1)*c2j7*c2j3+
561*c3j7*c2j3+343*c4j7*c2j3+469*s3j8*c2j3-242*c3j8*c2j3-
205*sj7*s3j3+262*s3j7*s3j3+208*cj7*c3j3-271*c3j7*c3j3-
382*c3j7*s4j3-376*s3j7*c4j3;
*ds *= 1e-7;
dp = (7.7108e-2+(7.186e-3-1.533e-3*j1)*j1)*sj5-7.075e-3*sj7+
(4.5803e-2-(1.4766e-2+5.36e-4*j1)*j1)*cj5-7.2586e-2*cj3-
7.5825e-2*sj7*sj3-2.4839e-2*s2j7*sj3-8.631e-3*s3j7*sj3-
1.50383e-1*cj7*cj3+2.6897e-2*c2j7*cj3+1.0053e-2*c3j7*cj3-
(1.3597e-2+1.719e-3*j1)*sj7*s2j3+1.1981e-2*s2j7*c2j3;
dp += -(7.742e-3-1.517e-3*j1)*cj7*s2j3+
(1.3586e-2-1.375e-3*j1)*c2j7*c2j3-
(1.3667e-2-1.239e-3*j1)*sj7*c2j3+
(1.4861e-2+1.136e-3*j1)*cj7*c2j3-
(1.3064e-2+1.628e-3*j1)*c2j7*c2j3;
*dm = *dml-(degrad(dp)/s);
*da = 572*sj5-1590*s2j7*cj3+2933*cj5-647*s3j7*cj3+33629*cj7-
344*s4j7*cj3-3081*c2j7+2885*cj7*cj3-1423*c3j7+
(2172+102*j1)*c2j7*cj3-671*c4j7+296*c3j7*cj3-320*c5j7-
267*s2j7*s2j3+1098*sj3-778*cj7*s2j3-2812*sj7*sj3;
*da += 495*c2j7*s2j3+688*s2j7*sj3+250*c3j7*s2j3-393*s3j7*sj3-
856*sj7*c2j3-228*s4j7*sj3+441*s2j7*c2j3+2138*cj7*sj3+
296*c2j7*c2j3-999*c2j7*sj3+211*c3j7*c2j3-642*c3j7*sj3-
427*sj7*s3j3-325*c4j7*sj3+398*s3j7*s3j3-890*cj3+
344*cj7*c3j3+2206*sj7*cj3-427*c3j7*c3j3;
*da *= 1e-6;
*dhl = 7.47e-4*cj7*sj3+1.069e-3*cj7*cj3+2.108e-3*s2j7*s2j3+
1.261e-3*c2j7*s2j3+1.236e-3*s2j7*c2j3-2.075e-3*c2j7*c2j3;
*dhl = degrad(*dhl);
}
/* ....uranus */
static
p_uranus (t, s, dl, dr, dml, ds, dm, da, dhl)
double t, s;
double *dl, *dr, *dml, *ds, *dm, *da, *dhl;
{
double dp;
double j1, j2, j3, j4, j5, j6;
double j8, j9, j10, j11, j12;
double sj4, cj4, s2j4, c2j4;
double sj9, cj9, s2j9, c2j9;
double sj11, cj11;
aux_jsun (t, &j1, &j2, &j3, &j4, &j5, &j6);
j8 = mod2PI(1.46205+3.81337*t);
j9 = 2*j8-j4;
sj9 = sin(j9);
cj9 = cos(j9);
s2j9 = sin(2*j9);
c2j9 = cos(2*j9);
j10 = j4-j2;
j11 = j4-j3;
j12 = j8-j4;
*dml = (8.64319e-1-1.583e-3*j1)*sj9+(8.2222e-2-6.833e-3*j1)*cj9+
3.6017e-2*s2j9-3.019e-3*c2j9+8.122e-3*sin(j6);
*dml = degrad(*dml);
dp = 1.20303e-1*sj9+6.197e-3*s2j9+(1.9472e-2-9.47e-4*j1)*cj9;
*dm = *dml-(degrad(dp)/s);
*ds = (163*j1-3349)*sj9+20981*cj9+1311*c2j9;
*ds *= 1e-7;
*da = -3.825e-3*cj9;
*dl = (1.0122e-2-9.88e-4*j1)*sin(j4+j11)+
(-3.8581e-2+(2.031e-3-1.91e-3*j1)*j1)*cos(j4+j11)+
(3.4964e-2-(1.038e-3-8.68e-4*j1)*j1)*cos(2*j4+j11)+
5.594e-3*sin(j4+3*j12)-1.4808e-2*sin(j10)-
5.794e-3*sin(j11)+2.347e-3*cos(j11)+9.872e-3*sin(j12)+
8.803e-3*sin(2*j12)-4.308e-3*sin(3*j12);
sj11 = sin(j11);
cj11 = cos(j11);
sj4 = sin(j4);
cj4 = cos(j4);
s2j4 = sin(2*j4);
c2j4 = cos(2*j4);
*dhl = (4.58e-4*sj11-6.42e-4*cj11-5.17e-4*cos(4*j12))*sj4-
(3.47e-4*sj11+8.53e-4*cj11+5.17e-4*sin(4*j11))*cj4+
4.03e-4*(cos(2*j12)*s2j4+sin(2*j12)*c2j4);
*dhl = degrad(*dhl);
*dr = -25948+4985*cos(j10)-1230*cj4+3354*cos(j11)+904*cos(2*j12)+
894*(cos(j12)-cos(3*j12))+(5795*cj4-1165*sj4+1388*c2j4)*sj11+
(1351*cj4+5702*sj4+1388*s2j4)*cos(j11);
*dr *= 1e-6;
}
/* ....neptune */
static
p_neptune (t, s, dl, dr, dml, ds, dm, da, dhl)
double t, s;
double *dl, *dr, *dml, *ds, *dm, *da, *dhl;
{
double dp;
double j1, j2, j3, j4, j5, j6;
double j8, j9, j10, j11, j12;
double sj8, cj8;
double sj9, cj9, s2j9, c2j9;
double s2j12, c2j12;
aux_jsun (t, &j1, &j2, &j3, &j4, &j5, &j6);
j8 = mod2PI(1.46205+3.81337*t);
j9 = 2*j8-j4;
sj9 = sin(j9);
cj9 = cos(j9);
s2j9 = sin(2*j9);
c2j9 = cos(2*j9);
j10 = j8-j2;
j11 = j8-j3;
j12 = j8-j4;
*dml = (1.089e-3*j1-5.89833e-1)*sj9+(4.658e-3*j1-5.6094e-2)*cj9-
2.4286e-2*s2j9;
*dml = degrad(*dml);
dp = 2.4039e-2*sj9-2.5303e-2*cj9+6.206e-3*s2j9-5.992e-3*c2j9;
*dm = *dml-(degrad(dp)/s);
*ds = 4389*sj9+1129*s2j9+4262*cj9+1089*c2j9;
*ds *= 1e-7;
*da = 8189*cj9-817*sj9+781*c2j9;
*da *= 1e-6;
s2j12 = sin(2*j12);
c2j12 = cos(2*j12);
sj8 = sin(j8);
cj8 = cos(j8);
*dl = -9.556e-3*sin(j10)-5.178e-3*sin(j11)+2.572e-3*s2j12-
2.972e-3*c2j12*sj8-2.833e-3*s2j12*cj8;
*dhl = 3.36e-4*c2j12*sj8+3.64e-4*s2j12*cj8;
*dhl = degrad(*dhl);
*dr = -40596+4992*cos(j10)+2744*cos(j11)+2044*cos(j12)+1051*c2j12;
*dr *= 1e-6;
}
xXx
echo x plot.c
cat > plot.c << 'xXx'
/* code to support the plotting capabilities.
* idea is to let the operator name a plot file and mark some fields for
* logging. then after each screen update, the logged fields are written to
* the plot file. later, the file may be plotted (very simplistically by
* ephem, for now anyway, or by some other program entirely.).
*
* format of the plot file is one line per coordinate: label,x,y
* if z was specified, it is a fourth field.
* x,y,z are plotted using %g format.
*/
#include <stdio.h>
#include <math.h>
#include "screen.h"
extern errno;
extern char *sys_errlist[];
#define errsys (sys_errlist[errno])
#define TRACE(x) {FILE *fp = fopen("trace","a"); fprintf x; fclose(fp);}
#define MAXPLTLINES 10 /* max number of labeled lines we can track.
* note we can't store more than NFLOGS fields
* anyway (see flog.c).
*/
#define FNLEN (14+1) /* longest filename; plus 1 for \0 */
static char plt_filename[FNLEN] = "ephem.plt"; /* default plot file name */
static FILE *plt_fp; /* the plot file; == 0 means don't plot */
/* store the label and rcfpack()s for each line to track. */
typedef struct {
char pl_label;
int pl_rcpx, pl_rcpy, pl_rcpz;
} PltLine;
static PltLine pltlines[MAXPLTLINES];
static int npltlines; /* number of pltlines[] in actual use */
static int plt_in_polar; /*if true plot in polar coords, else cartesian*/
/* picked the Plot label:
* if on, just turn it off.
* if off, turn on, define fields or select name of file to plot and do it.
* TODO: more flexability, more relevance.
*/
plot_setup()
{
if (plt_fp) {
plt_turn_off();
plot_prstate (0);
} else {
static char *chcs[4] = {
"Select fields", "Display a plot file", (char *)0,
"Begin plotting"
};
int ff = 0;
ask:
chcs[2] = plt_in_polar ? "Polar coords" : "Cartesian coords";
switch (popup(chcs, ff, npltlines > 0 ? 4 : 3)) {
case 0: plt_select_fields(); break;
case 1: plt_file(); break;
case 2: plt_in_polar ^= 1; ff = 2; goto ask;
case 3: plt_turn_on(); break;
}
}
}
/* write the active plotfields to the current plot file, if one is open. */
plot()
{
if (plt_fp) {
PltLine *plp;
double x, y, z;
for (plp = &pltlines[npltlines]; --plp >= pltlines; ) {
if (flog_get (plp->pl_rcpx, &x) == 0
&& flog_get (plp->pl_rcpy, &y) == 0) {
fprintf (plt_fp, "%c,%.12g,%.12g", plp->pl_label, x, y);
if (flog_get (plp->pl_rcpz, &z) == 0)
fprintf (plt_fp, ",%.12g", z);
fprintf (plt_fp, "\n");
}
}
}
}
plot_prstate (force)
int force;
{
static last;
int this = plt_fp != 0;
if (force || this != last) {
f_string (R_PLOT, C_PLOTV, this ? " on" : "off");
last = this;
}
}
plot_ison()
{
return (plt_fp != 0);
}
static
plt_reset()
{
PltLine *plp;
for (plp = &pltlines[npltlines]; --plp >= pltlines; ) {
(void) flog_delete (plp->pl_rcpx);
(void) flog_delete (plp->pl_rcpy);
(void) flog_delete (plp->pl_rcpz);
plp->pl_rcpx = plp->pl_rcpy = plp->pl_rcpz = 0;
}
npltlines = 0;
}
/* let operator select the fields he wants to plot.
* register them with flog and keep rcfpack() in pltlines[] array.
*/
static
plt_select_fields()
{
static char hlp[] = "move and RETURN to select a field, or q to quit";
static char sry[] = "Sorry; can not log any more fields.";
int r = R_UT, c = C_UTV; /* TODO: start where main was? */
char buf[64];
int rcp;
int i;
plt_reset();
for (i = 0; i < MAXPLTLINES; i++) {
sprintf (buf, "select x field for line %d", i+1);
rcp = sel_fld (r, c, alt_menumask()|F_PLT, buf, hlp);
if (!rcp)
break;
pltlines[i].pl_rcpx = rcp;
if (flog_add (rcp) < 0) {
f_msg (sry);
break;
}
sprintf (buf, "select y field for line %d", i+1);
r = unpackr (rcp); c = unpackc (rcp);
rcp = sel_fld (r, c, alt_menumask()|F_PLT, buf, hlp);
if (!rcp) {
(void) flog_delete (pltlines[i].pl_rcpx);
break;
}
if (flog_add (rcp) < 0) {
(void) flog_delete (pltlines[i].pl_rcpx);
f_msg (sry);
break;
}
pltlines[i].pl_rcpy = rcp;
sprintf (buf, "select z field for line %d", i+1);
r = unpackr (rcp); c = unpackc (rcp);
rcp = sel_fld (r, c, alt_menumask()|F_PLT, buf, hlp);
if (rcp) {
if (flog_add (rcp) < 0) {
(void) flog_delete (pltlines[i].pl_rcpx);
(void) flog_delete (pltlines[i].pl_rcpy);
f_msg (sry);
break;
}
pltlines[i].pl_rcpz = rcp;
r = unpackr (rcp); c = unpackc (rcp);
}
do {
sprintf (buf, "enter a one-character label for line %d: ", i+1);
f_prompt (buf);
} while (read_line (buf, 1) != 1);
pltlines[i].pl_label = *buf;
}
npltlines = i;
}
static
plt_turn_off ()
{
fclose (plt_fp);
plt_fp = 0;
plot_prstate(0);
}
/* turn on plotting.
* establish a file to use (and thereby set plt_fp, the plotting_is_on flag).
* also check that there is a srch function if it is being plotted.
*/
static
plt_turn_on ()
{
int sf = rcfpack(R_SRCH, C_SRCH, 0);
char fn[FNLEN], fnq[64];
char *optype;
int n;
PltLine *plp;
/* insure there is a valid srch function if we are to plot it */
for (plp = &pltlines[npltlines]; --plp >= pltlines; )
if ((plp->pl_rcpx == sf || plp->pl_rcpy == sf || plp->pl_rcpz == sf)
&& !prog_isgood()) {
f_msg ("Plotting search function but it is not defined.");
return;
}
/* prompt for file name, giving current as default */
sprintf (fnq, "file to write <%s>: ", plt_filename);
f_prompt (fnq);
n = read_line (fn, sizeof(fn)-1);
/* leave plotting off if type END.
* reuse same fn if just type \n
*/
if (n < 0)
return;
if (n > 0)
strcpy (plt_filename, fn);
/* give option to append if file already exists */
optype = "w";
if (access (plt_filename, 2) == 0) {
while (1) {
f_prompt ("files exists; append or overwrite (a/o)?: ");
n = read_char();
if (n == 'a') {
optype = "a";
break;
}
if (n == 'o')
break;
}
}
/* plotting is on if file opens ok */
plt_fp = fopen (plt_filename, optype);
if (!plt_fp) {
char buf[NC];
sprintf (buf, "can not open %s: %s", plt_filename, errsys);
f_prompt (buf);
(void)read_char();
}
plot_prstate (0);
}
/* ask operator for a file to plot. if it's ok, do it.
*/
static
plt_file ()
{
char fn[FNLEN], fnq[64];
FILE *pfp;
int n;
/* prompt for file name, giving current as default */
sprintf (fnq, "file to read <%s>: ", plt_filename);
f_prompt (fnq);
n = read_line (fn, sizeof(fn)-1);
/* forget it if type END.
* reuse same fn if just type \n
*/
if (n < 0)
return;
if (n > 0)
strcpy (plt_filename, fn);
/* do the plot if file opens ok */
pfp = fopen (plt_filename, "r");
if (pfp) {
if (plt_in_polar)
plot_polar (pfp);
else
plot_cartesian (pfp);
fclose (pfp);
} else {
char buf[NC];
sprintf (buf, "can not open %s: %s", plt_filename, errsys);
f_prompt (buf);
(void)read_char();
}
}
/* plot the given file on the screen in cartesian coords.
* TODO: add z tags somehow
* N.B. do whatever you like but redraw the screen when done.
*/
static
plot_cartesian (pfp)
FILE *pfp;
{
static char fmt[] = "%c,%F,%F";
double x, y; /* N.B. be sure these match what your scanf's %F wants*/
double minx, maxx, miny, maxy;
char buf[128];
int npts = 0;
char c;
/* find ranges and number of points */
while (fgets (buf, sizeof(buf), pfp)) {
sscanf (buf, fmt, &c, &x, &y);
if (npts++ == 0) {
maxx = minx = x;
maxy = miny = y;
} else {
if (x > maxx) maxx = x;
else if (x < minx) minx = x;
if (y > maxy) maxy = y;
else if (y < miny) miny = y;
}
}
#define SMALL (1e-10)
if (npts < 2 || fabs(minx-maxx) < SMALL || fabs(miny-maxy) < SMALL)
f_prompt ("At least two different points required to plot.");
else {
/* read file again, this time plotting */
rewind (pfp);
c_erase();
while (fgets (buf, sizeof(buf), pfp)) {
int row, col;
sscanf (buf, fmt, &c, &x, &y);
row = NR-(int)((NR-1)*(y-miny)/(maxy-miny)+0.5);
col = 1+(int)((NC-1)*(x-minx)/(maxx-minx)+0.5);
if (row == NR && col == NC)
col--; /* avoid lower right scrolling corner */
f_char (row, col, c);
}
/* label axes */
f_double (1, 1, "%.2g", maxy);
f_double (NR-1, 1, "%.2g", miny);
f_double (NR, 1, "%.2g", minx);
f_double (NR, NC-10, "%.2g", maxx);
}
/* hit any key to resume... */
(void) read_char();
redraw_screen (2); /* full redraw */
}
/* plot the given file on the screen in polar coords.
* first numberic field in plot file is r, second is theta in degrees.
* TODO: add z tags somehow
* N.B. do whatever you like but redraw the screen when done.
*/
static
plot_polar (pfp)
FILE *pfp;
{
static char fmt[] = "%c,%F,%F";
double r, th; /* N.B. be sure these match what your scanf's %F wants*/
double maxr;
char buf[128];
int npts = 0;
char c;
/* find ranges and number of points */
while (fgets (buf, sizeof(buf), pfp)) {
sscanf (buf, fmt, &c, &r, &th);
if (npts++ == 0)
maxr = r;
else
if (r > maxr)
maxr = r;
}
if (npts < 2)
f_prompt ("At least two points required to plot.");
else {
/* read file again, this time plotting */
rewind (pfp);
c_erase();
while (fgets (buf, sizeof(buf), pfp)) {
int row, col;
double x, y;
sscanf (buf, fmt, &c, &r, &th);
x = r * cos(th/57.2958); /* degs to rads */
y = r * sin(th/57.2958);
row = NR-(int)((NR-1)*(y+maxr)/(2.0*maxr)+0.5);
col = 1+(int)((NC-1)*(x+maxr)/(2.0*maxr)/ASPECT+0.5);
if (row == NR && col == NC)
col--; /* avoid lower right scrolling corner */
f_char (row, col, c);
}
/* label radius */
f_double (NR/2, NC-10, "%.4g", maxr);
}
/* hit any key to resume... */
(void) read_char();
redraw_screen (2); /* full redraw */
}
xXx
echo x popup.c
cat > popup.c << 'xXx'
/* put up a one-line menu consisting of the given fields and let op move
* between them with the same methods as sel_fld().
* return index of which he picked, or -1 if hit END.
*/
#include <stdio.h>
#include "screen.h"
#define FLDGAP 2 /* inter-field gap */
#define MAXFLDS 32 /* max number of fields we can handle */
static char pup[] = "Select: ";
/* put up an array of strings on prompt line and let op pick one.
* start with field fn.
* N.B. we do not do much error/bounds checking.
*/
popup (fields, fn, nfields)
char *fields[];
int fn;
int nfields;
{
int fcols[MAXFLDS]; /* column to use for each field */
int i;
if (nfields > MAXFLDS)
return (-1);
again:
/* erase the prompt line; we are going to take it over */
c_pos (R_PROMPT, C_PROMPT);
c_eol();
/* compute starting column for each field */
fcols[0] = sizeof(pup);
for (i = 1; i < nfields; i++)
fcols[i] = fcols[i-1] + strlen (fields[i-1]) + FLDGAP;
/* draw each field, with comma after all but last */
c_pos (R_PROMPT, 1);
fputs (pup, stdout);
for (i = 0; i < nfields; i++) {
c_pos (R_PROMPT, fcols[i]);
printf (i < nfields-1 ? "%s," : "%s", fields[i]);
}
/* let op choose one now; begin at fn.
*/
i = fn;
while (1) {
c_pos (R_PROMPT, fcols[i]);
switch (read_char()) {
case END: return (-1);
case REDRAW: redraw_screen(2); goto again;
case VERSION: version(); goto again;
case '\r': return (i);
case 'h':
if (--i < 0)
i = nfields - 1;
break;
case 'l':
if (++i >= nfields)
i = 0;
break;
}
}
}
xXx
echo x precess.c
cat > precess.c << 'xXx'
#include <stdio.h>
#include <math.h>
#include "astro.h"
/* corrects ra and dec, both in radians, for precession from epoch 1 to epoch 2.
* the epochs are given by their modified JDs, mjd1 and mjd2, respectively.
* N.B. ra and dec are modifed IN PLACE.
*/
precess (mjd1, mjd2, ra, dec)
double mjd1, mjd2; /* initial and final epoch modified JDs */
double *ra, *dec; /* ra/dec for mjd1 in, for mjd2 out */
{
static double lastmjd1, lastmjd2;
static double m, n, nyrs;
double dra, ddec; /* ra and dec corrections */
if (mjd1 != lastmjd1 || mjd2 != lastmjd2) {
double t1, t2; /* Julian centuries of 36525 days since Jan .5 1900 */
double m1, n1; /* "constants" for t1 */
double m2, n2; /* "constants" for t2 */
t1 = mjd1/36525.;
t2 = mjd2/36525.;
m1 = 3.07234+(1.86e-3*t1);
n1 = 20.0468-(8.5e-3*t1);
m2 = 3.07234+(1.86e-3*t2);
n2 = 20.0468-(8.5e-3*t2);
m = (m1+m2)/2; /* average m for the two epochs */
n = (n1+n2)/2; /* average n for the two epochs */
nyrs = (mjd2-mjd1)/365.2425;
lastmjd1 = mjd1;
lastmjd2 = mjd2;
}
dra = (m+(n*sin(*ra)*tan(*dec)/15))*7.272205e-5*nyrs;
ddec = n*cos(*ra)*4.848137e-6*nyrs;
*ra += dra;
*dec += ddec;
range (ra, 2*PI);
}
xXx
echo x refract.c
cat > refract.c << 'xXx'
#include <stdio.h>
#include <math.h>
#include "astro.h"
/* correct the true altitude, ta, for refraction to the apparent altitude, aa,
* each in radians, given the local atmospheric pressure, pr, in mbars, and
* the temperature, tr, in degrees C.
*/
refract (pr, tr, ta, aa)
double pr, tr;
double ta;
double *aa;
{
double r; /* refraction correction*/
if (ta >= degrad(15.)) {
/* model for altitudes at least 15 degrees above horizon */
r = 7.888888e-5*pr/((273+tr)*tan(ta));
} else if (ta > degrad(-5.)) {
/* hairier model for altitudes at least -5 and below 15 degrees */
double a, b, tadeg = raddeg(ta);
a = ((2e-5*tadeg+1.96e-2)*tadeg+1.594e-1)*pr;
b = (273+tr)*((8.45e-2*tadeg+5.05e-1)*tadeg+1);
r = degrad(a/b);
} else {
/* do nothing if more than 5 degrees below horizon.
*/
r = 0;
}
*aa = ta + r;
}
/* correct the apparent altitude, aa, for refraction to the true altitude, ta,
* each in radians, given the local atmospheric pressure, pr, in mbars, and
* the temperature, tr, in degrees C.
*/
unrefract (pr, tr, aa, ta)
double pr, tr;
double aa;
double *ta;
{
double err;
double appar;
double true;
/* iterative solution: search for the true that refracts to the
* given apparent.
* since refract() is discontinuous at -5 degrees, there is a range
* of apparent altitudes between about -4.5 and -5 degrees that are
* not invertable (the graph of ap vs. true has a vertical step at
* true = -5). thus, the iteration just oscillates if it gets into
* this region. if this happens the iteration is forced to abort.
* of course, this makes unrefract() discontinuous too.
*/
true = aa;
do {
refract (pr, tr, true, &appar);
err = appar - aa;
true -= err;
} while (fabs(err) >= 1e-6 && true > degrad(-5));
*ta = true;
}
xXx
echo x riset.c
cat > riset.c << 'xXx'
#include <stdio.h>
#include <math.h>
#include "astro.h"
/* given the true geocentric ra and dec of an object, in radians, the observer's
* latitude in radians, and a horizon displacement correction, dis, in radians
* find the local sidereal times and azimuths of rising and setting, lstr/s
* and azr/s, in radians, respectively.
* dis is the vertical displacement from the true position of the horizon. it
* is positive if the apparent position is higher than the true position.
* said another way, it is positive if the shift causes the object to spend
* longer above the horizon. for example, atmospheric refraction is typically
* assumed to produce a vertical shift of 34 arc minutes at the horizon; dis
* would then take on the value +9.89e-3 (radians). On the other hand, if
* your horizon has hills such that your apparent horizon is, say, 1 degree
* above sea level, you would allow for this by setting dis to -1.75e-2
* (radians).
* status: 0: normal; 1: never rises; -1: circumpolar; 2: trouble.
*/
riset (ra, dec, lat, dis, lstr, lsts, azr, azs, status)
double ra, dec;
double lat, dis;
double *lstr, *lsts;
double *azr, *azs;
int *status;
{
static double lastlat = 0, slat = 0, clat = 1.0;
double dec1, sdec, cdec, tdec;
double psi, spsi, cpsi;
double h, dh, ch; /* hr angle, delta and cos */
/* avoid needless sin/cos since latitude doesn't change often */
if (lat != lastlat) {
clat = cos(lat);
slat = sin(lat);
lastlat = lat;
}
/* can't cope with objects very near the celestial pole nor if we
* are located very near the earth's poles.
*/
cdec = cos(dec);
if (fabs(cdec*clat) < 1e-10) {
/* trouble */
*status = 2;
return;
}
cpsi = slat/cdec;
if (cpsi>1.0) cpsi = 1.0;
else if (cpsi<-1.0) cpsi = -1.0;
psi = acos(cpsi);
spsi = sin(psi);
dh = dis*spsi;
dec1 = dec + (dis*cpsi);
sdec = sin(dec1);
tdec = tan(dec1);
ch = (-1*slat*tdec)/clat;
if (ch < -1) {
/* circumpolar */
*status = -1;
return;
}
if (ch > 1) {
/* never rises */
*status = 1;
return;
}
*status = 0;
h = acos(ch)+dh;
*lstr = 24+radhr(ra-h);
*lsts = radhr(ra+h);
range (lstr, 24.0);
range (lsts, 24.0);
*azr = acos(sdec/clat);
range (azr, 2*PI);
*azs = 2*PI - *azr;
}
xXx
echo x riset_c.c
cat > riset_c.c << 'xXx'
/* find rise and set circumstances, ie, riset_cir() and related functions. */
#include <stdio.h>
#include <math.h>
#include "astro.h"
#include "circum.h"
#include "screen.h" /* just for SUN and MOON */
#define TRACE(x) {FILE *fp = fopen("trace","a"); fprintf x; fclose(fp);}
#define STDREF degrad(34./60.) /* nominal horizon refraction amount */
#define TWIREF degrad(18.) /* twilight horizon displacement */
/* find where and when a body, p, will rise and set and
* its transit circumstances. all times are local, angles rads e of n.
* return 0 if just returned same stuff as previous call, else 1 if new.
* status is set from the RS_* #defines in circum.h.
* also used to find astro twilight by calling with dis of 18 degrees.
*/
riset_cir (p, np, hzn, ltr, lts, ltt, azr, azs, altt, status)
int p; /* one of the body defines in astro.h or screen.h */
Now *np;
int hzn; /* STDHZN or ADPHZN or TWILIGHT */
double *ltr, *lts; /* local rise and set times */
double *ltt; /* local transit time */
double *azr, *azs; /* local rise and set azimuths, rads e of n */
double *altt; /* local altitude at transit */
int *status; /* one or more of the RS_* defines */
{
typedef struct {
Now l_now;
double l_ltr, l_lts, l_ltt, l_azr, l_azs, l_altt;
int l_hzn;
int l_status;
} Last;
/* must be in same order as the astro.h/screen.h #define's */
static Last last[NOBJ];
Last *lp;
int new;
lp = last + p;
if (same_cir (np, &lp->l_now) && same_lday (np, &lp->l_now)
&& lp->l_hzn == hzn) {
*ltr = lp->l_ltr;
*lts = lp->l_lts;
*ltt = lp->l_ltt;
*azr = lp->l_azr;
*azs = lp->l_azs;
*altt = lp->l_altt;
*status = lp->l_status;
new = 0;
} else {
*status = 0;
iterative_riset (p, np, hzn, ltr, lts, ltt, azr, azs, altt, status);
lp->l_ltr = *ltr;
lp->l_lts = *lts;
lp->l_ltt = *ltt;
lp->l_azr = *azr;
lp->l_azs = *azs;
lp->l_altt = *altt;
lp->l_status = *status;
lp->l_hzn = hzn;
lp->l_now = *np;
new = 1;
}
return (new);
}
static
iterative_riset (p, np, hzn, ltr, lts, ltt, azr, azs, altt, status)
int p;
Now *np;
int hzn;
double *ltr, *lts, *ltt; /* local times of rise, set and transit */
double *azr, *azs, *altt;/* local azimuths of rise, set and transit altitude */
int *status;
{
#define MAXPASSES 6
double lstr, lsts, lstt; /* local sidereal times of rising/setting */
double mjd0; /* mjd estimates of rise/set event */
double lnoon; /* mjd of local noon */
double x; /* discarded tmp value */
Now n; /* just used to call now_lst() */
double lst; /* lst at local noon */
double diff, lastdiff; /* iterative improvement to mjd0 */
int pass;
int rss;
/* first approximation is to find rise/set times of a fixed object
* in its position at local noon.
*/
lnoon = mjd_day(mjd - tz/24.0) + (12.0 + tz)/24.0; /*mjd of local noon*/
n.n_mjd = lnoon;
n.n_lng = lng;
now_lst (&n, &lst); /* lst at local noon */
mjd0 = lnoon;
stationary_riset (p,mjd0,np,hzn,&lstr,&lsts,&lstt,&x,&x,&x,&rss);
chkrss:
switch (rss) {
case 0: break;
case 1: *status = RS_NEVERUP; return;
case -1: *status = RS_CIRCUMPOLAR; goto transit;
default: *status = RS_ERROR; return;
}
/* find a better approximation to the rising circumstances based on
* more passes, each using a "fixed" object at the location at
* previous approximation of the rise time.
*/
lastdiff = 1000.0;
for (pass = 1; pass < MAXPASSES; pass++) {
diff = (lstr - lst)*SIDRATE; /* next guess at rise time wrt noon */
if (diff > 12.0)
diff -= 24.0*SIDRATE; /* not tomorrow, today */
else if (diff < -12.0)
diff += 24.0*SIDRATE; /* not yesterday, today */
mjd0 = lnoon + diff/24.0; /* next guess at mjd of rise */
stationary_riset (p,mjd0,np,hzn,&lstr,&x,&x,azr,&x,&x,&rss);
if (rss != 0) goto chkrss;
if (fabs (diff - lastdiff) < 0.25/60.0)
break; /* converged to better than 15 secs */
lastdiff = diff;
}
if (pass == MAXPASSES)
*status |= RS_NORISE; /* didn't converge - no rise today */
else {
*ltr = 12.0 + diff;
if (p != MOON &&
(*ltr <= 24.0*(1.0-SIDRATE) || *ltr >= 24.0*SIDRATE))
*status |= RS_2RISES;
}
/* find a better approximation to the setting circumstances based on
* more passes, each using a "fixed" object at the location at
* previous approximation of the set time.
*/
lastdiff = 1000.0;
for (pass = 1; pass < MAXPASSES; pass++) {
diff = (lsts - lst)*SIDRATE; /* next guess at set time wrt noon */
if (diff > 12.0)
diff -= 24.0*SIDRATE; /* not tomorrow, today */
else if (diff < -12.0)
diff += 24.0*SIDRATE; /* not yesterday, today */
mjd0 = lnoon + diff/24.0; /* next guess at mjd of set */
stationary_riset (p,mjd0,np,hzn,&x,&lsts,&x,&x,azs,&x,&rss);
if (rss != 0) goto chkrss;
if (fabs (diff - lastdiff) < 0.25/60.0)
break; /* converged to better than 15 secs */
lastdiff = diff;
}
if (pass == MAXPASSES)
*status |= RS_NOSET; /* didn't converge - no set today */
else {
*lts = 12.0 + diff;
if (p != MOON &&
(*lts <= 24.0*(1.0-SIDRATE) || *lts >= 24.0*SIDRATE))
*status |= RS_2SETS;
}
transit:
/* find a better approximation to the transit circumstances based on
* more passes, each using a "fixed" object at the location at
* previous approximation of the transit time.
*/
lastdiff = 1000.0;
for (pass = 1; pass < MAXPASSES; pass++) {
diff = (lstt - lst)*SIDRATE; /*next guess at transit time wrt noon*/
if (diff > 12.0)
diff -= 24.0*SIDRATE; /* not tomorrow, today */
else if (diff < -12.0)
diff += 24.0*SIDRATE; /* not yesterday, today */
mjd0 = lnoon + diff/24.0; /* next guess at mjd of set */
stationary_riset (p,mjd0,np,hzn,&x,&x,&lstt,&x,&x,altt,&rss);
if (fabs (diff - lastdiff) < 0.25/60.0)
break; /* converged to better than 15 secs */
lastdiff = diff;
}
if (pass == MAXPASSES || rss != 0)
*status |= RS_NOTRANS; /* didn't converge - no transit today */
else {
*ltt = 12.0 + diff;
if (p != MOON &&
(*ltt <= 24.0*(1.0-SIDRATE) || *ltt >= 24.0*SIDRATE))
*status |= RS_2TRANS;
}
}
static
stationary_riset (p, mjd0, np, hzn, lstr, lsts, lstt, azr, azs, altt, status)
int p;
double mjd0;
Now *np;
int hzn;
double *lstr, *lsts, *lstt;
double *azr, *azs, *altt;
int *status;
{
double dis;
Now n;
Sky s;
switch (hzn) {
case STDHZN:
/* nominal atmospheric refraction.
* then add nominal moon or sun semi-diameter, as appropriate.
*/
dis = STDREF;
break;
case TWILIGHT:
if (p != SUN) {
c_erase();
printf ("BUG 1!\n\r");
exit(1);
}
dis = TWIREF;
break;
default:
/* horizon displacement is refraction plus object's semi-diameter */
unrefract (pressure, temp, 0.0, &dis);
dis = -dis;
n = *np;
n.n_mjd = mjd0;
(void) body_cir (p, 0.0, &n, &s);
dis += degrad(s.s_size/3600./2.0);
}
switch (p) {
case SUN:
if (hzn == STDHZN)
dis += degrad (32./60./2.);
fixedsunriset (mjd0,np,dis,lstr,lsts,lstt,azr,azs,altt,status);
break;
case MOON:
if (hzn == STDHZN)
dis += degrad (32./60./2.);
fixedmoonriset (mjd0,np,dis,lstr,lsts,lstt,azr,azs,altt,status);
break;
default:
if (hzn == STDHZN) {
/* assume planets have negligible diameters.
* also, need to set s if hzn was STDHZN.
*/
n = *np;
n.n_mjd = mjd0;
(void) body_cir (p, 0.0, &n, &s);
}
riset (s.s_ra, s.s_dec, lat, dis, lstr, lsts, azr, azs, status);
transit (s.s_ra, s.s_dec, np, lstt, altt);
}
}
/* find the local rise/set sidereal times and azimuths of an object fixed
* at the ra/dec of the sun on the given mjd time as seen from lat.
* times are for the upper limb. dis should account for refraction and
* sun's semi-diameter; we add corrections for nutation and aberration.
*/
static
fixedsunriset (mjd0, np, dis, lstr, lsts, lstt, azr, azs, altt, status)
double mjd0;
Now *np;
double dis;
double *lstr, *lsts, *lstt;
double *azr, *azs, *altt;
int *status;
{
double lsn, rsn;
double deps, dpsi;
double r, d;
/* find ecliptic position of sun at mjd */
sunpos (mjd0, &lsn, &rsn);
/* allow for 20.4" aberation and nutation */
nutation (mjd0, &deps, &dpsi);
lsn += dpsi - degrad(20.4/3600.0);
/* convert ecliptic to equatorial coords */
ecl_eq (mjd0, 0.0, lsn, &r, &d);
/* find circumstances for given horizon displacement */
riset (r, d, lat, dis, lstr, lsts, azr, azs, status);
transit (r, d, np, lstt, altt);
}
/* find the local sidereal rise/set times and azimuths of an object fixed
* at the ra/dec of the moon on the given mjd time as seen from lat.
* times are for the upper limb. dis should account for refraction and moon's
* semi-diameter; we add corrections for parallax, and nutation.
* accuracy is to nearest minute.
*/
static
fixedmoonriset (mjd0, np, dis, lstr, lsts, lstt, azr, azs, altt, status)
double mjd0;
Now *np;
double dis;
double *lstr, *lsts, *lstt;
double *azr, *azs, *altt;
int *status;
{
double lam, bet, hp;
double deps, dpsi;
double r, d;
double ha;
/* find geocentric ecliptic location and equatorial horizontal parallax
* of moon at mjd.
*/
moon (mjd0, &lam, &bet, &hp);
/* allow for nutation */
nutation (mjd0, &deps, &dpsi);
lam += dpsi;
/* convert ecliptic to equatorial coords */
ecl_eq (mjd0, bet, lam, &r, &d);
/* find local sidereal times of rise/set times/azimuths for given
* equatorial coords, allowing for
* .27249*sin(hp) parallax (hp is radius of earth from moon;
* .27249 is radius of moon from earth where
* the ratio is the dia_moon/dia_earth).
* hp nominal angular earth radius. subtract because
* tangential line-of-sight makes moon appear lower
* as move out from center of earth.
*/
dis += .27249*sin(hp) - hp;
riset (r, d, lat, dis, lstr, lsts, azr, azs, status);
/* account for parallax on the meridian (0 hour-angle).
* TODO: is this really right? other stuff too? better way?? help!
*/
ta_par (0.0, d, lat, height, hp, &ha, &d);
transit (r+ha, d, np, lstt, altt);
}
/* find when and how hi object at (r,d) is when it transits. */
static
transit (r, d, np, lstt, altt)
double r, d; /* ra and dec, rads */
Now *np; /* for refraction info */
double *lstt; /* local sidereal time of transit */
double *altt; /* local, refracted, altitude at time of transit */
{
*lstt = radhr(r);
*altt = PI/2 - lat + d;
if (*altt > PI/2)
*altt = PI - *altt;
refract (pressure, temp, *altt, altt);
}
xXx