downey@cs.umn.edu@dimed1.UUCP (08/31/90)
Posting-number: Volume 14, Issue 76
Submitted-by: downey@cs.umn.edu@dimed1.UUCP
Archive-name: ephem-4.21/part01
Please accept this latest posting of my interactive astronomical ephemeris
program, ephem, v4.21 August 21, 1990. This distribution consists of
six shar files and a manual for a total of seven installments.
These are the changes since my last posting (V4.13):
add g/k and H/G magnitude model options for elliptical objects.
put moon's geocentric long/lat under Hlong/Hlat columns.
allow entering negative decimal years.
init all static mjd storage to unlike times.
add USE_TERMIO option to io.c.
add listing feature, with 'L' hot-key.
add title for plot file (as well as listing file).
add some (void) casts for lint sake.
fix parabolic comet bug in objx.c (bad lam computation).
add 'c' hot-key to Menu field.
display full Dec precision for fixed objx setup.
increase pluto auscale in watch.c, and guard screen boundries.
add Pause field.
further improve rise/set and dawn/dusk times.
add MENU={DATA,RISET,SEP} config/arg option.
watch popup now allows changing formats without returning.
add 'w' hot-key to watch field.
improve labeling a bit in Dome display.
move setjmp() in main so it catches error in config file too.
maintain name of objx/y and generally clean up objx.c.
fix bug in circum.c related to phase of fixed objects.
add "Sky dome" watch format and improve labels.
remember last watch, plot and search popup selections.
fix bug watching object y by adding y to body_tags[] in watch.c.
add ! support (#ifdef BANG in sel_fld()).
add support for VMS via #ifdef VMS.
switch to EPHEMCFG environ variable (HOME no longer used).
fix phase so it works for objects out of the ecliptic.
Please feel free to contact me if you have any problems or suggestions.
Elwood Downey
downey@dimed.com
-----CUT--------------------------------------------------------------------
#! /bin/sh
# This is a shell archive. Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file". To overwrite existing
# files, type "sh file -c". You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g.. If this archive is complete, you
# will see the following message at the end:
# "End of archive 1 (of 6)."
# Contents: Makefile Manifest aa_hadec.c anomaly.c astro.h cal_mjd.c
# circum.h comet.c ephem.cfg eq_ecl.c flog.c moonnf.c nutation.c
# obliq.c parallax.c pelement.c popup.c precess.c reduce.c refract.c
# riset.c sex_dec.c sun.c time.c utc_gst.c
# Wrapped by allbery@uunet on Thu Aug 30 20:46:27 1990
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'Makefile' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'Makefile'\"
else
echo shar: Extracting \"'Makefile'\" \(1278 characters\)
sed "s/^X//" >'Makefile' <<'END_OF_FILE'
X# Makefile for ephem
X
XCFLAGS = -O
X
XEPHEM= aa_hadec.o altmenus.o anomaly.o cal_mjd.o circum.o comet.o compiler.o \
X eq_ecl.o flog.o formats.o io.o listing.o main.o mainmenu.o moon.o \
X moonnf.o nutation.o objx.o obliq.o parallax.o pelement.o plans.o \
X plot.o popup.o precess.o reduce.o refract.o riset.o riset_c.o \
X sel_fld.o sex_dec.o srch.o sun.o time.o utc_gst.o version.o watch.o
X
Xephem: $(EPHEM)
X cc -o $@ $(EPHEM) -ltermcap -lm
X
Xaa_hadec.o: astro.h
X
Xaltmenus.o: astro.h circum.h screen.h
X
Xanomaly.o: astro.h
X
Xcal_mjd.o: astro.h
X
Xcircum.o: astro.h circum.h screen.h
X
Xcomet.o: astro.h
X
Xcompiler.o: screen.h
X
Xeq_ecl.o: astro.h
X
Xflog.o: screen.h
X
Xformats.o: astro.h screen.h
X
Xio.o: screen.h
X
Xlisting.o: screen.h
X
Xmain.o: astro.h circum.h screen.h
X
Xmainmenu.o: astro.h circum.h screen.h
X
Xmoon.o: astro.h
X
Xmoonnf.o: astro.h
X
Xnutation.o: astro.h
X
Xobjx.o: astro.h circum.h screen.h
X
Xobliq.o: astro.h
X
Xparallax.o: astro.h
X
Xpelement.o: astro.h
X
Xplans.o: astro.h
X
Xplot.o: screen.h
X
Xpopup.o: screen.h
X
Xprecess.o: astro.h
X
Xreduce.o: astro.h
X
Xrefract.o: astro.h
X
Xriset.o: astro.h
X
Xriset_c.o: astro.h circum.h screen.h
X
Xsel_fld.o: screen.h
X
Xsrch.o: screen.h
X
Xsun.o: astro.h
X
Xtime.o: astro.h circum.h
X
Xutc_gst.o: astro.h
X
Xversion.o: screen.h
X
Xwatch.o: astro.h circum.h screen.h
END_OF_FILE
if test 1278 -ne `wc -c <'Makefile'`; then
echo shar: \"'Makefile'\" unpacked with wrong size!
fi
# end of 'Makefile'
fi
if test -f 'Manifest' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'Manifest'\"
else
echo shar: Extracting \"'Manifest'\" \(2465 characters\)
sed "s/^X//" >'Manifest' <<'END_OF_FILE'
XFiles shipped as the ephem v4.20 package:
X
XMakefile just say "make"
XMan.txt manual, ready for a very generic printer.
XManifest this file.
XReadme check here for hints before building.
Xephem.cfg sample configuration file. you'll want to edit this for your
X particular location and interests after reading the manual.
Xephem.db sample database file. ditto.
X
XThe source files:
Xaa_hadec.c convert between alt/az and hour angle/dec.
Xaltmenus.c draws the three alternate lower screens.
Xanomaly.c compute anomaly.
Xastro.h unit conversion macros and planet defines.
Xcal_mjd.c converters to and from modified julian date.
Xcircum.c main "astronomy" entry point that finds where anything is.
Xcircum.h defines Now and Sky structures.
Xcomet.c compute comet position from elements.
Xcompiler.c compile and execute general expressions with screen fields.
Xeq_ecl.c convert between equitorial and eclipitic coords.
Xflog.c logs and retrieves screen locations for logging purposes.
Xformats.c basic date, time, prompts, etc formats.
Xio.c basic cursor and screen io control.
Xlisting.c set fields for and manage listings.
Xmain.c main loop.
Xmainmenu.c draws the top screen.
Xmoon.c compute moon position.
Xmoonnf.c compute new and full moon dates.
Xnutation.c compute nutation correction.
Xobjx.c set and make use of all the object-x info.
Xobliq.c compute obliquity.
Xparallax.c functions to compute earth rim parallax correction.
Xpelement.c basic planet position polynomial coefficients.
Xplans.c use polynomials to find planet location at any certain time.
Xplot.c set fields for and manage plots.
Xpopup.c handle the one-liner "popup" menus at the top of the screen.
Xprecess.c compute precession correction.
Xreduce.c convert elliptical elements from one epoch to another.
Xrefract.c atmospheric refraction model.
Xriset.c find basic rise/set sideral times of a fixed object.
Xriset_c.c iteratively solve for local rise/set times of moving objects.
Xscreen.h define all field locations and codes, extra planet codes.
Xsel_fld.c handle cursor movement commands in the various screens.
Xsex_dec.c convert between sexagesimal and decimal notation.
Xsrch.c set and manage the various search functions.
Xsun.c compute location of sun at any time.
Xtime.c manage setting and getting the time from the os.
Xutc_gst.c convert between UT1 and Greenwich sidereal time.
Xversion.c current version notice, and revision history comments.
Xwatch.c manage the screen during sky and solar system displays.
END_OF_FILE
if test 2465 -ne `wc -c <'Manifest'`; then
echo shar: \"'Manifest'\" unpacked with wrong size!
fi
# end of 'Manifest'
fi
if test -f 'aa_hadec.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'aa_hadec.c'\"
else
echo shar: Extracting \"'aa_hadec.c'\" \(1922 characters\)
sed "s/^X//" >'aa_hadec.c' <<'END_OF_FILE'
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X
X/* given latitude (n+, radians), lat, altitude (up+, radians), alt, and
X * azimuth (angle round to the east from north+, radians),
X * return hour angle (radians), ha, and declination (radians), dec.
X */
Xaa_hadec (lat, alt, az, ha, dec)
Xdouble lat;
Xdouble alt, az;
Xdouble *ha, *dec;
X{
X aaha_aux (lat, az, alt, ha, dec);
X}
X
X/* given latitude (n+, radians), lat, hour angle (radians), ha, and declination
X * (radians), dec,
X * return altitude (up+, radians), alt, and
X * azimuth (angle round to the east from north+, radians),
X */
Xhadec_aa (lat, ha, dec, alt, az)
Xdouble lat;
Xdouble ha, dec;
Xdouble *alt, *az;
X{
X aaha_aux (lat, ha, dec, az, alt);
X}
X
X/* the actual formula is the same for both transformation directions so
X * do it here once for each way.
X * N.B. all arguments are in radians.
X */
Xstatic
Xaaha_aux (lat, x, y, p, q)
Xdouble lat;
Xdouble x, y;
Xdouble *p, *q;
X{
X static double lastlat = -1000.;
X static double sinlastlat, coslastlat;
X double sy, cy;
X double sx, cx;
X double sq, cq;
X double a;
X double cp;
X
X /* latitude doesn't change much, so try to reuse the sin and cos evals.
X */
X if (lat != lastlat) {
X sinlastlat = sin (lat);
X coslastlat = cos (lat);
X lastlat = lat;
X }
X
X sy = sin (y);
X cy = cos (y);
X sx = sin (x);
X cx = cos (x);
X
X/* define GOODATAN2 if atan2 returns full range -PI through +PI.
X */
X#ifdef GOODATAN2
X *q = asin ((sy*sinlastlat) + (cy*coslastlat*cx));
X *p = atan2 (-cy*sx, -cy*cx*sinlastlat + sy*coslastlat);
X#else
X#define EPS (1e-20)
X sq = (sy*sinlastlat) + (cy*coslastlat*cx);
X *q = asin (sq);
X cq = cos (*q);
X a = coslastlat*cq;
X if (a > -EPS && a < EPS)
X a = a < 0 ? -EPS : EPS; /* avoid / 0 */
X cp = (sy - (sinlastlat*sq))/a;
X if (cp >= 1.0) /* the /a can be slightly > 1 */
X *p = 0.0;
X else if (cp <= -1.0)
X *p = PI;
X else
X *p = acos ((sy - (sinlastlat*sq))/a);
X if (sx>0) *p = 2.0*PI - *p;
X#endif
X}
END_OF_FILE
if test 1922 -ne `wc -c <'aa_hadec.c'`; then
echo shar: \"'aa_hadec.c'\" unpacked with wrong size!
fi
# end of 'aa_hadec.c'
fi
if test -f 'anomaly.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'anomaly.c'\"
else
echo shar: Extracting \"'anomaly.c'\" \(557 characters\)
sed "s/^X//" >'anomaly.c' <<'END_OF_FILE'
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X
X#define TWOPI (2*PI)
X
X/* given the mean anomaly, ma, and the eccentricity, s, of elliptical motion,
X * find the true anomaly, *nu, and the eccentric anomaly, *ea.
X * all angles in radians.
X */
Xanomaly (ma, s, nu, ea)
Xdouble ma, s;
Xdouble *nu, *ea;
X{
X double m, dla, fea;
X
X m = ma-TWOPI*(long)(ma/TWOPI);
X fea = m;
X while (1) {
X dla = fea-(s*sin(fea))-m;
X if (fabs(dla)<1e-6)
X break;
X dla /= 1-(s*cos(fea));
X fea -= dla;
X }
X
X *nu = 2*atan(sqrt((1+s)/(1-s))*tan(fea/2));
X *ea = fea;
X}
END_OF_FILE
if test 557 -ne `wc -c <'anomaly.c'`; then
echo shar: \"'anomaly.c'\" unpacked with wrong size!
fi
# end of 'anomaly.c'
fi
if test -f 'astro.h' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'astro.h'\"
else
echo shar: Extracting \"'astro.h'\" \(620 characters\)
sed "s/^X//" >'astro.h' <<'END_OF_FILE'
X#ifndef PI
X#define PI 3.141592653589793
X#endif
X
X/* conversions among hours (of ra), degrees and radians. */
X#define degrad(x) ((x)*PI/180.)
X#define raddeg(x) ((x)*180./PI)
X#define hrdeg(x) ((x)*15.)
X#define deghr(x) ((x)/15.)
X#define hrrad(x) degrad(hrdeg(x))
X#define radhr(x) deghr(raddeg(x))
X
X/* ratio of from synodic (solar) to sidereal (stellar) rate */
X#define SIDRATE .9972695677
X
X/* manifest names for planets.
X * N.B. must cooincide with usage in pelement.c and plans.c.
X */
X#define MERCURY 0
X#define VENUS 1
X#define MARS 2
X#define JUPITER 3
X#define SATURN 4
X#define URANUS 5
X#define NEPTUNE 6
X#define PLUTO 7
END_OF_FILE
if test 620 -ne `wc -c <'astro.h'`; then
echo shar: \"'astro.h'\" unpacked with wrong size!
fi
# end of 'astro.h'
fi
if test -f 'cal_mjd.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'cal_mjd.c'\"
else
echo shar: Extracting \"'cal_mjd.c'\" \(3067 characters\)
sed "s/^X//" >'cal_mjd.c' <<'END_OF_FILE'
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X
X/* given a date in months, mn, days, dy, years, yr,
X * return the modified Julian date (number of days elapsed since 1900 jan 0.5),
X * *mjd.
X */
Xcal_mjd (mn, dy, yr, mjd)
Xint mn, yr;
Xdouble dy;
Xdouble *mjd;
X{
X int b, d, m, y;
X long c;
X
X m = mn;
X y = (yr < 0) ? yr + 1 : yr;
X if (mn < 3) {
X m += 12;
X y -= 1;
X }
X
X if (yr < 1582 || yr == 1582 && (mn < 10 || mn == 10 && dy < 15))
X b = 0;
X else {
X int a;
X a = y/100;
X b = 2 - a + a/4;
X }
X
X if (y < 0)
X c = (long)((365.25*y) - 0.75) - 694025L;
X else
X c = (long)(365.25*y) - 694025L;
X
X d = 30.6001*(m+1);
X
X *mjd = b + c + d + dy - 0.5;
X}
X
X/* given the modified Julian date (number of days elapsed since 1900 jan 0.5,),
X * mjd, return the calendar date in months, *mn, days, *dy, and years, *yr.
X */
Xmjd_cal (mjd, mn, dy, yr)
Xdouble mjd;
Xint *mn, *yr;
Xdouble *dy;
X{
X double d, f;
X double i, a, b, ce, g;
X
X d = mjd + 0.5;
X i = floor(d);
X f = d-i;
X if (f == 1) {
X f = 0;
X i += 1;
X }
X
X if (i > -115860.0) {
X a = floor((i/36524.25)+.9983573)+14;
X i += 1 + a - floor(a/4.0);
X }
X
X b = floor((i/365.25)+.802601);
X ce = i - floor((365.25*b)+.750001)+416;
X g = floor(ce/30.6001);
X *mn = g - 1;
X *dy = ce - floor(30.6001*g)+f;
X *yr = b + 1899;
X
X if (g > 13.5)
X *mn = g - 13;
X if (*mn < 2.5)
X *yr = b + 1900;
X if (*yr < 1)
X *yr -= 1;
X}
X
X/* given an mjd, set *dow to 0..6 according to which dayof the week it falls
X * on (0=sunday) or set it to -1 if can't figure it out.
X */
Xmjd_dow (mjd, dow)
Xdouble mjd;
Xint *dow;
X{
X /* cal_mjd() uses Gregorian dates on or after Oct 15, 1582.
X * (Pope Gregory XIII dropped 10 days, Oct 5..14, and improved the leap-
X * year algorithm). however, Great Britian and the colonies did not
X * adopt it until Sept 14, 1752 (they dropped 11 days, Sept 3-13,
X * due to additional accumulated error). leap years before 1752 thus
X * can not easily be accounted for from the cal_mjd() number...
X */
X if (mjd < -53798.5) {
X /* pre sept 14, 1752 too hard to correct */
X *dow = -1;
X return;
X }
X *dow = ((long)floor(mjd-.5) + 1) % 7;/* 1/1/1900 (mjd 0.5) is a Monday*/
X if (*dow < 0)
X *dow += 7;
X}
X
X/* given a mjd, return the the number of days in the month. */
Xmjd_dpm (mjd, ndays)
Xdouble mjd;
Xint *ndays;
X{
X static short dpm[] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
X int m, y;
X double d;
X
X mjd_cal (mjd, &m, &d, &y);
X *ndays = (m==2 && ((y%4==0 && y%100!=0)||y%400==0)) ? 29 : dpm[m-1];
X}
X
X
X/* given a mjd, return the year as a double. */
Xmjd_year (mjd, yr)
Xdouble mjd;
Xdouble *yr;
X{
X int m, y;
X double d;
X double e0, e1; /* mjd of start of this year, start of next year */
X
X mjd_cal (mjd, &m, &d, &y);
X cal_mjd (1, 1.0, y, &e0);
X cal_mjd (1, 1.0, y+1, &e1);
X *yr = y + (mjd - e0)/(e1 - e0);
X}
X
X/* given a decimal year, return mjd */
Xyear_mjd (y, mjd)
Xdouble y;
Xdouble *mjd;
X{
X double e0, e1; /* mjd of start of this year, start of next year */
X int yf = floor (y);
X
X cal_mjd (1, 1.0, yf, &e0);
X cal_mjd (1, 1.0, yf+1, &e1);
X *mjd = e0 + (y - yf)*(e1-e0);
X}
END_OF_FILE
if test 3067 -ne `wc -c <'cal_mjd.c'`; then
echo shar: \"'cal_mjd.c'\" unpacked with wrong size!
fi
# end of 'cal_mjd.c'
fi
if test -f 'circum.h' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'circum.h'\"
else
echo shar: Extracting \"'circum.h'\" \(2856 characters\)
sed "s/^X//" >'circum.h' <<'END_OF_FILE'
X#define SPD (24.0*3600.0) /* seconds per day */
X
X#define EOD (-9786) /* special epoch flag: use epoch of date */
X#define RTC (-1324) /* special tminc flag: use rt clock */
X#define NOMJD (-58631.) /* an unlikely mjd for initing static mjd's */
X#define NOHELIO (-2314) /* special s_hlong flag: means it and s_hlat are
X * undefined
X */
X
X#define STDHZN 0 /* rise/set times based on nominal conditions */
X#define ADPHZN 1 /* rise/set times based on exact current " */
X#define TWILIGHT 2 /* rise/set times for sun 18 degs below hor */
X
X/* info about our local observing circumstances */
Xtypedef struct {
X double n_mjd; /* modified Julian date, ie, days since
X * Jan 0.5 1900 (== 12 noon, Dec 30, 1899), utc.
X * enough precision to get well better than 1 second.
X * N.B. if not first member, must move NOMJD inits.
X */
X double n_lat; /* latitude, >0 north, rads */
X double n_lng; /* longitude, >0 east, rads */
X double n_tz; /* time zone, hrs behind UTC */
X double n_temp; /* atmospheric temp, degrees C */
X double n_pressure; /* atmospheric pressure, mBar */
X double n_height; /* height above sea level, earth radii */
X double n_epoch; /* desired precession display epoch as an mjd, or EOD */
X char n_tznm[4]; /* time zone name; 3 chars or less, always 0 at end */
X} Now;
Xextern double mjd_day(), mjd_hr();
X
X/* info about where and how we see something in the sky */
Xtypedef struct {
X double s_ra; /* ra, rads (precessed to n_epoch) */
X double s_dec; /* dec, rads (precessed to n_epoch) */
X double s_az; /* azimuth, >0 e of n, rads */
X double s_alt; /* altitude above topocentric horizon, rads */
X double s_sdist; /* dist from object to sun, au */
X double s_edist; /* dist from object to earth, au */
X double s_elong; /* angular sep between object and sun, >0 if east */
X double s_hlong; /* heliocentric longitude, rads */
X double s_hlat; /* heliocentric latitude, rads */
X double s_size; /* angular size, arc secs */
X double s_phase; /* phase, % */
X double s_mag; /* visual magnitude */
X} Sky;
X
X/* flags for riset_cir() status */
X#define RS_NORISE 0x001 /* object does not rise as such today */
X#define RS_2RISES 0x002 /* object rises more than once today */
X#define RS_NOSET 0x004 /* object does not set as such today */
X#define RS_2SETS 0x008 /* object sets more than once today */
X#define RS_CIRCUMPOLAR 0x010 /* object stays up all day today */
X#define RS_2TRANS 0x020 /* transits twice in one day */
X#define RS_NEVERUP 0x040 /* object never rises today */
X#define RS_NOTRANS 0x080 /* doesn't transit today */
X#define RS_ERROR 0x100 /* can't figure out times... */
X
X/* shorthands for fields a Now pointer, np */
X#define mjd np->n_mjd
X#define lat np->n_lat
X#define lng np->n_lng
X#define tz np->n_tz
X#define temp np->n_temp
X#define pressure np->n_pressure
X#define height np->n_height
X#define epoch np->n_epoch
X#define tznm np->n_tznm
END_OF_FILE
if test 2856 -ne `wc -c <'circum.h'`; then
echo shar: \"'circum.h'\" unpacked with wrong size!
fi
# end of 'circum.h'
fi
if test -f 'comet.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'comet.c'\"
else
echo shar: Extracting \"'comet.c'\" \(2390 characters\)
sed "s/^X//" >'comet.c' <<'END_OF_FILE'
X#include <math.h>
X#include "astro.h"
X
X/* given a modified Julian date, mjd, and a set of heliocentric parabolic
X * orbital elements referred to the epoch of date (mjd):
X * ep: epoch of perihelion,
X * inc: inclination,
X * ap: argument of perihelion (equals the longitude of perihelion minus the
X * longitude of ascending node)
X * qp: perihelion distance,
X * om: longitude of ascending node;
X * find:
X * lpd: heliocentric longitude,
X * psi: heliocentric latitude,
X * rp: distance from the sun to the planet,
X * rho: distance from the Earth to the planet,
X * lam: geocentric ecliptic longitude,
X * bet: geocentric ecliptic latitude,
X * none are corrected for light time, ie, they are the true values for
X * the given instant.
X *
X * all angles are in radians, all distances in AU.
X * mutual perturbation corrections with other solar system objects are not
X * applied. corrections for nutation and abberation must be made by the caller.
X * The RA and DEC calculated from the fully-corrected ecliptic coordinates are
X * then the apparent geocentric coordinates. Further corrections can be made,
X * if required, for atmospheric refraction and geocentric parallax.
X */
Xcomet (mjd, ep, inc, ap, qp, om, lpd, psi, rp, rho, lam, bet)
Xdouble mjd;
Xdouble ep, inc, ap, qp, om;
Xdouble *lpd, *psi, *rp, *rho, *lam, *bet;
X{
X double w, s, s2;
X double l, sl, cl, y;
X double spsi, cpsi;
X double rd, lsn, rsn;
X double lg, re, ll;
X double cll, sll;
X double nu;
X
X#define ERRLMT 0.0001
X w = ((mjd-ep)*3.649116e-02)/(qp*sqrt(qp));
X s = w/3;
X while (1) {
X double d;
X s2 = s*s;
X d = (s2+3)*s-w;
X if (fabs(d) <= ERRLMT)
X break;
X s = ((2*s*s2)+w)/(3*(s2+1));
X }
X
X nu = 2*atan(s);
X *rp = qp*(1+s2);
X l = nu+ap;
X sl = sin(l);
X cl = cos(l);
X spsi = sl*sin(inc);
X *psi = asin(spsi);
X y = sl*cos(inc);
X *lpd = atan(y/cl)+om;
X cpsi = cos(*psi);
X if (cl<0) *lpd += PI;
X range (lpd, 2*PI);
X rd = *rp * cpsi;
X sunpos (mjd, &lsn, &rsn);
X lg = lsn+PI;
X re = rsn;
X ll = *lpd - lg;
X cll = cos(ll);
X sll = sin(ll);
X *rho = sqrt((re * re)+(*rp * *rp)-(2*re*rd*cll));
X if (rd<re)
X *lam = atan((-1*rd*sll)/(re-(rd*cll)))+lg+PI;
X else
X *lam = atan((re*sll)/(rd-(re*cll)))+*lpd;
X range (lam, 2*PI);
X *bet = atan((rd*spsi*sin(*lam-*lpd))/(cpsi*re*sll));
X}
END_OF_FILE
if test 2390 -ne `wc -c <'comet.c'`; then
echo shar: \"'comet.c'\" unpacked with wrong size!
fi
# end of 'comet.c'
fi
if test -f 'ephem.cfg' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'ephem.cfg'\"
else
echo shar: Extracting \"'ephem.cfg'\" \(140 characters\)
sed "s/^X//" >'ephem.cfg' <<'END_OF_FILE'
XUT=Now
XLONG=93:42:8
XLAT=44:50:37
XHEIGHT=800
XTEMP=40
XPRES=29.5
XSTPSZ=RTC
XPROPTS=TSMevmjsunpxy
XEPOCH=2000
XNSTEP=1
XPAUSE=0
X
XOBJX=Levy
XOBJY=HMP
END_OF_FILE
if test 140 -ne `wc -c <'ephem.cfg'`; then
echo shar: \"'ephem.cfg'\" unpacked with wrong size!
fi
# end of 'ephem.cfg'
fi
if test -f 'eq_ecl.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'eq_ecl.c'\"
else
echo shar: Extracting \"'eq_ecl.c'\" \(1899 characters\)
sed "s/^X//" >'eq_ecl.c' <<'END_OF_FILE'
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X
X#define EQtoECL 1
X#define ECLtoEQ (-1)
X
X/* given the modified Julian date, mjd, and an equitorial ra and dec, each in
X * radians, find the corresponding geocentric ecliptic latitude, *lat, and
X * longititude, *lng, also each in radians.
X * correction for the effect on the angle of the obliquity due to nutation is
X * included.
X */
Xeq_ecl (mjd, ra, dec, lat, lng)
Xdouble mjd, ra, dec;
Xdouble *lat, *lng;
X{
X ecleq_aux (EQtoECL, mjd, ra, dec, lng, lat);
X}
X
X/* given the modified Julian date, mjd, and a geocentric ecliptic latitude,
X * *lat, and longititude, *lng, each in radians, find the corresponding
X * equitorial ra and dec, also each in radians.
X * correction for the effect on the angle of the obliquity due to nutation is
X * included.
X */
Xecl_eq (mjd, lat, lng, ra, dec)
Xdouble mjd, lat, lng;
Xdouble *ra, *dec;
X{
X ecleq_aux (ECLtoEQ, mjd, lng, lat, ra, dec);
X}
X
Xstatic
Xecleq_aux (sw, mjd, x, y, p, q)
Xint sw; /* +1 for eq to ecliptic, -1 for vv. */
Xdouble mjd, x, y; /* sw==1: x==ra, y==dec. sw==-1: x==lng, y==lat. */
Xdouble *p, *q; /* sw==1: p==lng, q==lat. sw==-1: p==ra, q==dec. */
X{
X static double lastmjd = -10000; /* last mjd calculated */
X static double seps, ceps; /* sin and cos of mean obliquity */
X double sx, cx, sy, cy, ty;
X
X if (mjd != lastmjd) {
X double eps;
X double deps, dpsi;
X obliquity (mjd, &eps); /* mean obliquity for date */
X nutation (mjd, &deps, &dpsi);
X eps += deps;
X seps = sin(eps);
X ceps = cos(eps);
X lastmjd = mjd;
X }
X
X sy = sin(y);
X cy = cos(y); /* always non-negative */
X if (fabs(cy)<1e-20) cy = 1e-20; /* insure > 0 */
X ty = sy/cy;
X cx = cos(x);
X sx = sin(x);
X *q = asin((sy*ceps)-(cy*seps*sx*sw));
X *p = atan(((sx*ceps)+(ty*seps*sw))/cx);
X if (cx<0) *p += PI; /* account for atan quad ambiguity */
X range (p, 2*PI);
X}
END_OF_FILE
if test 1899 -ne `wc -c <'eq_ecl.c'`; then
echo shar: \"'eq_ecl.c'\" unpacked with wrong size!
fi
# end of 'eq_ecl.c'
fi
if test -f 'flog.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'flog.c'\"
else
echo shar: Extracting \"'flog.c'\" \(3241 characters\)
sed "s/^X//" >'flog.c' <<'END_OF_FILE'
X/* this is a simple little package to manage the saving and retrieving of
X * field values, which we call field logging or "flogs". a flog consists of a
X * field location, ala rcfpack(), its value as a double and its value as
X * a string (ie, however it was printed). you can reset the list of flogs, add
X * to and remove from the list of registered fields and log a field if it has
X * been registered.
X *
X * this is used by the plotting and searching facilities of ephem to maintain
X * the values of the fields that are being plotted or used in search
X * expressions. it is used by the listing facility to generate listing files.
X *
X * a field can be in use for more than one
X * thing at a time (eg, all the X plot values may the same time field, or
X * searching and plotting might be on at one time using the same field) so
X * we consider the field to be in use as long a usage count is > 0.
X */
X
X#include "screen.h"
X
Xextern char *strcpy(), *strncpy();
X
X#define NFLOGS 32 /* max number of distinct simultaneous logged
X * fields
X */
X
Xtypedef struct {
X int fl_usagecnt; /* number of "users" logging to this field */
X int fl_fld; /* an rcfpack(r,c,0) */
X double fl_val; /* stored value as a double */
X char fl_str[16]; /* stored value as a formatted string.
X * N.B.: never overwrite last char: keep as \0
X */
X} FLog;
X
Xstatic FLog flog[NFLOGS];
X
X/* add fld to the list. if already there, just increment usage count.
X * return 0 if ok, else -1 if no more room.
X */
Xflog_add (fld)
Xint fld;
X{
X FLog *flp, *unusedflp = 0;
X
X /* scan for fld already in list, or find an unused one along the way */
X for (flp = &flog[NFLOGS]; --flp >= flog; ) {
X if (flp->fl_usagecnt > 0) {
X if (flp->fl_fld == fld) {
X flp->fl_usagecnt++;
X return (0);
X }
X } else
X unusedflp = flp;
X }
X if (unusedflp) {
X unusedflp->fl_fld = fld;
X unusedflp->fl_usagecnt = 1;
X return (0);
X }
X return (-1);
X}
X
X/* decrement usage count for flog for fld. if goes to 0 take it out of list.
X * ok if not in list i guess...
X */
Xflog_delete (fld)
Xint fld;
X{
X FLog *flp;
X
X for (flp = &flog[NFLOGS]; --flp >= flog; )
X if (flp->fl_fld == fld && flp->fl_usagecnt > 0) {
X if (--flp->fl_usagecnt <= 0) {
X flp->fl_usagecnt = 0;
X }
X break;
X }
X}
X
X/* if plotting, listing or searching is active then
X * if rcfpack(r,c,0) is in the fld list, set its value to val.
X * return 0 if ok, else -1 if not in list.
X */
Xflog_log (r, c, val, str)
Xint r, c;
Xdouble val;
Xchar *str;
X{
X if (plot_ison() || listing_ison() || srch_ison()) {
X FLog *flp;
X int fld = rcfpack (r, c, 0);
X for (flp = &flog[NFLOGS]; --flp >= flog; )
X if (flp->fl_fld == fld && flp->fl_usagecnt > 0) {
X flp->fl_val = val;
X (void) strncpy (flp->fl_str, str, sizeof(flp->fl_str)-1);
X return(0);
X }
X return (-1);
X } else
X return (0);
X}
X
X/* search for fld in list. if find it, return its value and str, if str.
X * return 0 if found it, else -1 if not in list.
X */
Xflog_get (fld, vp, str)
Xint fld;
Xdouble *vp;
Xchar *str;
X{
X FLog *flp;
X
X for (flp = &flog[NFLOGS]; --flp >= flog; )
X if (flp->fl_fld == fld && flp->fl_usagecnt > 0) {
X *vp = flp->fl_val;
X if (str)
X (void) strcpy (str, flp->fl_str);
X return (0);
X }
X return (-1);
X}
END_OF_FILE
if test 3241 -ne `wc -c <'flog.c'`; then
echo shar: \"'flog.c'\" unpacked with wrong size!
fi
# end of 'flog.c'
fi
if test -f 'moonnf.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'moonnf.c'\"
else
echo shar: Extracting \"'moonnf.c'\" \(1557 characters\)
sed "s/^X//" >'moonnf.c' <<'END_OF_FILE'
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X
X#define unw(w,z) ((w)-floor((w)/(z))*(z))
X
X/* given a modified Julian date, mjd, return the mjd of the new
X * and full moons about then, mjdn and mjdf.
X * TODO: exactly which ones does it find? eg:
X * 5/28/1988 yields 5/15 and 5/31
X * 5/29 6/14 6/29
X */
Xmoonnf (mjd, mjdn, mjdf)
Xdouble mjd;
Xdouble *mjdn, *mjdf;
X{
X int mo, yr;
X double dy;
X double mjd0;
X double k, tn, tf, t;
X
X mjd_cal (mjd, &mo, &dy, &yr);
X cal_mjd (1, 0., yr, &mjd0);
X k = (yr-1900+((mjd-mjd0)/365))*12.3685;
X k = floor(k+0.5);
X tn = k/1236.85;
X tf = (k+0.5)/1236.85;
X t = tn;
X m (t, k, mjdn);
X t = tf;
X k += 0.5;
X m (t, k, mjdf);
X}
X
Xstatic
Xm (t, k, mjd)
Xdouble t, k;
Xdouble *mjd;
X{
X double t2, a, a1, b, b1, c, ms, mm, f, ddjd;
X
X t2 = t*t;
X a = 29.53*k;
X c = degrad(166.56+(132.87-9.173e-3*t)*t);
X b = 5.8868e-4*k+(1.178e-4-1.55e-7*t)*t2+3.3e-4*sin(c)+7.5933E-1;
X ms = 359.2242+360*unw(k/1.236886e1,1)-(3.33e-5+3.47e-6*t)*t2;
X mm = 306.0253+360*unw(k/9.330851e-1,1)+(1.07306e-2+1.236e-5*t)*t2;
X f = 21.2964+360*unw(k/9.214926e-1,1)-(1.6528e-3+2.39e-6*t)*t2;
X ms = unw(ms,360);
X mm = unw(mm,360);
X f = unw(f,360);
X ms = degrad(ms);
X mm = degrad(mm);
X f = degrad(f);
X ddjd = (1.734e-1-3.93e-4*t)*sin(ms)+2.1e-3*sin(2*ms)
X -4.068e-1*sin(mm)+1.61e-2*sin(2*mm)-4e-4*sin(3*mm)
X +1.04e-2*sin(2*f)-5.1e-3*sin(ms+mm)-7.4e-3*sin(ms-mm)
X +4e-4*sin(2*f+ms)-4e-4*sin(2*f-ms)-6e-4*sin(2*f+mm)
X +1e-3*sin(2*f-mm)+5e-4*sin(ms+2*mm);
X a1 = (long)a;
X b = b+ddjd+(a-a1);
X b1 = (long)b;
X a = a1+b1;
X b = b-b1;
X *mjd = a + b;
X}
END_OF_FILE
if test 1557 -ne `wc -c <'moonnf.c'`; then
echo shar: \"'moonnf.c'\" unpacked with wrong size!
fi
# end of 'moonnf.c'
fi
if test -f 'nutation.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'nutation.c'\"
else
echo shar: Extracting \"'nutation.c'\" \(2019 characters\)
sed "s/^X//" >'nutation.c' <<'END_OF_FILE'
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X
X/* given the modified JD, mjd, find the nutation in obliquity, *deps, and
X * the nutation in longitude, *dpsi, each in radians.
X */
Xnutation (mjd, deps, dpsi)
Xdouble mjd;
Xdouble *deps, *dpsi;
X{
X static double lastmjd = -10000, lastdeps, lastdpsi;
X double ls, ld; /* sun's mean longitude, moon's mean longitude */
X double ms, md; /* sun's mean anomaly, moon's mean anomaly */
X double nm; /* longitude of moon's ascending node */
X double t, t2; /* number of Julian centuries of 36525 days since
X * Jan 0.5 1900.
X */
X double tls, tnm, tld; /* twice above */
X double a, b; /* temps */
X
X if (mjd == lastmjd) {
X *deps = lastdeps;
X *dpsi = lastdpsi;
X return;
X }
X
X t = mjd/36525.;
X t2 = t*t;
X
X a = 100.0021358*t;
X b = 360.*(a-(long)a);
X ls = 279.697+.000303*t2+b;
X
X a = 1336.855231*t;
X b = 360.*(a-(long)a);
X ld = 270.434-.001133*t2+b;
X
X a = 99.99736056000026*t;
X b = 360.*(a-(long)a);
X ms = 358.476-.00015*t2+b;
X
X a = 13255523.59*t;
X b = 360.*(a-(long)a);
X md = 296.105+.009192*t2+b;
X
X a = 5.372616667*t;
X b = 360.*(a-(long)a);
X nm = 259.183+.002078*t2-b;
X
X /* convert to radian forms for use with trig functions.
X */
X tls = 2*degrad(ls);
X nm = degrad(nm);
X tnm = 2*degrad(nm);
X ms = degrad(ms);
X tld = 2*degrad(ld);
X md = degrad(md);
X
X /* find delta psi and eps, in arcseconds.
X */
X lastdpsi = (-17.2327-.01737*t)*sin(nm)+(-1.2729-.00013*t)*sin(tls)
X +.2088*sin(tnm)-.2037*sin(tld)+(.1261-.00031*t)*sin(ms)
X +.0675*sin(md)-(.0497-.00012*t)*sin(tls+ms)
X -.0342*sin(tld-nm)-.0261*sin(tld+md)+.0214*sin(tls-ms)
X -.0149*sin(tls-tld+md)+.0124*sin(tls-nm)+.0114*sin(tld-md);
X lastdeps = (9.21+.00091*t)*cos(nm)+(.5522-.00029*t)*cos(tls)
X -.0904*cos(tnm)+.0884*cos(tld)+.0216*cos(tls+ms)
X +.0183*cos(tld-nm)+.0113*cos(tld+md)-.0093*cos(tls-ms)
X -.0066*cos(tls-nm);
X
X /* convert to radians.
X */
X lastdpsi = degrad(lastdpsi/3600);
X lastdeps = degrad(lastdeps/3600);
X
X lastmjd = mjd;
X *deps = lastdeps;
X *dpsi = lastdpsi;
X}
END_OF_FILE
if test 2019 -ne `wc -c <'nutation.c'`; then
echo shar: \"'nutation.c'\" unpacked with wrong size!
fi
# end of 'nutation.c'
fi
if test -f 'obliq.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'obliq.c'\"
else
echo shar: Extracting \"'obliq.c'\" \(421 characters\)
sed "s/^X//" >'obliq.c' <<'END_OF_FILE'
X#include <stdio.h>
X#include "astro.h"
X
X/* given the modified Julian date, mjd, find the obliquity of the
X * ecliptic, *eps, in radians.
X */
Xobliquity (mjd, eps)
Xdouble mjd;
Xdouble *eps;
X{
X static double lastmjd = -10000, lasteps;
X
X if (mjd != lastmjd) {
X double t;
X t = mjd/36525.;
X lasteps = degrad(2.345229444E1
X - ((((-1.81E-3*t)+5.9E-3)*t+4.6845E1)*t)/3600.0);
X lastmjd = mjd;
X }
X *eps = lasteps;
X}
END_OF_FILE
if test 421 -ne `wc -c <'obliq.c'`; then
echo shar: \"'obliq.c'\" unpacked with wrong size!
fi
# end of 'obliq.c'
fi
if test -f 'parallax.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'parallax.c'\"
else
echo shar: Extracting \"'parallax.c'\" \(2280 characters\)
sed "s/^X//" >'parallax.c' <<'END_OF_FILE'
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X
X/* given true ha and dec, tha and tdec, the geographical latitude, phi, the
X * height above sea-level (as a fraction of the earths radius, 6378.16km),
X * ht, and the equatorial horizontal parallax, ehp, find the apparent
X * ha and dec, aha and adec allowing for parallax.
X * all angles in radians. ehp is the angle subtended at the body by the
X * earth's equator.
X */
Xta_par (tha, tdec, phi, ht, ehp, aha, adec)
Xdouble tha, tdec, phi, ht, ehp;
Xdouble *aha, *adec;
X{
X static double last_phi, last_ht, rsp, rcp;
X double rp; /* distance to object in Earth radii */
X double ctha;
X double stdec, ctdec;
X double tdtha, dtha;
X double caha;
X
X /* avoid calcs involving the same phi and ht */
X if (phi != last_phi || ht != last_ht) {
X double cphi, sphi, u;
X cphi = cos(phi);
X sphi = sin(phi);
X u = atan(9.96647e-1*sphi/cphi);
X rsp = (9.96647e-1*sin(u))+(ht*sphi);
X rcp = cos(u)+(ht*cphi);
X last_phi = phi;
X last_ht = ht;
X }
X
X rp = 1/sin(ehp);
X
X ctha = cos(tha);
X stdec = sin(tdec);
X ctdec = cos(tdec);
X tdtha = (rcp*sin(tha))/((rp*ctdec)-(rcp*ctha));
X dtha = atan(tdtha);
X *aha = tha+dtha;
X caha = cos(*aha);
X range (aha, 2*PI);
X *adec = atan(caha*(rp*stdec-rsp)/(rp*ctdec*ctha-rcp));
X}
X
X/* given the apparent ha and dec, aha and adec, the geographical latitude, phi,
X * the height above sea-level (as a fraction of the earths radius, 6378.16km),
X * ht, and the equatorial horizontal parallax, ehp, find the true ha and dec,
X * tha and tdec allowing for parallax.
X * all angles in radians. ehp is the angle subtended at the body by the
X * earth's equator.
X * uses ta_par() iteratively: find a set of true ha/dec that converts back
X * to the given apparent ha/dec.
X */
Xat_par (aha, adec, phi, ht, ehp, tha, tdec)
Xdouble aha, adec, phi, ht, ehp;
Xdouble *tha, *tdec;
X{
X double nha, ndec; /* ha/dec corres. to current true guesses */
X double eha, edec; /* error in ha/dec */
X
X /* first guess for true is just the apparent */
X *tha = aha;
X *tdec = adec;
X
X while (1) {
X ta_par (*tha, *tdec, phi, ht, ehp, &nha, &ndec);
X eha = aha - nha;
X edec = adec - ndec;
X if (fabs(eha)<1e-6 && fabs(edec)<1e-6)
X break;
X *tha += eha;
X *tdec += edec;
X }
X}
END_OF_FILE
if test 2280 -ne `wc -c <'parallax.c'`; then
echo shar: \"'parallax.c'\" unpacked with wrong size!
fi
# end of 'parallax.c'
fi
if test -f 'pelement.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'pelement.c'\"
else
echo shar: Extracting \"'pelement.c'\" \(4797 characters\)
sed "s/^X//" >'pelement.c' <<'END_OF_FILE'
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X
X/* this array contains polynomial coefficients to find the various orbital
X * elements for the mean orbit at any instant in time for each major planet.
X * the first five elements are in the form a0 + a1*t + a2*t**2 + a3*t**3,
X * where t is the number of Julian centuries of 36525 Julian days since 1900
X * Jan 0.5. the last three elements are constants.
X *
X * the orbital element (column) indeces are:
X * [ 0- 3]: coefficients for mean longitude, in degrees;
X * [ 4- 7]: coefficients for longitude of the perihelion, in degrees;
X * [ 8-11]: coefficients for eccentricity;
X * [12-15]: coefficients for inclination, in degrees;
X * [16-19]: coefficients for longitude of the ascending node, in degrees;
X * [20]: semi-major axis, in AU;
X * [21]: angular diameter at 1 AU, in arcsec;
X * [22]: standard visual magnitude, ie, the visual magnitude of the planet
X * when at a distance of 1 AU from both the Sun and the Earth and
X * with zero phase angle.
X *
X * the planent (row) indeces are:
X * [0]: Mercury; [1]: Venus; [2]: Mars; [3]: Jupiter; [4]: Saturn;
X * [5]: Uranus; [6]: Neptune; [7]: Pluto.
X */
X#define NPELE (5*4 + 3) /* 4 coeffs for ea of 5 elems, + 3 constants */
Xstatic double elements[8][NPELE] = {
X
X { /* mercury... */
X
X 178.179078, 415.2057519, 3.011e-4, 0.0,
X 75.899697, 1.5554889, 2.947e-4, 0.0,
X .20561421, 2.046e-5, 3e-8, 0.0,
X 7.002881, 1.8608e-3, -1.83e-5, 0.0,
X 47.145944, 1.1852083, 1.739e-4, 0.0,
X .3870986, 6.74, -0.42
X },
X
X { /* venus... */
X
X 342.767053, 162.5533664, 3.097e-4, 0.0,
X 130.163833, 1.4080361, -9.764e-4, 0.0,
X 6.82069e-3, -4.774e-5, 9.1e-8, 0.0,
X 3.393631, 1.0058e-3, -1e-6, 0.0,
X 75.779647, .89985, 4.1e-4, 0.0,
X .7233316, 16.92, -4.4
X },
X
X { /* mars... */
X
X 293.737334, 53.17137642, 3.107e-4, 0.0,
X 3.34218203e2, 1.8407584, 1.299e-4, -1.19e-6,
X 9.33129e-2, 9.2064e-5, 7.7e-8, 0.0,
X 1.850333, -6.75e-4, 1.26e-5, 0.0,
X 48.786442, .7709917, -1.4e-6, -5.33e-6,
X 1.5236883, 9.36, -1.52
X },
X
X { /* jupiter... */
X
X 238.049257, 8.434172183, 3.347e-4, -1.65e-6,
X 1.2720972e1, 1.6099617, 1.05627e-3, -3.43e-6,
X 4.833475e-2, 1.6418e-4, -4.676e-7, -1.7e-9,
X 1.308736, -5.6961e-3, 3.9e-6, 0.0,
X 99.443414, 1.01053, 3.5222e-4, -8.51e-6,
X 5.202561, 196.74, -9.4
X },
X
X { /* saturn... */
X
X 266.564377, 3.398638567, 3.245e-4, -5.8e-6,
X 9.1098214e1, 1.9584158, 8.2636e-4, 4.61e-6,
X 5.589232e-2, -3.455e-4, -7.28e-7, 7.4e-10,
X 2.492519, -3.9189e-3, -1.549e-5, 4e-8,
X 112.790414, .8731951, -1.5218e-4, -5.31e-6,
X 9.554747, 165.6, -8.88
X },
X
X { /* uranus... */
X
X 244.19747, 1.194065406, 3.16e-4, -6e-7,
X 1.71548692e2, 1.4844328, 2.372e-4, -6.1e-7,
X 4.63444e-2, -2.658e-5, 7.7e-8, 0.0,
X .772464, 6.253e-4, 3.95e-5, 0.0,
X 73.477111, .4986678, 1.3117e-3, 0.0,
X 19.21814, 65.8, -7.19
X },
X
X { /* neptune... */
X
X 84.457994, .6107942056, 3.205e-4, -6e-7,
X 4.6727364e1, 1.4245744, 3.9082e-4, -6.05e-7,
X 8.99704e-3, 6.33e-6, -2e-9, 0.0,
X 1.779242, -9.5436e-3, -9.1e-6, 0.0,
X 130.681389, 1.098935, 2.4987e-4, -4.718e-6,
X 30.10957, 62.2, -6.87
X },
X
X { /* pluto...(osculating 1984 jan 21) */
X
X 95.3113544, .3980332167, 0.0, 0.0,
X 224.017, 0.0, 0.0, 0.0,
X .25515, 0.0, 0.0, 0.0,
X 17.1329, 0.0, 0.0, 0.0,
X 110.191, 0.0, 0.0, 0.0,
X 39.8151, 8.2, -1.0
X }
X};
X
X/* given a modified Julian date, mjd, return the elements for the mean orbit
X * at that instant of all the major planets, together with their
X * mean daily motions in longitude, angular diameter and standard visual
X * magnitude.
X * plan[i][j] contains all the values for all the planets at mjd, such that
X * i = 0..7: mercury, venus, mars, jupiter, saturn, unranus, neptune, pluto;
X * j = 0..8: mean longitude, mean daily motion in longitude, longitude of
X * the perihelion, eccentricity, inclination, longitude of the ascending
X * node, length of the semi-major axis, angular diameter from 1 AU, and
X * the standard visual magnitude (see elements[][] comment, above).
X */
Xpelement (mjd, plan)
Xdouble mjd;
Xdouble plan[8][9];
X{
X register double *ep, *pp;
X register double t = mjd/36525.;
X double aa;
X int planet, i;
X
X for (planet = 0; planet < 8; planet++) {
X ep = elements[planet];
X pp = plan[planet];
X aa = ep[1]*t;
X pp[0] = ep[0] + 360.*(aa-(long)aa) + (ep[3]*t + ep[2])*t*t;
X range (pp, 360.);
X pp[1] = (ep[1]*9.856263e-3) + (ep[2] + ep[3])/36525;
X
X for (i = 4; i < 20; i += 4)
X pp[i/4+1] = ((ep[i+3]*t + ep[i+2])*t + ep[i+1])*t + ep[i+0];
X
X pp[6] = ep[20];
X pp[7] = ep[21];
X pp[8] = ep[22];
X }
X}
END_OF_FILE
if test 4797 -ne `wc -c <'pelement.c'`; then
echo shar: \"'pelement.c'\" unpacked with wrong size!
fi
# end of 'pelement.c'
fi
if test -f 'popup.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'popup.c'\"
else
echo shar: Extracting \"'popup.c'\" \(1724 characters\)
sed "s/^X//" >'popup.c' <<'END_OF_FILE'
X/* put up a one-line menu consisting of the given fields and let op move
X * between them with the same methods as sel_fld().
X * return index of which he picked, or -1 if hit END.
X */
X
X#include <stdio.h>
X#include "screen.h"
X
Xextern void bye();
X
X#define FLDGAP 2 /* inter-field gap */
X#define MAXFLDS 32 /* max number of fields we can handle */
X
Xstatic char pup[] = "Select: ";
X
X/* put up an array of strings on prompt line and let op pick one.
X * start with field fn.
X * N.B. we do not do much error/bounds checking.
X */
Xpopup (fields, fn, nfields)
Xchar *fields[];
Xint fn;
Xint nfields;
X{
X int fcols[MAXFLDS]; /* column to use for each field */
X int i;
X
X if (nfields > MAXFLDS)
X return (-1);
X
X again:
X /* erase the prompt line; we are going to take it over */
X c_pos (R_PROMPT, C_PROMPT);
X c_eol();
X
X /* compute starting column for each field */
X fcols[0] = sizeof(pup);
X for (i = 1; i < nfields; i++)
X fcols[i] = fcols[i-1] + strlen (fields[i-1]) + FLDGAP;
X
X /* draw each field, with comma after all but last */
X c_pos (R_PROMPT, 1);
X (void) fputs (pup, stdout);
X for (i = 0; i < nfields; i++) {
X c_pos (R_PROMPT, fcols[i]);
X printf (i < nfields-1 ? "%s," : "%s", fields[i]);
X }
X
X /* let op choose one now; begin at fn.
X */
X while (1) {
X c_pos (R_PROMPT, fcols[fn]);
X switch (read_char()) {
X case END: return (-1);
X case QUIT:
X f_prompt ("Exit ephem? (y) ");
X if (read_char() == 'y')
X bye(); /* never returns */
X goto again;
X case REDRAW: redraw_screen(2); goto again;
X case VERSION: version(); goto again;
X case '\r': return (fn);
X case 'h':
X if (--fn < 0)
X fn = nfields - 1;
X break;
X case 'l':
X if (++fn >= nfields)
X fn = 0;
X break;
X }
X }
X}
END_OF_FILE
if test 1724 -ne `wc -c <'popup.c'`; then
echo shar: \"'popup.c'\" unpacked with wrong size!
fi
# end of 'popup.c'
fi
if test -f 'precess.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'precess.c'\"
else
echo shar: Extracting \"'precess.c'\" \(1469 characters\)
sed "s/^X//" >'precess.c' <<'END_OF_FILE'
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X
X/* corrects ra and dec, both in radians, for precession from epoch 1 to epoch 2.
X * the epochs are given by their modified JDs, mjd1 and mjd2, respectively.
X * N.B. ra and dec are modifed IN PLACE.
X * TODO: find a better algorithm; this one is not even symmetric.
X */
Xprecess (mjd1, mjd2, ra, dec)
Xdouble mjd1, mjd2; /* initial and final epoch modified JDs */
Xdouble *ra, *dec; /* ra/dec for mjd1 in, for mjd2 out */
X{
X static double lastmjd1 = -10000, lastmjd2 = -10000;
X static double m, n, nyrs;
X double dra, ddec; /* ra and dec corrections */
X
X if (mjd1 != lastmjd1 || mjd2 != lastmjd2) {
X double t1, t2; /* Julian centuries of 36525 days since Jan .5 1900*/
X double m1, n1; /* "constants" for t1 */
X double m2, n2; /* "constants" for t2 */
X t1 = mjd1/36525.;
X t2 = mjd2/36525.;
X m1 = 3.07234+(1.86e-3*t1);
X n1 = 20.0468-(8.5e-3*t1);
X m2 = 3.07234+(1.86e-3*t2);
X n2 = 20.0468-(8.5e-3*t2);
X m = (m1+m2)/2; /* average m for the two epochs */
X n = (n1+n2)/2; /* average n for the two epochs */
X nyrs = (mjd2-mjd1)/365.2425;
X lastmjd1 = mjd1;
X lastmjd2 = mjd2;
X }
X
X dra = (m+(n*sin(*ra)*tan(*dec)/15))*7.272205e-5*nyrs;
X ddec = n*cos(*ra)*4.848137e-6*nyrs;
X *ra += dra;
X *dec += ddec;
X /* added by ECD */
X if (*dec > PI/2) {
X *dec = PI - *dec;
X *ra += PI;
X } else if (*dec < -PI/2) {
X *dec = -PI - *dec;
X *ra += PI;
X }
X range (ra, 2*PI);
X}
END_OF_FILE
if test 1469 -ne `wc -c <'precess.c'`; then
echo shar: \"'precess.c'\" unpacked with wrong size!
fi
# end of 'precess.c'
fi
if test -f 'reduce.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'reduce.c'\"
else
echo shar: Extracting \"'reduce.c'\" \(1691 characters\)
sed "s/^X//" >'reduce.c' <<'END_OF_FILE'
X#include <math.h>
X#include "astro.h"
X
X/* convert those orbital elements that change from epoch mjd0 to epoch mjd.
X */
Xreduce_elements (mjd0, mjd, inc0, ap0, om0, inc, ap, om)
Xdouble mjd0; /* initial epoch */
Xdouble mjd; /* desired epoch */
Xdouble inc0; /* initial inclination, rads */
Xdouble ap0; /* initial argument of perihelion, as an mjd */
Xdouble om0; /* initial long of ascending node, rads */
Xdouble *inc; /* desired inclination, rads */
Xdouble *ap; /* desired epoch of perihelion, as an mjd */
Xdouble *om; /* desired long of ascending node, rads */
X{
X double t0, t1;
X double tt, tt2, t02, tt3;
X double eta, th, th0;
X double a, b;
X double dap;
X double cinc, sinc;
X double ot, sot, cot, ot1;
X double seta, ceta;
X
X t0 = mjd0/365250.0;
X t1 = mjd/365250.0;
X
X tt = t1-t0;
X tt2 = tt*tt;
X t02 = t0*t0;
X tt3 = tt*tt2;
X eta = (471.07-6.75*t0+.57*t02)*tt+(.57*t0-3.37)*tt2+.05*tt3;
X th0 = 32869.0*t0+56*t02-(8694+55*t0)*tt+3*tt2;
X eta = degrad(eta/3600.0);
X th0 = degrad((th0/3600.0)+173.950833);
X th = (50256.41+222.29*t0+.26*t02)*tt+(111.15+.26*t0)*tt2+.1*tt3;
X th = th0+degrad(th/3600.0);
X cinc = cos(inc0);
X sinc = sin(inc0);
X ot = om0-th0;
X sot = sin(ot);
X cot = cos(ot);
X seta = sin(eta);
X ceta = cos(eta);
X a = sinc*sot;
X b = ceta*sinc*cot-seta*cinc;
X ot1 = atan(a/b);
X if (b<0) ot1 += PI;
X b = sinc*ceta-cinc*seta*cot;
X a = -1*seta*sot;
X dap = atan(a/b);
X if (b<0) dap += PI;
X
X *ap = ap0+dap;
X range (ap, 2*PI);
X *om = ot1+th;
X range (om, 2*PI);
X
X if (inc0<.175)
X *inc = asin(a/sin(dap));
X else
X *inc = 1.570796327-asin((cinc*ceta)+(sinc*seta*cot));
X}
END_OF_FILE
if test 1691 -ne `wc -c <'reduce.c'`; then
echo shar: \"'reduce.c'\" unpacked with wrong size!
fi
# end of 'reduce.c'
fi
if test -f 'refract.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'refract.c'\"
else
echo shar: Extracting \"'refract.c'\" \(1857 characters\)
sed "s/^X//" >'refract.c' <<'END_OF_FILE'
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X
X/* correct the true altitude, ta, for refraction to the apparent altitude, aa,
X * each in radians, given the local atmospheric pressure, pr, in mbars, and
X * the temperature, tr, in degrees C.
X */
Xrefract (pr, tr, ta, aa)
Xdouble pr, tr;
Xdouble ta;
Xdouble *aa;
X{
X double r; /* refraction correction*/
X
X if (ta >= degrad(15.)) {
X /* model for altitudes at least 15 degrees above horizon */
X r = 7.888888e-5*pr/((273+tr)*tan(ta));
X } else if (ta > degrad(-5.)) {
X /* hairier model for altitudes at least -5 and below 15 degrees */
X double a, b, tadeg = raddeg(ta);
X a = ((2e-5*tadeg+1.96e-2)*tadeg+1.594e-1)*pr;
X b = (273+tr)*((8.45e-2*tadeg+5.05e-1)*tadeg+1);
X r = degrad(a/b);
X } else {
X /* do nothing if more than 5 degrees below horizon.
X */
X r = 0;
X }
X
X *aa = ta + r;
X}
X
X/* correct the apparent altitude, aa, for refraction to the true altitude, ta,
X * each in radians, given the local atmospheric pressure, pr, in mbars, and
X * the temperature, tr, in degrees C.
X */
Xunrefract (pr, tr, aa, ta)
Xdouble pr, tr;
Xdouble aa;
Xdouble *ta;
X{
X double err;
X double appar;
X double true;
X
X /* iterative solution: search for the true that refracts to the
X * given apparent.
X * since refract() is discontinuous at -5 degrees, there is a range
X * of apparent altitudes between about -4.5 and -5 degrees that are
X * not invertable (the graph of ap vs. true has a vertical step at
X * true = -5). thus, the iteration just oscillates if it gets into
X * this region. if this happens the iteration is forced to abort.
X * of course, this makes unrefract() discontinuous too.
X */
X true = aa;
X do {
X refract (pr, tr, true, &appar);
X err = appar - aa;
X true -= err;
X } while (fabs(err) >= 1e-6 && true > degrad(-5));
X
X *ta = true;
X}
END_OF_FILE
if test 1857 -ne `wc -c <'refract.c'`; then
echo shar: \"'refract.c'\" unpacked with wrong size!
fi
# end of 'refract.c'
fi
if test -f 'riset.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'riset.c'\"
else
echo shar: Extracting \"'riset.c'\" \(4591 characters\)
sed "s/^X//" >'riset.c' <<'END_OF_FILE'
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X
X/* given the true geocentric ra and dec of an object, the observer's latitude,
X * lat, and a horizon displacement correction, dis, all in radians, find the
X * local sidereal times and azimuths of rising and setting, lstr/s
X * and azr/s, also all in radians, respectively.
X * dis is the vertical displacement from the true position of the horizon. it
X * is positive if the apparent position is higher than the true position.
X * said another way, it is positive if the shift causes the object to spend
X * longer above the horizon. for example, atmospheric refraction is typically
X * assumed to produce a vertical shift of 34 arc minutes at the horizon; dis
X * would then take on the value +9.89e-3 (radians). On the other hand, if
X * your horizon has hills such that your apparent horizon is, say, 1 degree
X * above sea level, you would allow for this by setting dis to -1.75e-2
X * (radians).
X *
X * algorithm:
X * the situation is described by two spherical triangles with two equal angles
X * (the right horizon intercepts, and the common horizon transverse):
X * given lat, d(=d1+d2), and dis find z(=z1+z2) and rho, where /| eq pole
X * lat = latitude, / |
X * dis = horizon displacement (>0 is below ideal) / rho|
X * d = angle from pole = PI/2 - declination / |
X * z = azimuth east of north / |
X * rho = polar rotation from down = PI - hour angle / |
X * solve simultaneous equations for d1 and d2: / |
X * 1) cos(d) = cos(d1+d2) / d2 | lat
X * = cos(d1)cos(d2) - sin(d1)sin(d2) / |
X * 2) sin(d2) = sin(lat)sin(d1)/sin(dis) / |
X * then can solve for z1, z2 and rho, taking / |
X * care to preserve quadrant information. / -|
X * z1 / z2 | |
X * ideal horizon ------------/--------------------|
X * | | / N
X * -| / d1
X * dis | /
X * |/
X * apparent horizon ---------------------------------
X *
X * note that when lat=0 this all breaks down (because d2 and z2 degenerate to 0)
X * but fortunately then we can solve for z and rho directly.
X *
X * status: 0: normal; 1: never rises; -1: circumpolar; 2: trouble.
X */
Xriset (ra, dec, lat, dis, lstr, lsts, azr, azs, status)
Xdouble ra, dec;
Xdouble lat, dis;
Xdouble *lstr, *lsts;
Xdouble *azr, *azs;
Xint *status;
X{
X#define EPS (1e-6) /* math rounding fudge - always the way, eh? */
X double d; /* angle from pole */
X double h; /* hour angle */
X double crho; /* cos hour-angle complement */
X int shemi; /* flag for southern hemisphere reflection */
X
X d = PI/2 - dec;
X
X /* reflect if in southern hemisphere.
X * (then reflect azimuth back after computation.)
X */
X if (shemi = lat < 0) {
X lat = -lat;
X d = PI - d;
X }
X
X /* do the easy ones (and avoid violated assumptions) if d arc never
X * meets horizon.
X */
X if (d <= lat + dis + EPS) {
X *status = -1; /* never sets */
X return;
X }
X if (d >= PI - lat + dis - EPS) {
X *status = 1; /* never rises */
X return;
X }
X
X /* find rising azimuth and cosine of hour-angle complement */
X if (lat > EPS) {
X double d2, d1; /* polr arc to ideal hzn, and corrctn for apparent */
X double z2, z1; /* azimuth to ideal horizon, and " */
X double a; /* intermediate temp */
X double sdis, slat, clat, cz2, cd2; /* trig temps */
X sdis = sin(dis);
X slat = sin(lat);
X a = sdis*sdis + slat*slat + 2*cos(d)*sdis*slat;
X if (a <= 0) {
X *status = 2; /* can't happen - hah! */
X return;
X }
X d1 = asin (sin(d) * sdis / sqrt(a));
X d2 = d - d1;
X cd2 = cos(d2);
X clat = cos(lat);
X cz2 = cd2/clat;
X z2 = acos (cz2);
X z1 = acos (cos(d1)/cos(dis));
X if (dis < 0)
X z1 = -z1;
X *azr = z1 + z2;
X range (azr, PI);
X crho = (cz2 - cd2*clat)/(sin(d2)*slat);
X } else {
X *azr = acos (cos(d)/cos(dis));
X crho = sin(dis)/sin(d);
X }
X
X if (shemi)
X *azr = PI - *azr;
X *azs = 2*PI - *azr;
X
X /* find hour angle */
X h = PI - acos (crho);
X *lstr = radhr(ra-h);
X *lsts = radhr(ra+h);
X range (lstr, 24.0);
X range (lsts, 24.0);
X
X *status = 0;
X}
END_OF_FILE
if test 4591 -ne `wc -c <'riset.c'`; then
echo shar: \"'riset.c'\" unpacked with wrong size!
fi
# end of 'riset.c'
fi
if test -f 'sex_dec.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'sex_dec.c'\"
else
echo shar: Extracting \"'sex_dec.c'\" \(1194 characters\)
sed "s/^X//" >'sex_dec.c' <<'END_OF_FILE'
X/* given hours (or degrees), hd, minutes, m, and seconds, s,
X * return decimal hours (or degrees), *d.
X * in the case of hours (angles) < 0, only the first non-zero element should
X * be negative.
X */
Xsex_dec (hd, m, s, d)
Xint hd, m, s;
Xdouble *d;
X{
X int sign = 1;
X
X if (hd < 0) {
X sign = -1;
X hd = -hd;
X } else if (m < 0) {
X sign = -1;
X m = -m;
X } else if (s < 0) {
X sign = -1;
X s = -s;
X }
X
X *d = (((double)s/60.0 + (double)m)/60.0 + (double)hd) * sign;
X}
X
X/* given decimal hours (or degrees), d.
X * return nearest hours (or degrees), *hd, minutes, *m, and seconds, *s,
X * each always non-negative; *isneg is set to 1 if d is < 0, else to 0.
X */
Xdec_sex (d, hd, m, s, isneg)
Xdouble d;
Xint *hd, *m, *s, *isneg;
X{
X double min;
X
X if (d < 0) {
X *isneg = 1;
X d = -d;
X } else
X *isneg = 0;
X
X *hd = (int)d;
X min = (d - *hd)*60.;
X *m = (int)min;
X *s = (int)((min - *m)*60. + 0.5);
X
X if (*s == 60) {
X if ((*m += 1) == 60) {
X *hd += 1;
X *m = 0;
X }
X *s = 0;
X }
X /* no negative 0's */
X if (*hd == 0 && *m == 0 && *s == 0)
X *isneg = 0;
X}
X
X/* insure 0 <= *v < r.
X */
Xrange (v, r)
Xdouble *v, r;
X{
X while (*v < 0) *v += r;
X while (*v >= r) *v -= r;
X}
END_OF_FILE
if test 1194 -ne `wc -c <'sex_dec.c'`; then
echo shar: \"'sex_dec.c'\" unpacked with wrong size!
fi
# end of 'sex_dec.c'
fi
if test -f 'sun.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'sun.c'\"
else
echo shar: Extracting \"'sun.c'\" \(1760 characters\)
sed "s/^X//" >'sun.c' <<'END_OF_FILE'
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X
X/* given the modified JD, mjd, return the true geocentric ecliptic longitude
X * of the sun for the mean equinox of the date, *lsn, in radians, and the
X * sun-earth distance, *rsn, in AU. (the true ecliptic latitude is never more
X * than 1.2 arc seconds and so may be taken to be a constant 0.)
X * if the APPARENT ecliptic longitude is required, correct the longitude for
X * nutation to the true equinox of date and for aberration (light travel time,
X * approximately -9.27e7/186000/(3600*24*365)*2*pi = -9.93e-5 radians).
X */
Xsunpos (mjd, lsn, rsn)
Xdouble mjd;
Xdouble *lsn, *rsn;
X{
X double t, t2;
X double ls, ms; /* mean longitude and mean anomoay */
X double s, nu, ea; /* eccentricity, true anomaly, eccentric anomaly */
X double a, b, a1, b1, c1, d1, e1, h1, dl, dr;
X
X t = mjd/36525.;
X t2 = t*t;
X a = 100.0021359*t;
X b = 360.*(a-(long)a);
X ls = 279.69668+.0003025*t2+b;
X a = 99.99736042000039*t;
X b = 360*(a-(long)a);
X ms = 358.47583-(.00015+.0000033*t)*t2+b;
X s = .016751-.0000418*t-1.26e-07*t2;
X anomaly (degrad(ms), s, &nu, &ea);
X a = 62.55209472000015*t;
X b = 360*(a-(long)a);
X a1 = degrad(153.23+b);
X a = 125.1041894*t;
X b = 360*(a-(long)a);
X b1 = degrad(216.57+b);
X a = 91.56766028*t;
X b = 360*(a-(long)a);
X c1 = degrad(312.69+b);
X a = 1236.853095*t;
X b = 360*(a-(long)a);
X d1 = degrad(350.74-.00144*t2+b);
X e1 = degrad(231.19+20.2*t);
X a = 183.1353208*t;
X b = 360*(a-(long)a);
X h1 = degrad(353.4+b);
X dl = .00134*cos(a1)+.00154*cos(b1)+.002*cos(c1)+.00179*sin(d1)+
X .00178*sin(e1);
X dr = 5.43e-06*sin(a1)+1.575e-05*sin(b1)+1.627e-05*sin(c1)+
X 3.076e-05*cos(d1)+9.27e-06*sin(h1);
X *lsn = nu+degrad(ls-ms+dl);
X *rsn = 1.0000002*(1-s*cos(ea))+dr;
X range (lsn, 2*PI);
X}
END_OF_FILE
if test 1760 -ne `wc -c <'sun.c'`; then
echo shar: \"'sun.c'\" unpacked with wrong size!
fi
# end of 'sun.c'
fi
if test -f 'time.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'time.c'\"
else
echo shar: Extracting \"'time.c'\" \(2869 characters\)
sed "s/^X//" >'time.c' <<'END_OF_FILE'
X/* get the time from the os.
X *
X * here are two methods I was able to verify; pick one for your system and
X * define exactly one of TZA or TZB:
X * TZA works on our ibm-pc/turbo-c and at&t systems,
X * TZB works on our 4.2 BSD vax.
X *
X * I'm told that on Sun OS 4.0.3 (BSD 4.3?) and Apollo SR 10.1 TZB works if
X * you use <sys/time.h> in place of <time.h>.
X *
X * On VMS, you DON'T want to define EITHER TZA nor TZB since it can't handle
X * time zones, period. time_fromsys() will detect that fact based on gmtime()
X * returning 0.
X */
X
X#define TZA
X
X#ifdef VMS
X#undef TZA
X#undef TZB
X#endif
X
X#include <stdio.h>
X#include <time.h>
X
X#include "astro.h"
X#include "circum.h"
X
Xextern char *strncpy();
Xextern long time();
X
Xstatic long c0;
Xstatic double mjd0;
X
X/* save current mjd and corresponding system clock for use by inc_mjd().
X * this establishes the base correspondence between the mjd and system clock.
X */
Xset_t0 (np)
XNow *np;
X{
X mjd0 = mjd;
X (void) time (&c0);
X}
X
X/* fill in n_mjd/tz/tznm from system clock.
X */
Xtime_fromsys (np)
XNow *np;
X{
X extern struct tm *gmtime(), *localtime();
X struct tm *tp;
X long c;
X double day, hr;
X
X (void) time (&c);
X
X tp = gmtime (&c);
X if (tp) {
X cal_mjd (tp->tm_mon+1, (double)tp->tm_mday, tp->tm_year+1900, &day);
X sex_dec (tp->tm_hour, tp->tm_min, tp->tm_sec, &hr);
X mjd = day + hr/24.0;
X tp = localtime (&c);
X settzstuff (tp->tm_isdst ? 1 : 0, np);
X } else {
X /* if gmtime() doesn't work, we assume the timezone stuff won't
X * either, so we just use what it is and leave it alone. Some
X * systems (like VMS) do not know about time zones, so this is the
X * best guess in that case.
X */
X tp = localtime (&c);
X cal_mjd (tp->tm_mon+1, (double)tp->tm_mday, tp->tm_year+1900, &day);
X sex_dec (tp->tm_hour, tp->tm_min, tp->tm_sec, &hr);
X mjd = day + hr/24.0 + tz/24.0;
X }
X}
X
X/* given whether dst is now in effect (must be strictly 0 or 1), fill in
X * tzname and tz within np.
X */
Xstatic
Xsettzstuff (dst, np)
Xint dst;
XNow *np;
X{
X#ifdef TZA
X extern long timezone;
X extern char *tzname[2];
X
X tzset();
X tz = timezone/3600;
X if (dst)
X tz -= 1.0;
X (void) strncpy (tznm, tzname[dst], sizeof(tznm)-1);
X#endif
X#ifdef TZB
X extern char *timezone();
X struct timeval timev;
X struct timezone timez;
X
X gettimeofday (&timev, &timez);
X tz = timez.tz_minuteswest/60;
X if (dst)
X tz -= 1.0;
X (void) strncpy (tznm, timezone(timez.tz_minuteswest, dst),
X sizeof(tznm)-1);
X#endif
X tznm[sizeof(tznm)-1] = '\0'; /* insure string is terminated */
X}
X
Xinc_mjd (np, inc)
XNow *np;
Xdouble inc;
X{
X if (inc == RTC) {
X long c;
X (void) time (&c);
X mjd = mjd0 + (c - c0)/SPD;
X } else
X mjd += inc/24.0;
X
X /* round to nearest whole second.
X * without this, you can get fractional days so close to .5 but
X * not quite there that mjd_hr() can return 24.0
X */
X rnd_second (&mjd);
X}
END_OF_FILE
if test 2869 -ne `wc -c <'time.c'`; then
echo shar: \"'time.c'\" unpacked with wrong size!
fi
# end of 'time.c'
fi
if test -f 'utc_gst.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'utc_gst.c'\"
else
echo shar: Extracting \"'utc_gst.c'\" \(1189 characters\)
sed "s/^X//" >'utc_gst.c' <<'END_OF_FILE'
X#include "astro.h"
X
X/* given a modified julian date, mjd, and a universally coordinated time, utc,
X * return greenwich mean siderial time, *gst.
X */
Xutc_gst (mjd, utc, gst)
Xdouble mjd;
Xdouble utc;
Xdouble *gst;
X{
X double tnaught();
X static double lastmjd = -10000;
X static double t0;
X
X if (mjd != lastmjd) {
X t0 = tnaught (mjd);
X lastmjd = mjd;
X }
X *gst = (1.0/SIDRATE)*utc + t0;
X range (gst, 24.0);
X}
X
X/* given a modified julian date, mjd, and a greenwich mean siderial time, gst,
X * return universally coordinated time, *utc.
X */
Xgst_utc (mjd, gst, utc)
Xdouble mjd;
Xdouble gst;
Xdouble *utc;
X{
X double tnaught();
X static double lastmjd = -10000;
X static double t0;
X
X if (mjd != lastmjd) {
X t0 = tnaught (mjd);
X range (&t0, 24.0);
X lastmjd = mjd;
X }
X *utc = gst - t0;
X range (utc, 24.0);
X *utc *= SIDRATE;
X}
X
Xstatic double
Xtnaught (mjd)
Xdouble mjd; /* julian days since 1900 jan 0.5 */
X{
X double dmjd;
X int m, y;
X double d;
X double t, t0;
X
X mjd_cal (mjd, &m, &d, &y);
X cal_mjd (1, 0., y, &dmjd);
X t = dmjd/36525;
X t0 = 6.57098e-2 * (mjd - dmjd) -
X (24 - (6.6460656 + (5.1262e-2 + (t * 2.581e-5))*t) -
X (2400 * (t - (((double)y - 1900)/100))));
X return (t0);
X}
END_OF_FILE
if test 1189 -ne `wc -c <'utc_gst.c'`; then
echo shar: \"'utc_gst.c'\" unpacked with wrong size!
fi
# end of 'utc_gst.c'
fi
echo shar: End of archive 1 \(of 6\).
cp /dev/null ark1isdone
MISSING=""
for I in 1 2 3 4 5 6 ; do
if test ! -f ark${I}isdone ; then
MISSING="${MISSING} ${I}"
fi
done
if test "${MISSING}" = "" ; then
echo You have unpacked all 6 archives.
rm -f ark[1-9]isdone
else
echo You still need to unpack the following archives:
echo " " ${MISSING}
fi
## End of shell archive.
exit 0