allbery@uunet.UU.NET (Brandon S. Allbery - comp.sources.misc) (06/04/89)
Posting-number: Volume 7, Issue 11
Submitted-by: ecd@ncs-med.UUCP (Elwood C. Downey)
Archive-name: ephem/part03
[BTW, not putting "guard characters" at the beginnings of lines in a shar'ed
file is asking for trouble on some systems. Especially with idiot UUCP mail
systems -- some, perhaps all, of this group *does* get sent over mail links.
One "From" will produce a mangled transmission. And Bitnet seems to know of
other ways to destroy such files. ++bsa]
#!/bin/sh
echo mkdir lib
mkdir lib
echo cd lib
cd lib
echo extracting Makefile
cat > Makefile << 'xXx'
CLNFLAGS=$(CLNF)
LNFLAGS=$(CLNFLAGS) $(LNF)
CFLAGS=$(CLNFLAGS) $(CF)
LINTFLAGS=$(CFLAGS) $(LINTF)
LINTLIBS=
.c.a:
-$(CC) -c $(CFLAGS) $<
.PRECIOUS: lib.a
lib.a: lib.a(aa_hadec.o) \
lib.a(anomaly.o) \
lib.a(cal_mjd.o) \
lib.a(eq_ecl.o) \
lib.a(moon.o) \
lib.a(moonnf.o) \
lib.a(moonrs.o) \
lib.a(nutation.o) \
lib.a(obliq.o) \
lib.a(parallax.o) \
lib.a(pelement.o) \
lib.a(plans.o) \
lib.a(precess.o) \
lib.a(refract.o) \
lib.a(riset.o) \
lib.a(sex_dec.o) \
lib.a(sun.o) \
lib.a(sunrs.o) \
lib.a(utc_gst.o)
ar rv $@ $?
ranlib $@
rm -f $?
xXx
echo extracting aa_hadec.c
cat > aa_hadec.c << 'xXx'
#include <stdio.h>
#include <math.h>
#include "astro.h"
/* given latitude (n+, radians), lat, altitude (up+, radians), alt, and
* azimuth (angle round to the east from north+, radians),
* return hour angle (radians), ha, and declination (radians), dec.
*/
aa_hadec (lat, alt, az, ha, dec)
float lat;
float alt, az;
float *ha, *dec;
{
aaha_aux (lat, az, alt, ha, dec);
}
/* given latitude (n+, radians), lat, hour angle (radians), ha, and declination
* (radians), dec,
* return altitude (up+, radians), alt, and
* azimuth (angle round to the east from north+, radians),
*/
hadec_aa (lat, ha, dec, alt, az)
float lat;
float ha, dec;
float *alt, *az;
{
aaha_aux (lat, ha, dec, az, alt);
}
/* the actual formula is the same for both transformation directions so
* do it here once for each way.
* N.B. all arguments are in radians.
*/
static
aaha_aux (lat, x, y, p, q)
float lat;
float x, y;
float *p, *q;
{
static float lastlat = -1000.;
static float sinlastlat, coslastlat;
float sy, cy;
float sx, cx;
float sq, cq;
float a;
float cp;
/* latitude doesn't change much, so try to reuse the sin and cos evals.
*/
if (lat != lastlat) {
sinlastlat = sin (lat);
coslastlat = cos (lat);
lastlat = lat;
}
sy = sin (y);
cy = cos (y);
sx = sin (x);
cx = cos (x);
/* define GOODATAN2 if atan2 returns full range -PI through +PI.
*/
#ifdef GOODATAN2
*q = asin ((sy*sinlastlat) + (cy*coslastlat*cx));
*p = atan2 (-cy*sx, -cy*cx*sinlastlat + sy*coslastlat);
#else
#define EPS (1e-20)
sq = (sy*sinlastlat) + (cy*coslastlat*cx);
*q = asin (sq);
cq = cos (*q);
a = coslastlat*cq;
if (a > -EPS && a < EPS)
a = a < 0 ? -EPS : EPS; /* avoid / 0 */
cp = (sy - (sinlastlat*sq))/a;
if (cp >= 1.0) /* the /a can be slightly > 1 */
*p = 0.0;
else if (cp <= -1.0)
*p = PI;
else
*p = acos ((sy - (sinlastlat*sq))/a);
if (sx>0) *p = 2.0*PI - *p;
#endif
}
xXx
echo extracting anomaly.c
cat > anomaly.c << 'xXx'
#include <stdio.h>
#include <math.h>
#include "astro.h"
#define TWOPI (2*PI)
/* given the mean anomaly, ma, and the eccentricity, s, of elliptical motion,
* find the true anomaly, *nu, and the eccentric anomaly, *ea.
* all angles in radians.
*/
anomaly (ma, s, nu, ea)
float ma, s;
float *nu, *ea;
{
float m, dla, fea;
m = ma-TWOPI*(int)(ma/TWOPI);
fea = m;
while (1) {
dla = fea-(s*sin(fea))-m;
if (fabs(dla)<1e-6)
break;
dla /= 1-(s*cos(fea));
fea -= dla;
}
*nu = 2*atan(sqrt((1+s)/(1-s))*tan(fea/2));
*ea = fea;
}
xXx
echo extracting astro.h
cat > astro.h << 'xXx'
#ifndef PI
#define PI 3.141592653589793
#endif
/* conversions among hours (of ra), degrees and radians. */
#define degrad(x) ((x)*PI/180.)
#define raddeg(x) ((x)*180./PI)
#define hrdeg(x) ((x)*15.)
#define deghr(x) ((x)/15.)
#define hrrad(x) degrad(hrdeg(x))
#define radhr(x) deghr(raddeg(x))
/* manifest names for planets.
* N.B. must cooincide with usage in pelement.c and plans.c.
*/
#define MERCURY 0
#define VENUS 1
#define MARS 2
#define JUPITER 3
#define SATURN 4
#define URANUS 5
#define NEPTUNE 6
#define PLUTO 7
xXx
echo extracting cal_mjd.c
cat > cal_mjd.c << 'xXx'
#include <stdio.h>
#include <math.h>
#include "astro.h"
/* given a date in months, mn, days, dy, years, yr,
* return the modified Julian date (number of days elapsed since 1900 jan 0.5),
* *mjd.
*/
cal_mjd (mn, dy, yr, mjd)
int mn, yr;
float dy;
float *mjd;
{
int b, d, m, y;
long c;
m = mn;
y = (yr < 0) ? yr + 1 : yr;
if (mn < 3) {
m += 12;
y -= 1;
}
if (yr < 1582 || yr == 1582 && (mn < 10 || mn == 10 && dy < 15))
b = 0;
else {
int a;
a = y/100;
b = 2 - a + a/4;
}
if (y < 0)
c = (long)((365.25*y) - 0.75) - 694025L;
else
c = (long)(365.25*y) - 694025L;
d = 30.6001*(m+1);
*mjd = b + c + d + dy - 0.5;
}
/* given the modified Julian date (number of days elapsed since 1900 jan 0.5,),
* mjd, return the calendar date in months, *mn, days, *dy, and years, *yr.
*/
mjd_cal (mjd, mn, dy, yr)
float mjd;
int *mn, *yr;
float *dy;
{
float d, f;
float i, a, b, ce, g;
d = mjd + 0.5;
i = floor(d);
f = d-i;
if (f == 1) {
f = 0;
i += 1;
}
if (i > -115860.0) {
a = floor((i/36524.25)+.9983573)+14;
i += 1 + a - floor(a/4.0);
}
b = floor((i/365.25)+.802601);
ce = i - floor((365.25*b)+.750001)+416;
g = floor(ce/30.6001);
*mn = g - 1;
*dy = ce - floor(30.6001*g)+f;
*yr = b + 1899;
if (g > 13.5)
*mn = g - 13;
if (*mn < 2.5)
*yr = b + 1900;
if (*yr < 1)
*yr -= 1;
}
/* given an mjd, set *dow to 0..6 according to which dayof the week it falls
* on (0=sunday) or set it to -1 if can't figure it out.
*/
mjd_dow (mjd, dow)
float mjd;
int *dow;
{
/* cal_mjd() uses Gregorian dates on or after Oct 15, 1582.
* (Pope Gregory XIII dropped 10 days, Oct 5..14, and improved the leap-
* year algorithm). however, Great Britian and the colonies did not
* adopt it until Sept 14, 1752 (they dropped 11 days, Sept 3-13,
* due to additional accumulated error). leap years before 1752 thus
* can not easily be accounted for from the cal_mjd() number...
*/
if (mjd < -53798.5) {
/* pre sept 14, 1752 too hard to correct */
*dow = -1;
return;
}
*dow = ((int)floor(mjd-.5) + 1) % 7; /* 1/1/1900 (mjd 0.5) is a Monday*/
if (*dow < 0)
*dow += 7;
}
/* given a mjd, return the the number of days in the month. */
mjd_dpm (mjd, ndays)
float mjd;
int *ndays;
{
static short dpm[] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
int m, y;
float d;
mjd_cal (mjd, &m, &d, &y);
*ndays = (m==2 && ((y%4==0 && y%100!=0)||y%400==0)) ? 29 : dpm[m-1];
}
xXx
echo extracting eq_ecl.c
cat > eq_ecl.c << 'xXx'
#include <stdio.h>
#include <math.h>
#include "astro.h"
#define EQtoECL 1
#define ECLtoEQ (-1)
/* given the modified Julian date, mjd, and an equitorial ra and dec, each in
* radians, find the corresponding geocentric ecliptic latitude, *lat, and
* longititude, *lng, also each in radians.
* the effects of nutation and the angle of the obliquity are included.
*/
eq_ecl (mjd, ra, dec, lat, lng)
float mjd, ra, dec;
float *lat, *lng;
{
ecleq_aux (EQtoECL, mjd, ra, dec, lng, lat);
}
/* given the modified Julian date, mjd, and a geocentric ecliptic latitude,
* *lat, and longititude, *lng, each in radians, find the corresponding
* equitorial ra and dec, also each in radians.
* the effects of nutation and the angle of the obliquity are included.
*/
ecl_eq (mjd, lat, lng, ra, dec)
float mjd, lat, lng;
float *ra, *dec;
{
ecleq_aux (ECLtoEQ, mjd, lng, lat, ra, dec);
}
static
ecleq_aux (sw, mjd, x, y, p, q)
int sw; /* +1 for eq to ecliptic, -1 for vv. */
float mjd, x, y; /* sw==1: x==ra, y==dec. sw==-1: x==lng, y==lat. */
float *p, *q; /* sw==1: p==lng, q==lat. sw==-1: p==ra, q==dec. */
{
static float lastmjd; /* last mjd calculated */
static float seps, ceps; /* sin and cos of mean obliquity */
float sx, cx, sy, cy, ty;
if (mjd != lastmjd) {
float deps, dpsi, eps;
obliquity (mjd, &eps); /* mean obliquity for date */
#ifdef NONUTATION
#define NONUTATION
deps = 0; /* checkout handler did not
* include nutation correction.
*/
#else
nutation (mjd, &deps, &dpsi); /* correctin due to nutation */
#endif
eps += deps; /* true obliquity for date */
seps = sin(eps);
ceps = cos(eps);
lastmjd = mjd;
}
sy = sin(y);
cy = cos(y); /* always non-negative */
if (fabs(cy)<1e-20) cy = 1e-20; /* insure > 0 */
ty = sy/cy;
cx = cos(x);
sx = sin(x);
*q = asin((sy*ceps)-(cy*seps*sx*sw));
*p = atan(((sx*ceps)+(ty*seps*sw))/cx);
if (cx<0) *p += PI; /* account for atan quad ambiguity */
range (p, 2*PI);
}
xXx
echo extracting moon.c
cat > moon.c << 'xXx'
#include <stdio.h>
#include <math.h>
#include "astro.h"
/* given the mjd, find the geocentric ecliptic longitude, lam, and latitude,
* bet, and horizontal parallax, hp for the moon.
* N.B. series for long and lat are good to about 10 and 3 arcseconds. however,
* math errors cause up to 100 and 30 arcseconds error, even if use double.
* why?? suspect highly sensitive nature of difference used to get m1..6.
* N.B. still need to correct for nutation. then for topocentric location
* further correct for parallax and refraction.
*/
moon (mjd, lam, bet, hp)
float mjd;
float *lam, *bet, *hp;
{
float t, t2;
float ld;
float ms;
float md;
float de;
float f;
float n;
float a, sa, sn, b, sb, c, sc, e, e2, l, g, w1, w2;
float m1, m2, m3, m4, m5, m6;
t = mjd/36525.;
t2 = t*t;
m1 = mjd/27.32158213;
m1 = 360.0*(m1-(int)m1);
m2 = mjd/365.2596407;
m2 = 360.0*(m2-(int)m2);
m3 = mjd/27.55455094;
m3 = 360.0*(m3-(int)m3);
m4 = mjd/29.53058868;
m4 = 360.0*(m4-(int)m4);
m5 = mjd/27.21222039;
m5 = 360.0*(m5-(int)m5);
m6 = mjd/6798.363307;
m6 = 360.0*(m6-(int)m6);
ld = 270.434164+m1-(.001133-.0000019*t)*t2;
ms = 358.475833+m2-(.00015+.0000033*t)*t2;
md = 296.104608+m3+(.009192+.0000144*t)*t2;
de = 350.737486+m4-(.001436-.0000019*t)*t2;
f = 11.250889+m5-(.003211+.0000003*t)*t2;
n = 259.183275-m6+(.002078+.000022*t)*t2;
a = degrad(51.2+20.2*t);
sa = sin(a);
sn = sin(degrad(n));
b = 346.56+(132.87-.0091731*t)*t;
sb = .003964*sin(degrad(b));
c = degrad(n+275.05-2.3*t);
sc = sin(c);
ld = ld+.000233*sa+sb+.001964*sn;
ms = ms-.001778*sa;
md = md+.000817*sa+sb+.002541*sn;
f = f+sb-.024691*sn-.004328*sc;
de = de+.002011*sa+sb+.001964*sn;
e = 1-(.002495+7.52e-06*t)*t;
e2 = e*e;
ld = degrad(ld);
ms = degrad(ms);
n = degrad(n);
de = degrad(de);
f = degrad(f);
md = degrad(md);
l = 6.28875*sin(md)+1.27402*sin(2*de-md)+.658309*sin(2*de)+
.213616*sin(2*md)-e*.185596*sin(ms)-.114336*sin(2*f)+
.058793*sin(2*(de-md))+.057212*e*sin(2*de-ms-md)+
.05332*sin(2*de+md)+.045874*e*sin(2*de-ms)+.041024*e*sin(md-ms);
l = l-.034718*sin(de)-e*.030465*sin(ms+md)+.015326*sin(2*(de-f))-
.012528*sin(2*f+md)-.01098*sin(2*f-md)+.010674*sin(4*de-md)+
.010034*sin(3*md)+.008548*sin(4*de-2*md)-e*.00791*sin(ms-md+2*de)-
e*.006783*sin(2*de+ms);
l = l+.005162*sin(md-de)+e*.005*sin(ms+de)+.003862*sin(4*de)+
e*.004049*sin(md-ms+2*de)+.003996*sin(2*(md+de))+
.003665*sin(2*de-3*md)+e*.002695*sin(2*md-ms)+
.002602*sin(md-2*(f+de))+e*.002396*sin(2*(de-md)-ms)-
.002349*sin(md+de);
l = l+e2*.002249*sin(2*(de-ms))-e*.002125*sin(2*md+ms)-
e2*.002079*sin(2*ms)+e2*.002059*sin(2*(de-ms)-md)-
.001773*sin(md+2*(de-f))-.001595*sin(2*(f+de))+
e*.00122*sin(4*de-ms-md)-.00111*sin(2*(md+f))+.000892*sin(md-3*de);
l = l-e*.000811*sin(ms+md+2*de)+e*.000761*sin(4*de-ms-2*md)+
e2*.000704*sin(md-2*(ms+de))+e*.000693*sin(ms-2*(md-de))+
e*.000598*sin(2*(de-f)-ms)+.00055*sin(md+4*de)+.000538*sin(4*md)+
e*.000521*sin(4*de-ms)+.000486*sin(2*md-de);
l = l+e2*.000717*sin(md-2*ms);
*lam = ld+degrad(l);
range (lam, 2*PI);
g = 5.12819*sin(f)+.280606*sin(md+f)+.277693*sin(md-f)+
.173238*sin(2*de-f)+.055413*sin(2*de+f-md)+.046272*sin(2*de-f-md)+
.032573*sin(2*de+f)+.017198*sin(2*md+f)+.009267*sin(2*de+md-f)+
.008823*sin(2*md-f)+e*.008247*sin(2*de-ms-f);
g = g+.004323*sin(2*(de-md)-f)+.0042*sin(2*de+f+md)+
e*.003372*sin(f-ms-2*de)+e*.002472*sin(2*de+f-ms-md)+
e*.002222*sin(2*de+f-ms)+e*.002072*sin(2*de-f-ms-md)+
e*.001877*sin(f-ms+md)+.001828*sin(4*de-f-md)-e*.001803*sin(f+ms)-
.00175*sin(3*f);
g = g+e*.00157*sin(md-ms-f)-.001487*sin(f+de)-e*.001481*sin(f+ms+md)+
e*.001417*sin(f-ms-md)+e*.00135*sin(f-ms)+.00133*sin(f-de)+
.001106*sin(f+3*md)+.00102*sin(4*de-f)+.000833*sin(f+4*de-md)+
.000781*sin(md-3*f)+.00067*sin(f+4*de-2*md);
g = g+.000606*sin(2*de-3*f)+.000597*sin(2*(de+md)-f)+
e*.000492*sin(2*de+md-ms-f)+.00045*sin(2*(md-de)-f)+
.000439*sin(3*md-f)+.000423*sin(f+2*(de+md))+
.000422*sin(2*de-f-3*md)-e*.000367*sin(ms+f+2*de-md)-
e*.000353*sin(ms+f+2*de)+.000331*sin(f+4*de);
g = g+e*.000317*sin(2*de+f-ms+md)+e2*.000306*sin(2*(de-ms)-f)-
.000283*sin(md+3*f);
w1 = .0004664*cos(n);
w2 = .0000754*cos(c);
*bet = degrad(g)*(1-w1-w2);
*hp = .950724+.051818*cos(md)+.009531*cos(2*de-md)+.007843*cos(2*de)+
.002824*cos(2*md)+.000857*cos(2*de+md)+e*.000533*cos(2*de-ms)+
e*.000401*cos(2*de-md-ms)+e*.00032*cos(md-ms)-.000271*cos(de)-
e*.000264*cos(ms+md)-.000198*cos(2*f-md);
*hp = *hp+.000173*cos(3*md)+.000167*cos(4*de-md)-e*.000111*cos(ms)+
.000103*cos(4*de-2*md)-.000084*cos(2*md-2*de)-
e*.000083*cos(2*de+ms)+.000079*cos(2*de+2*md)+.000072*cos(4*de)+
e*.000064*cos(2*de-ms+md)-e*.000063*cos(2*de+ms-md)+
e*.000041*cos(ms+de);
*hp = *hp+e*.000035*cos(2*md-ms)-.000033*cos(3*md-2*de)-
.00003*cos(md+de)-.000029*cos(2*(f-de))-e*.000029*cos(2*md+ms)+
e2*.000026*cos(2*(de-ms))-.000023*cos(2*(f-de)+md)+
e*.000019*cos(4*de-ms-md);
*hp = degrad(*hp);
}
xXx
echo extracting moonnf.c
cat > moonnf.c << 'xXx'
#include <stdio.h>
#include <math.h>
#include "astro.h"
#define unw(w,z) ((w)-floor((w)/(z))*(z))
/* given a modified Julian date, mjd, return the mjd of the new
* and full moons about then, mjdn and mjdf.
* TODO: exactly which ones does it find? eg:
* 5/28/1988 yields 5/15 and 5/31
* 5/29 6/14 6/29
*/
moonnf (mjd, mjdn, mjdf)
float mjd;
float *mjdn, *mjdf;
{
int mo, yr;
float dy;
float mjd0;
float k, tn, tf, t;
mjd_cal (mjd, &mo, &dy, &yr);
cal_mjd (1, 0., yr, &mjd0);
k = (yr-1900+((mjd-mjd0)/365))*12.3685;
k = floor(k+0.5);
tn = k/1236.85;
tf = (k+0.5)/1236.85;
t = tn;
m (t, k, mjdn);
t = tf;
k += 0.5;
m (t, k, mjdf);
}
static
m (t, k, mjd)
float t, k;
float *mjd;
{
float t2, a, a1, b, b1, c, ms, mm, f, ddjd;
t2 = t*t;
a = 29.53*k;
c = degrad(166.56+(132.87-9.173e-3*t)*t);
b = 5.8868e-4*k+(1.178e-4-1.55e-7*t)*t2+3.3e-4*sin(c)+7.5933E-1;
ms = 359.2242+360*unw(k/1.236886e1,1)-(3.33e-5+3.47e-6*t)*t2;
mm = 306.0253+360*unw(k/9.330851e-1,1)+(1.07306e-2+1.236e-5*t)*t2;
f = 21.2964+360*unw(k/9.214926e-1,1)-(1.6528e-3+2.39e-6*t)*t2;
ms = unw(ms,360);
mm = unw(mm,360);
f = unw(f,360);
ms = degrad(ms);
mm = degrad(mm);
f = degrad(f);
ddjd = (1.734e-1-3.93e-4*t)*sin(ms)+2.1e-3*sin(2*ms)
-4.068e-1*sin(mm)+1.61e-2*sin(2*mm)-4e-4*sin(3*mm)
+1.04e-2*sin(2*f)-5.1e-3*sin(ms+mm)-7.4e-3*sin(ms-mm)
+4e-4*sin(2*f+ms)-4e-4*sin(2*f-ms)-6e-4*sin(2*f+mm)
+1e-3*sin(2*f-mm)+5e-4*sin(ms+2*mm);
a1 = (int)a;
b = b+ddjd+(a-a1);
b1 = (int)b;
a = a1+b1;
b = b-b1;
*mjd = a + b;
}
xXx
echo extracting moonrs.c
cat > moonrs.c << 'xXx'
#include <stdio.h>
#include <math.h>
#include "astro.h"
#define PASSES 2
/* given the mjd and geographical latitude and longitude, lat and lng,
* find the ut and azimuth of moonrise, utcr and azr, and moonset, utcs and
* azs for the day of mjd. all angles are in radians.
* times are for the upper limb, including effects of refraction (a constant 34
* arcminutes), parallax, and nutation. accuracy is to nearest minute.
* the calculated times can be in error if they occur within about a half hour
* of utc midnight due to the algorithm's confusion about the date. the best
* way to verify it is to calculate the times for several days either side of
* the date to see whether the problem values conform to a smooth trend set
* by the others.
* possible values of status:
* 2: can't cope (such as geographical latitude very near +-90).
* 1: moon never rises/sets on date.
* 0: all ok.
* -1: moon is circumpolar.
* -2: possible error in near-midnight moonrise time; see above.
* -3: possible error in near-midnight moonset time; see above.
*/
moonrs (mjd, lat, lng, utcr, utcs, azr, azs, status)
float mjd, lat, lng;
float *utcr, *utcs;
float *azr, *azs;
int *status;
{
float lstr, lsts; /* local sidereal times of rising/setting */
float al = radhr(lng); /* time correction for longitude */
int midnight = 0; /* midnight caution */
float mjd0; /* start of mjd day */
float x; /* discarded tmp value */
float ut;
int pass;
mjd0 = floor(mjd+0.5)-0.5;
/* first approximation is to find rise/set times of a fixed object
* in the position of the moon at local midday.
*/
fixedmoonriset (mjd0+(12.0-al)/24.0, lat, &lstr, &lsts, &x, &x, status);
if (*status != 0) return;
/* find a better approximation to the rising circumstances based on two
* more passes, each using a "fixed" object at the moons location at
* previous approximation of the rise time.
*/
pass = 0;
while (1) {
gst_utc (mjd0, lstr, &ut);
if (++pass > PASSES)
break;
fixedmoonriset (mjd0+(ut-al)/24., lat, &lstr, &x, azr, &x, status);
if (*status != 0) return;
}
if (ut > 23.5 || ut < 0.5)
midnight = -2; /* moonrise caution near midnight */
*utcr = ut - (al * .9972696); /* allow for sidereal shift */
/* find a better approximation to the setting circumstances based on two
* more passes, each using a "fixed" object at the moons location at
* previous approximation of the setting time.
*/
pass = 0;
while (1) {
gst_utc (mjd0, lsts, &ut);
if (++pass > PASSES)
break;
fixedmoonriset (mjd0+(ut-al)/24., lat, &x, &lsts, &x, azs, status);
if (*status != 0) return;
}
if (ut > 23.5 || ut < 0.5)
midnight = -3; /* moonset caution near midnight */
*utcs = ut - (al * .9972696); /* allow for sidereal shift */
if (midnight)
*status = midnight; /* report caution if near midnight */
}
/* find the local rise/set sidereal times and azimuths of an object fixed
* at the ra/dec of the moon on the given mjd time as seen from lat.
*/
static
fixedmoonriset (mjd, lat, lstr, lsts, azr, azs, status)
float mjd, lat;
float *lstr, *lsts;
float *azr, *azs;
int *status;
{
float lam, bet, hp;
float deps, dpsi;
float dis;
float ra, dec;
/* find geocentric ecliptic location and equatorial horizontal parallax
* of moon at mjd.
*/
moon (mjd, &lam, &bet, &hp);
/* allow for nutation */
nutation (mjd, &deps, &dpsi);
lam += dpsi;
/* convert ecliptic to equatorial coords */
ecl_eq (mjd, bet, lam, &ra, &dec);
/* 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).
* 34' (9.8902e-3) nominal refraction
* hp nominal angular earth radius. subtract because
* tangential line-of-sight makes moon appear lower
* as move out from center of earth.
* TODO: do better at refraction; see fixedsunriset() in sunrs.c.
*/
dis = .27249*sin(hp) + 9.8902e-3 - hp;
riset (ra, dec, lat, dis, lstr, lsts, azr, azs, status);
}
xXx
echo extracting nutation.c
cat > nutation.c << 'xXx'
#include <stdio.h>
#include <math.h>
#include "astro.h"
/* given the modified JD, mjd, find the nutation in obliquity, *deps, and
* the nutation in longitude, *dpsi, each in radians.
*/
nutation (mjd, deps, dpsi)
float mjd;
float *deps, *dpsi;
{
static float lastmjd, lastdeps, lastdpsi;
float ls, ld; /* sun's mean longitude, moon's mean longitude */
float ms, md; /* sun's mean anomaly, moon's mean anomaly */
float nm; /* longitude of moon's ascending node */
float t, t2; /* number of Julian centuries of 36525 days since
* Jan 0.5 1900.
*/
float tls, tnm, tld; /* twice above */
double a, b; /* temps */
if (mjd == lastmjd) {
*deps = lastdeps;
*dpsi = lastdpsi;
return;
}
t = mjd/36525.;
t2 = t*t;
a = 100.0021358*t;
b = 360.*(a-(int)a);
ls = 279.697+.000303*t2+b;
a = 1336.855231*t;
b = 360.*(a-(int)a);
ld = 270.434-.001133*t2+b;
a = 99.99736056000026*t;
b = 360.*(a-(int)a);
ms = 358.476-.00015*t2+b;
a = 13255523.59*t;
b = 360.*(a-(int)a);
md = 296.105+.009192*t2+b;
a = 5.372616667*t;
b = 360.*(a-(int)a);
nm = 259.183+.002078*t2-b;
/* convert to radian forms for use with trig functions.
*/
tls = 2*degrad(ls);
nm = degrad(nm);
tnm = 2*degrad(nm);
ms = degrad(ms);
tld = 2*degrad(ld);
md = degrad(md);
/* find delta psi and eps, in arcseconds.
*/
lastdpsi = (-17.2327-.01737*t)*sin(nm)+(-1.2729-.00013*t)*sin(tls)
+.2088*sin(tnm)-.2037*sin(tld)+(.1261-.00031*t)*sin(ms)
+.0675*sin(md)-(.0497-.00012*t)*sin(tls+ms)
-.0342*sin(tld-nm)-.0261*sin(tld+md)+.0214*sin(tls-ms)
-.0149*sin(tls-tld+md)+.0124*sin(tls-nm)+.0114*sin(tld-md);
lastdeps = (9.21+.00091*t)*cos(nm)+(.5522-.00029*t)*cos(tls)
-.0904*cos(tnm)+.0884*cos(tld)+.0216*cos(tls+ms)
+.0183*cos(tld-nm)+.0113*cos(tld+md)-.0093*cos(tls-ms)
-.0066*cos(tls-nm);
/* convert to radians.
*/
lastdpsi = degrad(lastdpsi/3600);
lastdeps = degrad(lastdeps/3600);
lastmjd = mjd;
*deps = lastdeps;
*dpsi = lastdpsi;
}
xXx
echo extracting obliq.c
cat > obliq.c << 'xXx'
#include <stdio.h>
#include "astro.h"
/* given the modified Julian date, mjd, find the obliquity of the
* ecliptic, *eps, in radians.
*/
obliquity (mjd, eps)
float mjd;
float *eps;
{
static float lastmjd, lasteps;
if (mjd != lastmjd) {
float t;
t = mjd/36525.;
lasteps = degrad(2.345229444E1 - ((((-1.81E-3*t)+5.9E-3)*t+4.6845E1)*t)/3600.0);
lastmjd = mjd;
}
*eps = lasteps;
}
xXx
echo extracting parallax.c
cat > parallax.c << 'xXx'
#include <stdio.h>
#include <math.h>
#include "astro.h"
/* given true ha and dec, tha and tdec, the geographical latitude, phi, the
* height above sea-level (as a fraction of the earths radius, 6378.16km),
* ht, and the equatorial horizontal parallax, ehp, find the apparent
* ha and dec, aha and adec allowing for parallax.
* all angles in radians. ehp is the angle subtended at the body by the
* earth's equator.
*/
ta_par (tha, tdec, phi, ht, ehp, aha, adec)
float tha, tdec, phi, ht, ehp;
float *aha, *adec;
{
static float last_phi, last_ht, rsp, rcp;
float rp; /* distance to object in Earth radii */
float ctha;
float stdec, ctdec;
float tdtha, dtha;
float caha;
/* avoid calcs involving the same phi and ht */
if (phi != last_phi || ht != last_ht) {
float cphi, sphi, u;
cphi = cos(phi);
sphi = sin(phi);
u = atan(9.96647e-1*sphi/cphi);
rsp = (9.96647e-1*sin(u))+(ht*sphi);
rcp = cos(u)+(ht*cphi);
last_phi = phi;
last_ht = ht;
}
rp = 1/sin(ehp);
ctha = cos(tha);
stdec = sin(tdec);
ctdec = cos(tdec);
tdtha = (rcp*sin(tha))/((rp*ctdec)-(rcp*ctha));
dtha = atan(tdtha);
*aha = tha+dtha;
caha = cos(*aha);
range (aha, 2*PI);
*adec = atan(caha*(rp*stdec-rsp)/(rp*ctdec*ctha-rcp));
}
/* given the apparent ha and dec, aha and adec, the geographical latitude, phi,
* the height above sea-level (as a fraction of the earths radius, 6378.16km),
* ht, and the equatorial horizontal parallax, ehp, find the true ha and dec,
* tha and tdec allowing for parallax.
* all angles in radians. ehp is the angle subtended at the body by the
* earth's equator.
* uses ta_par() iteratively: find a set of true ha/dec that converts back
* to the given apparent ha/dec.
*/
at_par (aha, adec, phi, ht, ehp, tha, tdec)
float aha, adec, phi, ht, ehp;
float *tha, *tdec;
{
float nha, ndec; /* ha/dec corres. to current true guesses */
float eha, edec; /* error in ha/dec */
/* first guess for true is just the apparent */
*tha = aha;
*tdec = adec;
while (1) {
ta_par (*tha, *tdec, phi, ht, ehp, &nha, &ndec);
eha = aha - nha;
edec = adec - ndec;
if (fabs(eha)<1e-6 && fabs(edec)<1e-6)
break;
*tha += eha;
*tdec += edec;
}
}
xXx
echo extracting pelement.c
cat > pelement.c << 'xXx'
#include <stdio.h>
#include <math.h>
#include "astro.h"
/* this array contains polynomial coefficients to find the various orbital
* elements for the mean orbit at any instant in time for each major planet.
* the first five elements are in the form a0 + a1*t + a2*t**2 + a3*t**3,
* where t is the number of Julian centuries of 36525 Julian days since 1900
* Jan 0.5. the last three elements are constants.
*
* the orbital element (column) indeces are:
* [ 0- 3]: coefficients for mean longitude, in degrees;
* [ 4- 7]: coefficients for longitude of the perihelion, in degrees;
* [ 8-11]: coefficients for eccentricity;
* [12-15]: coefficients for inclination, in degrees;
* [16-19]: coefficients for longitude of the ascending node, in degrees;
* [20]: semi-major axis, in AU;
* [21]: angular diameter at 1 AU, in arcsec;
* [22]: standard visual magnitude, ie, the visual magnitude of the planet
* when at a distance of 1 AU from both the Sun and the Earth and
* with zero phase angle.
*
* the planent (row) indeces are:
* [0]: Mercury; [1]: Venus; [2]: Mars; [3]: Jupiter; [4]: Saturn;
* [5]: Uranus; [6]: Neptune; [7]: Pluto.
*/
#define NPELE (5*4 + 3) /* 4 coeffs for ea of 5 elems, + 3 constants */
static float elements[8][NPELE] = {
{ /* mercury... */
178.179078, 415.2057519, 3.011e-4, 0.0,
75.899697, 1.5554889, 2.947e-4, 0.0,
.20561421, 2.046e-5, 3e-8, 0.0,
7.002881, 1.8608e-3, -1.83e-5, 0.0,
47.145944, 1.1852083, 1.739e-4, 0.0,
.3870986, 6.74, -0.42
},
{ /* venus... */
342.767053, 162.5533664, 3.097e-4, 0.0,
130.163833, 1.4080361, -9.764e-4, 0.0,
6.82069e-3, -4.774e-5, 9.1e-8, 0.0,
3.393631, 1.0058e-3, -1e-6, 0.0,
75.779647, .89985, 4.1e-4, 0.0,
.7233316, 16.92, -4.4
},
{ /* mars... */
293.737334, 53.17137642, 3.107e-4, 0.0,
3.34218203e2, 1.8407584, 1.299e-4, -1.19e-6,
9.33129e-2, 9.2064e-5, 7.7e-8, 0.0,
1.850333, -6.75e-4, 1.26e-5, 0.0,
48.786442, .7709917, -1.4e-6, -5.33e-6,
1.5236883, 9.36, -1.52
},
{ /* jupiter... */
238.049257, 8.434172183, 3.347e-4, -1.65e-6,
1.2720972e1, 1.6099617, 1.05627e-3, -3.43e-6,
4.833475e-2, 1.6418e-4, -4.676e-7, -1.7e-9,
1.308736, -5.6961e-3, 3.9e-6, 0.0,
99.443414, 1.01053, 3.5222e-4, -8.51e-6,
5.202561, 196.74, -9.4
},
{ /* saturn... */
266.564377, 3.398638567, 3.245e-4, -5.8e-6,
9.1098214e1, 1.9584158, 8.2636e-4, 4.61e-6,
5.589232e-2, -3.455e-4, -7.28e-7, 7.4e-10,
2.492519, -3.9189e-3, -1.549e-5, 4e-8,
112.790414, .8731951, -1.5218e-4, -5.31e-6,
9.554747, 165.6, -8.88
},
{ /* uranus... */
244.19747, 1.194065406, 3.16e-4, -6e-7,
1.71548692e2, 1.4844328, 2.372e-4, -6.1e-7,
4.63444e-2, -2.658e-5, 7.7e-8, 0.0,
.772464, 6.253e-4, 3.95e-5, 0.0,
73.477111, .4986678, 1.3117e-3, 0.0,
19.21814, 65.8, -7.19
},
{ /* neptune... */
84.457994, .6107942056, 3.205e-4, -6e-7,
4.6727364e1, 1.4245744, 3.9082e-4, -6.05e-7,
8.99704e-3, 6.33e-6, -2e-9, 0.0,
1.779242, -9.5436e-3, -9.1e-6, 0.0,
130.681389, 1.098935, 2.4987e-4, -4.718e-6,
30.10957, 62.2, -6.87
},
{ /* pluto...(osculating 1984 jan 21) */
95.3113544, .3980332167, 0.0, 0.0,
224.017, 0.0, 0.0, 0.0,
.25515, 0.0, 0.0, 0.0,
17.1329, 0.0, 0.0, 0.0,
110.191, 0.0, 0.0, 0.0,
39.8151, 8.2, -1.0
}
};
/* given a modified Julian date, mjd, return the elements for the mean orbit
* at that instant of all the major planets, together with their
* mean daily motions in longitude, angular diameter and standard visual
* magnitude.
* plan[i][j] contains all the values for all the planets at mjd, such that
* i = 0..7: mercury, venus, mars, jupiter, saturn, unranus, neptune, pluto;
* j = 0..8: mean longitude, mean daily motion in longitude, longitude of
* the perihelion, eccentricity, inclination, longitude of the ascending
* node, length of the semi-major axis, angular diameter from 1 AU, and
* the standard visual magnitude (see elements[][] comment, above).
*/
pelement (mjd, plan)
float mjd;
float plan[8][9];
{
register float *ep, *pp;
register float t = mjd/36525.;
float aa;
int planet, i;
for (planet = 0; planet < 8; planet++) {
ep = elements[planet];
pp = plan[planet];
aa = ep[1]*t;
pp[0] = ep[0] + 360.*(aa-(int)aa) + (ep[3]*t + ep[2])*t*t;
range (pp, 360.);
pp[1] = (ep[1]*9.856263e-3) + (ep[2] + ep[3])/36525;
for (i = 4; i < 20; i += 4)
pp[i/4+1] = ((ep[i+3]*t + ep[i+2])*t + ep[i+1])*t + ep[i+0];
pp[6] = ep[20];
pp[7] = ep[21];
pp[8] = ep[22];
}
}
xXx
echo extracting 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)
float mjd;
int p;
float *lpd0, *psi0, *rp0, *rho0, *lam, *bet, *dia, *mag;
{
static float plan[8][9];
static lastmjd = -10000;
float dl; /* perturbation correction for longitude */
float dr; /* " orbital radius */
float dml; /* " mean longitude */
float ds; /* " eccentricity */
float dm; /* " mean anomaly */
float da; /* " semi-major axis */
float dhl; /* " heliocentric longitude */
float lsn, rsn; /* true geocentric longitude of sun and sun-earth rad */
float mas; /* mean anomaly of the sun */
float re; /* radius of earth's orbit */
float lg; /* longitude of earth */
float map[8]; /* array of mean anomalies for each planet */
float lpd, psi, rp, rho;
float ll, sll, cll;
float t;
float dt;
int pass;
int j;
float s, ma;
float nu, ea;
float lp, om;
float lo, slo, clo;
float inc, y;
float spsi, cpsi;
float 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.;
sun (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)
float t;
float *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)
float mjd;
float *mas;
{
float t, t2;
float 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)
float map[];
float *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)
float t, mas, map[];
float *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)
float mas, map[];
float *dl, *dr, *dml, *dm;
{
float 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)+
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)
float t, s;
float *dml, *ds, *dm, *da;
{
float dp;
float j1, j2, j3, j4, j5, j6, j7;
float sj3, cj3, s2j3, c2j3;
float sj5, cj5, s2j5;
float sj6;
float 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-
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-
(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)
float t, s;
float *dml, *ds, *dm, *da, *dhl;
{
float dp;
float j1, j2, j3, j4, j5, j6, j7, j8;
float sj3, cj3, s2j3, c2j3, s3j3, c3j3, s4j3, c4j3;
float sj5, cj5, s2j5, c2j5;
float sj6;
float sj7, cj7, s2j7, c2j7, s3j7, c3j7, s4j7, c4j7, c5j7, s5j7;
float 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+
(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-
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+
(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-
(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+
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)
float t, s;
float *dl, *dr, *dml, *ds, *dm, *da, *dhl;
{
float dp;
float j1, j2, j3, j4, j5, j6;
float j8, j9, j10, j11, j12;
float sj4, cj4, s2j4, c2j4;
float sj9, cj9, s2j9, c2j9;
float 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)
float t, s;
float *dl, *dr, *dml, *ds, *dm, *da, *dhl;
{
float dp;
float j1, j2, j3, j4, j5, j6;
float j8, j9, j10, j11, j12;
float sj8, cj8;
float sj9, cj9, s2j9, c2j9;
float 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 extracting 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)
float mjd1, mjd2; /* initial and final epoch modified JDs */
float *ra, *dec; /* ra/dec for mjd1 in, for mjd2 out */
{
static float lastmjd1, lastmjd2;
static float m, n, nyrs;
float dra, ddec; /* ra and dec corrections */
if (mjd1 != lastmjd1 || mjd2 != lastmjd2) {
float t1, t2; /* Julian centuries of 36525 days since Jan .5 1900 */
float m1, n1; /* "constants" for t1 */
float 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 extracting 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)
float pr, tr;
float ta;
float *aa;
{
float 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 */
float 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)
float pr, tr;
float aa;
float *ta;
{
float err;
float appar;
float 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 extracting 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)
float ra, dec;
float lat, dis;
float *lstr, *lsts;
float *azr, *azs;
int *status;
{
static float lastlat = 0, slat = 0, clat = 1.0;
float dec1, sdec, cdec, tdec;
float psi, spsi, cpsi;
float 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 extracting sex_dec.c
cat > sex_dec.c << 'xXx'
/* given hours (or degrees), hd, minutes, m, and seconds, s,
* return decimal hours (or degrees), *d.
* in the case of hours (angles) < 0, only the first non-zero element should
* be negative.
*/
sex_dec (hd, m, s, d)
int hd, m, s;
float *d;
{
int sign = 1;
if (hd < 0) {
sign = -1;
hd = -hd;
} else if (m < 0) {
sign = -1;
m = -m;
} else if (s < 0) {
sign = -1;
s = -s;
}
*d = ((s/60.0 + m)/60.0 + hd) * sign;
}
/* given decimal hours (or degrees), d.
* return nearest hours (or degrees), *hd, minutes, *m, and seconds, *s,
* each always non-negative; *isneg is set to 1 if d is < 0, else to 0.
*/
dec_sex (d, hd, m, s, isneg)
float d;
int *hd, *m, *s, *isneg;
{
float min;
if (d < 0) {
*isneg = 1;
d = -d;
} else
*isneg = 0;
*hd = (int)d;
min = (d - *hd)*60.;
*m = (int)min;
*s = (int)((min - *m)*60. + 0.5);
if (*s == 60) {
if ((*m += 1) == 60) {
*hd += 1;
*m = 0;
}
*s = 0;
}
/* no negative 0's */
if (*hd == 0 && *m == 0 && *s == 0)
*isneg = 0;
}
/* insure 0 <= *v < r.
*/
range (v, r)
float *v, r;
{
while (*v < 0) *v += r;
while (*v >= r) *v -= r;
}
xXx
echo extracting sun.c
cat > sun.c << 'xXx'
#include <stdio.h>
#include <math.h>
#include "astro.h"
/* given the modified JD, mjd, return the true geocentric ecliptic longitude
* of the sun for the mean equinox of the date, *lsn, in radians, and the
* sun-earth distance, *rsn, in AU. (the true ecliptic latitude is never more
* than 1.2 arc seconds and so may be taken to be a constant 0.)
* if the APPARENT ecliptic longitude is required, correct the longitude for
* nutation to the true equinox of date and for aberration (light travel time,
* approximately -9.27e7/186000/(3600*24*365)*2*pi = -9.93e-5 radians).
*/
sun (mjd, lsn, rsn)
float mjd;
float *lsn, *rsn;
{
float t, t2;
float ls, ms; /* mean longitude and mean anomoay */
float s, nu, ea; /* eccentricity, true anomaly, eccentric anomaly */
double a, b, a1, b1, c1, d1, e1, h1, dl, dr;
t = mjd/36525.;
t2 = t*t;
a = 100.0021359*t;
b = 360.*(a-(int)a);
ls = 279.69668+.0003025*t2+b;
a = 99.99736042000039*t;
b = 360*(a-(int)a);
ms = 358.47583-(.00015+.0000033*t)*t2+b;
s = .016751-.0000418*t-1.26e-07*t2;
anomaly (degrad(ms), s, &nu, &ea);
a = 62.55209472000015*t;
b = 360*(a-(int)a);
a1 = degrad(153.23+b);
a = 125.1041894*t;
b = 360*(a-(int)a);
b1 = degrad(216.57+b);
a = 91.56766028*t;
b = 360*(a = (int)a);
c1 = degrad(312.69+b);
a = 1236.853095*t;
b = 360*(a-(int)a);
d1 = degrad(350.74-.00144*t2+b);
e1 = degrad(231.19+20.2*t);
a = 183.1353208*t;
b = 360*(a-(int)a);
h1 = degrad(353.4+b);
dl = .00134*cos(a1)+.00154*cos(b1)+.002*cos(c1)+.00179*sin(d1)+
.00178*sin(e1);
dr = 5.43e-06*sin(a1)+1.575e-05*sin(b1)+1.627e-05*sin(c1)+
3.076e-05*cos(d1)+9.27e-06*sin(h1);
*lsn = nu+degrad(ls-ms+dl);
*rsn = 1.0000002*(1-s*cos(ea))+dr;
range (lsn, 2*PI);
}
xXx
echo extracting sunrs.c
cat > sunrs.c << 'xXx'
#include <stdio.h>
#include <math.h>
#include "astro.h"
/* given the modified julian date, mjd, geographical latitude, lat, and
* longitude, lng, each in radians with west lng < 0, and a horizon
* displacement, dis, find the utc times of sunrise, utcr, and sunset, utcs,
* and the azimuths of sunrise, azr, and sunset, azs. all angles in radians;
* azimuths are E of N.
* times are for the upper limb, including effects of nutation and aberration.
* displacement is additional amount added to suns altitude compared to its
* true location; see riset.c for more discussion. it can be used to
* account for irregular horizons, various refraction models, even times
* of twilight. eg, refraction to a level horizon is often taken to be about
* 32/2(sun's semi-diameter) + 34(nominal refraction) = 50 arc minutes.
* better is:
* unrefract (pre, temp, 0.0, &a); /* true alt of upper limb
* a -= SUN_DIAM; /* true alt of sun center
* also used to find astro twilight by calling with dis of 18 degrees.
* status:
* 2: error
* 1: never rises
* 0: normal
* -1: circumpolar
*/
sunrs (mjd, lat, lng, dis, utcr, utcs, azr, azs, status)
float mjd, lat, lng, dis;
float *utcr, *utcs;
float *azr, *azs;
int *status;
{
float lstr, lsts; /* local sidereal times of rising/setting */
float al = radhr(lng); /* time correction for longitude */
float tmp1, tmp2; /* tmp */
float mjd0, t;
mjd0 = floor(mjd+0.5)-0.5; /* insure mjd0 is greenwich start of day */
/* first approximation is to find rise/set times of a fixed object
* in the position of the sun at local midday. the sun is not
* fixed but moves about a degree per day so these are then refined.
*/
fixedsunriset(mjd0+(12.-al)/24.,lat,dis,&lstr,&lsts,&tmp1,&tmp2,status);
if (*status != 0) return;
/* find a better approximation to the rising circumstances based on a
* fixed object at the suns location at the first approximation of the
* rise time.
* N.B. more iterations are less than sp float precision.
*/
gst_utc (mjd0, lstr, &t);
fixedsunriset(mjd0+(t-al)/24.,lat,dis,&lstr,&tmp1,azr,&tmp2,status);
if (*status != 0) return;
gst_utc (mjd0, lstr, utcr);
*utcr -= al*.9972696; /* allow for sideral shift */
/* find a better approximation to the setting circumstances based on a
* fixed object at the suns location at the first approximation of the
* setting time.
*/
gst_utc (mjd0, lsts, &t);
fixedsunriset (mjd0+(t-al)/24.,lat,dis,&tmp1,&lsts,&tmp2,azs,status);
if (*status != 0) return;
gst_utc (mjd0, lsts, utcs);
*utcs -= al*.9972696; /* allow for sideral shift */
}
/* 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.
*/
static
fixedsunriset (mjd, lat, dis, lstr, lsts, azr, azs, status)
float mjd, lat, dis;
float *lstr, *lsts;
float *azr, *azs;
int *status;
{
float lsn, rsn;
float deps, dpsi;
float ra, dec;
/* find ecliptic position of sun at mjd */
sun (mjd, &lsn, &rsn);
/* allow for 20.4" aberation and nutation */
nutation (mjd, &deps, &dpsi);
lsn += dpsi - degrad(20.4/3600.0);
/* convert ecliptic to equatorial coords */
ecl_eq (mjd, 0.0, lsn, &ra, &dec);
/* find circumstances for given horizon displacement */
riset (ra, dec, lat, dis, lstr, lsts, azr, azs, status);
}
xXx
echo extracting utc_gst.c
cat > utc_gst.c << 'xXx'
/* given a modified julian date, mjd, and a universally coordinated time, utc,
* return greenwich mean siderial time, *gst.
*/
utc_gst (mjd, utc, gst)
float mjd;
float utc;
float *gst;
{
float tnaught();
static float lastmjd;
static float t0;
if (mjd != lastmjd) {
t0 = tnaught (mjd);
lastmjd = mjd;
}
*gst = 1.002737908*utc + t0;
range (gst, 24.0);
}
/* given a modified julian date, mjd, and a greenwich mean siderial time, gst,
* return universally coordinated time, *utc.
*/
gst_utc (mjd, gst, utc)
float mjd;
float gst;
float *utc;
{
float tnaught();
static float lastmjd;
static float t0;
if (mjd != lastmjd) {
t0 = tnaught (mjd);
range (&t0, 24.0);
lastmjd = mjd;
}
*utc = gst - t0;
range (utc, 24.0);
*utc *= .9972695677;
}
static float
tnaught (mjd)
float mjd; /* julian days since 1900 jan 0.5 */
{
float dmjd;
int m, y;
float d;
float t, t0;
mjd_cal (mjd, &m, &d, &y);
cal_mjd (1, 0., y, &dmjd);
t = dmjd/36525;
t0 = 6.57098e-2 * (mjd - dmjd) -
(24 - (6.6460656 + (5.1262e-2 + (t * 2.581e-5))*t) -
(2400 * (t - (((float)y - 1900)/100))));
return (t0);
}
xXx