Source Code 1994 March/Source_Code_CD-ROM_Walnut_Creek_March_1994.iso/compsrcs/x/volume19/xephem/part21

Text File | 1993-05-15 | 57.9 KB | 1,918 lines

Newsgroups: comp.sources.x
From: [e-mail] (Elwood Downey)
Subject: v19i109:  xephem - astronomical ephemeris program, Part21/21
Message-ID: <[e-mail]>
X-Md4-Signature: c2ebeab684497b5418ce2dfbd923d04e
Date: Mon, 10 May 1993 22:14:58 GMT
Approved: [e-mail]

Submitted-by: [e-mail] (Elwood Downey)
Posting-number: Volume 19, Issue 109
Archive-name: xephem/part21
Environment: X11r4, OSF/Motif
Supersedes: xephem: Volume 16, Issue 112-134

#! /bin/sh
# This is a shell archive.  Remove anything before this line, then feed it
# into a shell via "sh file" or similar.  To overwrite existing files,
# type "sh file -c".
# The tool that generated this appeared in the comp.sources.unix newsgroup;
# send mail to [e-mail] if you want that tool.
# Contents:  Imakefile Makefile.simple Manifest aa_hadec.c astro.h
#   eq_ecl.c marsmap.c moon.c moonnf.c msgmenu.c nutation.c obliq.c
#   patchlevel.h pelement.c reduce.c refract.c riset.c sun.c utc_gst.c
# Wrapped by chris@nova on Mon May 10 16:41:55 1993
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
echo If this archive is complete, you will see the following message:
echo '          "shar: End of archive 21 (of 21)."'
if test -f 'Imakefile' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'Imakefile'\"
else
  echo shar: Extracting \"'Imakefile'\" \(1456 characters\)
  sed "s/^X//" >'Imakefile' <<'END_OF_FILE'
XSRCS =    aa_hadec.c anomaly.c cal_mjd.c calmenu.c circum.c comet.c compiler.c \
X    constel.c datamenu.c db.c dbmenu.c earthmap.c earthmenu.c eq_ecl.c \
X    formats.c helpmenu.c jupmenu.c listmenu.c mainmenu.c marsmap.c \
X    marsmenu.c misc.c moon.c moonmenu.c moonnf.c msgmenu.c nutation.c \
X    objmenu.c obliq.c parallax.c pelement.c plans.c plot_aux.c plotmenu.c \
X    precess.c preferences.c query.c reduce.c refract.c riset.c \
X    riset_cir.c satmenu.c sex_dec.c skyfiltmenu.c skyviewmenu.c \
X    solsysmenu.c srchmenu.c sun.c time.c utc_gst.c versionmenu.c xephem.c
X
XOBJS =    aa_hadec.o anomaly.o cal_mjd.o calmenu.o circum.o comet.o compiler.o \
X    constel.o datamenu.o db.o dbmenu.o earthmap.o earthmenu.o eq_ecl.o \
X    formats.o helpmenu.o jupmenu.o listmenu.o mainmenu.o marsmap.o \
X    marsmenu.o misc.o moon.o moonmenu.o moonnf.o msgmenu.o nutation.o \
X    objmenu.o obliq.o parallax.o pelement.o plans.o plot_aux.o plotmenu.o \
X    precess.o preferences.o query.o reduce.o refract.o riset.o \
X    riset_cir.o satmenu.o sex_dec.o skyfiltmenu.o skyviewmenu.o \
X    solsysmenu.o srchmenu.o sun.o time.o utc_gst.o versionmenu.o xephem.o
X
X       DEPXMLIB = $(USRLIBDIR)/libXm.a
X          XMLIB = -lXm
XLOCAL_LIBRARIES = $(XMLIB) $(XTOOLLIB) $(XLIB)
X        DEPLIBS = $(DEPXMLIB) $(DEPXTOOLLIB) $(DEPXLIB)
X  SYS_LIBRARIES = -lm
X       PROGRAMS = xephem
X
XComplexProgramTarget(xephem)
X
X# especially worth having this dependency so the displayed version is current.
Xversionmenu.o:    patchlevel.h
END_OF_FILE
  if test 1456 -ne `wc -c <'Imakefile'`; then
    echo shar: \"'Imakefile'\" unpacked with wrong size!
  fi
  # end of 'Imakefile'
fi
if test -f 'Makefile.simple' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'Makefile.simple'\"
else
  echo shar: Extracting \"'Makefile.simple'\" \(1493 characters\)
  sed "s/^X//" >'Makefile.simple' <<'END_OF_FILE'
X# Simple Makefile for xephem v2.4b.
X# We include sample compile and link flags for a few popular systems.
X
X# "stock" X systems
XCLDFLAGS =
XCFLAGS = $(CLDFLAGS) -O -D_NO_PROTO
XLDFLAGS = $(CLDFLAGS)
XLIBS = -lXm -lXt -lX11 -lm
X
X# SVR4
X# CLDFLAGS =
X# CFLAGS = $(CLDFLAGS) -O
X# LDFLAGS = $(CLDFLAGS)
X# LIBS = -lXm -lXt -lX11 -lsocket -lnsl -lc -lm /usr/ucblib/libucb.a
X
X# SVR3.2 with Metrolink X
X# CLDFLAGS =
X# CFLAGS = $(CLDFLAGS) -D_NO_PROTO -DSYSV -O
X# LDFLAGS = $(CLDFLAGS) -L/usr/lib/X11/libs
X# LIBS = -lXm -lXt -lX11 -lm -lpt -lsocket -lnet -lnsl_s -lc -lgnu -lPW
X
X# try the following CFLAGS and LIBS if set WANT_EDITRES in xephem.c
X# CLDFLAGS =
X# CFLAGS = $(CLDFLAGS) -DWANT_EDITRES -O
X# LDFLAGS = $(CLDFLAGS)
X# LIBS= -lXm -lXt -lXmu -lXext -lX11 -lm
X
XOBJS =    aa_hadec.o anomaly.o cal_mjd.o calmenu.o circum.o comet.o compiler.o \
X    constel.o datamenu.o db.o dbmenu.o earthmap.o earthmenu.o eq_ecl.o \
X    formats.o helpmenu.o jupmenu.o listmenu.o mainmenu.o marsmap.o \
X    marsmenu.o misc.o moon.o moonmenu.o moonnf.o msgmenu.o nutation.o \
X    objmenu.o obliq.o parallax.o pelement.o plans.o plot_aux.o plotmenu.o \
X    precess.o preferences.o query.o reduce.o refract.o riset.o \
X    riset_cir.o satmenu.o sex_dec.o skyfiltmenu.o skyviewmenu.o \
X    solsysmenu.o srchmenu.o sun.o time.o utc_gst.o versionmenu.o xephem.o
X
X.PRECIOUS:    xephem
X
Xxephem: $(OBJS)
X    $(CC) $(LDFLAGS) -o $@ $(OBJS) $(LIBS)
X
X# especially worth having this dependency so the displayed version is current.
Xversionmenu.o:    patchlevel.h
END_OF_FILE
  if test 1493 -ne `wc -c <'Makefile.simple'`; then
    echo shar: \"'Makefile.simple'\" unpacked with wrong size!
  fi
  # end of 'Makefile.simple'
fi
if test -f 'Manifest' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'Manifest'\"
else
  echo shar: Extracting \"'Manifest'\" \(3608 characters\)
  sed "s/^X//" >'Manifest' <<'END_OF_FILE'
XCopyright    notice of ownership, permissions and limitations
XImakefile    sample Imakefile for xmkmf
XMakefile.simple    simple sample make file.
XManifest    this file.
XREADME        check here for hints before building.
XXEphem.ad    defaults file (basically the same as fallback ones in xephem.c)
Xephem.db    sample database file.
Xephem.plt    sample plot file (an analemma)
Xxephem.hlp    help text file.
Xxephem.man    manual page, written with the nroff -man macros.
X
Xaa_hadec.c      convert between alt/az and hour angle/dec.
Xanomaly.c       compute anomaly.
Xastro.h         unit conversion macros and planet defines.
Xcal_mjd.c       converters to and from modified julian date.
Xcalmenu.c    control the calendar on the main menu.
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.
Xconstel.c       handles determing and printing constellation info.
Xdatamenu.c    manage the general data menu.
Xdb.c        read and process the database files, and maintain in-memory db.
Xdbmenu.c    manage the database menu.
Xearthmap.c    coordinates for the simple earth map.
Xearthmenu.c    manage the sunlit earth menu.
Xeq_ecl.c        convert between equitorial and eclipitic coords.
Xformats.c       basic date, time, prompts, etc formats.
Xhelpmenu.c    manage the help menu and supporting text file.
Xjupmenu.c    manage the jupiter moon map menu.
Xlistmenu.c    manage the listing control menu.
Xmainmenu.c    manage the main menu.
Xmap.h        typedefs useful for maps consisting of sets of lines.
Xmarsmap.c    coordincates for the simple mars map.
Xmarsmenu.c    manage the mars central meridian longitude menu.
Xmisc.c        everyone needs one of these!
Xmoon.c          compute moon position.
Xmoonmenu.c    manage the moon display menu.
Xmoonnf.c        compute new and full moon dates.
Xmsgmenu.c    manage the message alert mechanism.
Xnutation.c      compute nutation correction.
Xobjmenu.c    manage the objx/y menu.
Xobliq.c         compute obliquity.
Xparallax.c      functions to compute earth rim parallax correction.
Xpatchlevel.h    defines PATCHLEVEL which establishes the official release level.
Xpelement.c      basic planet position polynomial coefficients.
Xplans.c         use polynomials to find planet location at any certain time.
Xplot_aux.c    manage the actual plotting menu and read supporting files.
Xplotmenu.c    manage the plot control menu and write supporting files.
Xprecess.c       compute precession correction.
Xpreferences.c    code to maintain the setting of the user preference settings.
Xpreferences.h    header file for use by modules setting or getting preferences.
Xquery.c        general purpose menu query tool.
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_cir.c     iteratively solve for local rise/set times of moving objects.
Xsatmenu.c    manage the Saturn menu.
Xsex_dec.c       convert between sexagesimal and decimal notation.
Xskyfiltmenu.c    manage the sky view type filter menu.
Xskyviewmenu.c    manage the circular sky view menu.
Xsmallfm.xbm    small moon bitmap.
Xsolsysmenu.c    manage the solar system menu.
Xsrchmenu.c    manage the search control menu.
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
Xversionmenu.c   current version notice, and revision history comments.
Xxephem.c    main() and misc minor support utility functions.
END_OF_FILE
  if test 3608 -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'\" \(2283 characters\)
  sed "s/^X//" >'aa_hadec.c' <<'END_OF_FILE'
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X
X#if defined(__STDC__) || defined(__cplusplus)
X#define P_(s) s
X#else
X#define P_(s) ()
X#endif
X
Xvoid aa_hadec P_((double lat, double alt, double az, double *ha, double *dec));
Xvoid hadec_aa P_((double lat, double ha, double dec, double *alt, double *az));
Xstatic void aaha_aux P_((double lat, double x, double y, double *p, double *q));
X
X#undef P_
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 */
Xvoid
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 */
Xvoid
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 void
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 2283 -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 'astro.h' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'astro.h'\"
else
  echo shar: Extracting \"'astro.h'\" \(885 characters\)
  sed "s/^X//" >'astro.h' <<'END_OF_FILE'
X/* defined derived strictly from the Duffett-Smith functions.
X */
X
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 coincide with usage in pelement.c and plans.c.
X * N.B. only the first 8 are valid for use with plans().
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
X
X/* a few more handy ones */
X#define    SUN    8
X#define    MOON    9
X#define    OBJX    10
X#define    OBJY    11
X#define    NOBJ    12    /* total number of basic objects */
END_OF_FILE
  if test 885 -ne `wc -c <'astro.h'`; then
    echo shar: \"'astro.h'\" unpacked with wrong size!
  fi
  # end of 'astro.h'
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'\" \(2433 characters\)
  sed "s/^X//" >'eq_ecl.c' <<'END_OF_FILE'
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X
X#if defined(__STDC__) || defined(__cplusplus)
X#define P_(s) s
X#else
X#define P_(s) ()
X#endif
X
Xextern void nutation P_((double Mjd, double *deps, double *dpsi));
Xextern void obliquity P_((double Mjd, double *eps));
Xextern void range P_((double *v, double r));
X
Xvoid eq_ecl P_((double mjd, double ra, double dec, double *lat, double *lng));
Xvoid ecl_eq P_((double mjd, double lat, double lng, double *ra, double *dec));
Xstatic void ecleq_aux P_((int sw, double mjd, double x, double y, double *p, double *q));
X
X#undef P_
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 */
Xvoid
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 */
Xvoid
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 void
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 2433 -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 'marsmap.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'marsmap.c'\"
else
  echo shar: Extracting \"'marsmap.c'\" \(1639 characters\)
  sed "s/^X//" >'marsmap.c' <<'END_OF_FILE'
X#include "map.h"
X
Xstatic MCoord mc0[] = {
X    {325,-11}, {9,-10}, {10,5}, {5,2}, {0,0}, {352,2}, {348,-5}, {330,6}
X};
Xstatic MCoord mc1[] = {
X    {160,-40}, {200,-40}, {200,-35}, {240,-15}, {235,-4}, {195,-32}, {170,-35}
X};
Xstatic MCoord mc2[] = {
X    {220,-25},{228,-35},{240,-35},{250,-26},{260,-30},{280,-25},{282,-11},
X    {280,-3},{260,-7},{245,-20},{240,-15}
X};
Xstatic MCoord mc3[] = {
X    {281,5}, {285,1}, {298,2}, {290,22}, {283,22}
X};
Xstatic MCoord mc4[] = {
X    {120,-32}, {140,-35}, {150,-33}, {165,-25}, {165,-22}, {162,-23},
X    {142,-31}, {138,-30}
X};
Xstatic MCoord mc5[] = {
X    {40,21}, {52,30}, {55,20}, {61,19}, {61,26}, {68,29}, {68,30}, {42,50},
X    {22,48}, {10,40}, {15,35}, {20,25}
X};
Xstatic MCoord mc6[] = {
X    {72,-28}, {82,-29}, {90,-34}, {100,-22}, {90,-17}, {80,-20}, {72,-20}
X};
Xstatic MCoord mc7[] = {
X    {58,-28}, {68,-24}, {64,-19}, {56,-20}
X};
Xstatic MCoord mc8[] = {
X    {20,-31}, {30,-35}, {55,-26}, {53,-20}, {30,-20}
X};
Xstatic MCoord mc9[] = {
X    {30,-20}, {53,-20}, {65,-15}, {57,-6}, {38,-5}, {29,-16}
X};
Xstatic MCoord mc10[] = {
X    {5,-26}, {38,-8}, {21,0}
X};
X
X#define    ASIZ(a)        (sizeof(a)/sizeof(a[0]))
X
XMRegion mreg[] = {
X    {"Meridiani Sinus",        mc0, ASIZ(mc0)},
X    {"Cimmerium Mare",        mc1, ASIZ(mc1)},
X    {"Tyrrhenum Mare",        mc2, ASIZ(mc2)},
X    {"Syrtis Major",        mc3, ASIZ(mc3)},
X    {"Sirenum Mare",        mc4, ASIZ(mc4)},
X    {"Niliacus Lacus",        mc5, ASIZ(mc5)},
X    {"Solis Lacus",        mc6, ASIZ(mc6)},
X    {"Nectar",            mc7, ASIZ(mc7)},
X    {"Erythraeum Mare",        mc8, ASIZ(mc8)},
X    {"Aurorae Sinus",        mc9, ASIZ(mc9)},
X    {"Margaritifer Sinus",    mc10, ASIZ(mc10)},
X};
X
Xint nmreg = ASIZ(mreg);
END_OF_FILE
  if test 1639 -ne `wc -c <'marsmap.c'`; then
    echo shar: \"'marsmap.c'\" unpacked with wrong size!
  fi
  # end of 'marsmap.c'
fi
if test -f 'moon.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'moon.c'\"
else
  echo shar: Extracting \"'moon.c'\" \(5365 characters\)
  sed "s/^X//" >'moon.c' <<'END_OF_FILE'
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X
X#if defined(__STDC__) || defined(__cplusplus)
X#define P_(s) s
X#else
X#define P_(s) ()
X#endif
X
Xextern void range P_((double *v, double r));
X
Xvoid moon P_((double mjd, double *lam, double *bet, double *hp));
X
X#undef P_
X
X/* given the mjd, find the geocentric ecliptic longitude, lam, and latitude,
X * bet, and horizontal parallax, hp for the moon.
X * N.B. series for long and lat are good to about 10 and 3 arcseconds. however,
X *   math errors cause up to 100 and 30 arcseconds error, even if use double.
X *   why?? suspect highly sensitive nature of difference used to get m1..6.
X * N.B. still need to correct for nutation. then for topocentric location
X *   further correct for parallax and refraction.
X */
Xvoid
Xmoon (mjd, lam, bet, hp)
Xdouble mjd;
Xdouble *lam, *bet, *hp;
X{
X    double t, t2;
X    double ld;
X    double ms;
X    double md;
X    double de;
X    double f;
X    double n;
X    double a, sa, sn, b, sb, c, sc, e, e2, l, g, w1, w2;
X    double m1, m2, m3, m4, m5, m6;
X
X    t = mjd/36525.;
X    t2 = t*t;
X
X    m1 = mjd/27.32158213;
X    m1 = 360.0*(m1-(long)m1);
X    m2 = mjd/365.2596407;
X    m2 = 360.0*(m2-(long)m2);
X    m3 = mjd/27.55455094;
X    m3 = 360.0*(m3-(long)m3);
X    m4 = mjd/29.53058868;
X    m4 = 360.0*(m4-(long)m4);
X    m5 = mjd/27.21222039;
X    m5 = 360.0*(m5-(long)m5);
X    m6 = mjd/6798.363307;
X    m6 = 360.0*(m6-(long)m6);
X
X    ld = 270.434164+m1-(.001133-.0000019*t)*t2;
X    ms = 358.475833+m2-(.00015+.0000033*t)*t2;
X    md = 296.104608+m3+(.009192+.0000144*t)*t2;
X    de = 350.737486+m4-(.001436-.0000019*t)*t2;
X    f = 11.250889+m5-(.003211+.0000003*t)*t2;
X    n = 259.183275-m6+(.002078+.000022*t)*t2;
X
X    a = degrad(51.2+20.2*t);
X    sa = sin(a);
X    sn = sin(degrad(n));
X    b = 346.56+(132.87-.0091731*t)*t;
X    sb = .003964*sin(degrad(b));
X    c = degrad(n+275.05-2.3*t);
X    sc = sin(c);
X    ld = ld+.000233*sa+sb+.001964*sn;
X    ms = ms-.001778*sa;
X    md = md+.000817*sa+sb+.002541*sn;
X    f = f+sb-.024691*sn-.004328*sc;
X    de = de+.002011*sa+sb+.001964*sn;
X    e = 1-(.002495+7.52e-06*t)*t;
X    e2 = e*e;
X
X    ld = degrad(ld);
X    ms = degrad(ms);
X    n = degrad(n);
X    de = degrad(de);
X    f = degrad(f);
X    md = degrad(md);
X
X    l = 6.28875*sin(md)+1.27402*sin(2*de-md)+.658309*sin(2*de)+
X        .213616*sin(2*md)-e*.185596*sin(ms)-.114336*sin(2*f)+
X        .058793*sin(2*(de-md))+.057212*e*sin(2*de-ms-md)+
X        .05332*sin(2*de+md)+.045874*e*sin(2*de-ms)+.041024*e*sin(md-ms);
X    l = l-.034718*sin(de)-e*.030465*sin(ms+md)+.015326*sin(2*(de-f))-
X        .012528*sin(2*f+md)-.01098*sin(2*f-md)+.010674*sin(4*de-md)+
X        .010034*sin(3*md)+.008548*sin(4*de-2*md)-e*.00791*sin(ms-md+2*de)-
X        e*.006783*sin(2*de+ms);
X    l = l+.005162*sin(md-de)+e*.005*sin(ms+de)+.003862*sin(4*de)+
X        e*.004049*sin(md-ms+2*de)+.003996*sin(2*(md+de))+
X        .003665*sin(2*de-3*md)+e*.002695*sin(2*md-ms)+
X        .002602*sin(md-2*(f+de))+e*.002396*sin(2*(de-md)-ms)-
X        .002349*sin(md+de);
X    l = l+e2*.002249*sin(2*(de-ms))-e*.002125*sin(2*md+ms)-
X        e2*.002079*sin(2*ms)+e2*.002059*sin(2*(de-ms)-md)-
X        .001773*sin(md+2*(de-f))-.001595*sin(2*(f+de))+
X        e*.00122*sin(4*de-ms-md)-.00111*sin(2*(md+f))+.000892*sin(md-3*de);
X    l = l-e*.000811*sin(ms+md+2*de)+e*.000761*sin(4*de-ms-2*md)+
X         e2*.000704*sin(md-2*(ms+de))+e*.000693*sin(ms-2*(md-de))+
X         e*.000598*sin(2*(de-f)-ms)+.00055*sin(md+4*de)+.000538*sin(4*md)+
X         e*.000521*sin(4*de-ms)+.000486*sin(2*md-de);
X    l = l+e2*.000717*sin(md-2*ms);
X    *lam = ld+degrad(l);
X    range (lam, 2*PI);
X
X    g = 5.12819*sin(f)+.280606*sin(md+f)+.277693*sin(md-f)+
X        .173238*sin(2*de-f)+.055413*sin(2*de+f-md)+.046272*sin(2*de-f-md)+
X        .032573*sin(2*de+f)+.017198*sin(2*md+f)+.009267*sin(2*de+md-f)+
X        .008823*sin(2*md-f)+e*.008247*sin(2*de-ms-f);
X    g = g+.004323*sin(2*(de-md)-f)+.0042*sin(2*de+f+md)+
X        e*.003372*sin(f-ms-2*de)+e*.002472*sin(2*de+f-ms-md)+
X        e*.002222*sin(2*de+f-ms)+e*.002072*sin(2*de-f-ms-md)+
X        e*.001877*sin(f-ms+md)+.001828*sin(4*de-f-md)-e*.001803*sin(f+ms)-
X        .00175*sin(3*f);
X    g = g+e*.00157*sin(md-ms-f)-.001487*sin(f+de)-e*.001481*sin(f+ms+md)+
X         e*.001417*sin(f-ms-md)+e*.00135*sin(f-ms)+.00133*sin(f-de)+
X         .001106*sin(f+3*md)+.00102*sin(4*de-f)+.000833*sin(f+4*de-md)+
X         .000781*sin(md-3*f)+.00067*sin(f+4*de-2*md);
X    g = g+.000606*sin(2*de-3*f)+.000597*sin(2*(de+md)-f)+
X        e*.000492*sin(2*de+md-ms-f)+.00045*sin(2*(md-de)-f)+
X        .000439*sin(3*md-f)+.000423*sin(f+2*(de+md))+
X        .000422*sin(2*de-f-3*md)-e*.000367*sin(ms+f+2*de-md)-
X        e*.000353*sin(ms+f+2*de)+.000331*sin(f+4*de);
X    g = g+e*.000317*sin(2*de+f-ms+md)+e2*.000306*sin(2*(de-ms)-f)-
X        .000283*sin(md+3*f);
X    w1 = .0004664*cos(n);
X    w2 = .0000754*cos(c);
X    *bet = degrad(g)*(1-w1-w2);
X
X    *hp = .950724+.051818*cos(md)+.009531*cos(2*de-md)+.007843*cos(2*de)+
X          .002824*cos(2*md)+.000857*cos(2*de+md)+e*.000533*cos(2*de-ms)+
X          e*.000401*cos(2*de-md-ms)+e*.00032*cos(md-ms)-.000271*cos(de)-
X          e*.000264*cos(ms+md)-.000198*cos(2*f-md);
X    *hp = *hp+.000173*cos(3*md)+.000167*cos(4*de-md)-e*.000111*cos(ms)+
X         .000103*cos(4*de-2*md)-.000084*cos(2*md-2*de)-
X         e*.000083*cos(2*de+ms)+.000079*cos(2*de+2*md)+.000072*cos(4*de)+
X         e*.000064*cos(2*de-ms+md)-e*.000063*cos(2*de+ms-md)+
X         e*.000041*cos(ms+de);
X    *hp = *hp+e*.000035*cos(2*md-ms)-.000033*cos(3*md-2*de)-
X         .00003*cos(md+de)-.000029*cos(2*(f-de))-e*.000029*cos(2*md+ms)+
X         e2*.000026*cos(2*(de-ms))-.000023*cos(2*(f-de)+md)+
X         e*.000019*cos(4*de-ms-md);
X    *hp = degrad(*hp);
X}
END_OF_FILE
  if test 5365 -ne `wc -c <'moon.c'`; then
    echo shar: \"'moon.c'\" unpacked with wrong size!
  fi
  # end of 'moon.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'\" \(1918 characters\)
  sed "s/^X//" >'moonnf.c' <<'END_OF_FILE'
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X
X#if defined(__STDC__) || defined(__cplusplus)
X#define P_(s) s
X#else
X#define P_(s) ()
X#endif
X
Xextern void cal_mjd P_((int mn, double dy, int yr, double *Mjd));
Xextern void mjd_cal P_((double Mjd, int *mn, double *dy, int *yr));
X
Xvoid moonnf P_((double mjd, double *mjdn, double *mjdf));
Xstatic void m P_((double t, double k, double *mjd));
X
X#undef P_
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 */
Xvoid
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 void
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 1918 -ne `wc -c <'moonnf.c'`; then
    echo shar: \"'moonnf.c'\" unpacked with wrong size!
  fi
  # end of 'moonnf.c'
fi
if test -f 'msgmenu.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'msgmenu.c'\"
else
  echo shar: Extracting \"'msgmenu.c'\" \(7013 characters\)
  sed "s/^X//" >'msgmenu.c' <<'END_OF_FILE'
X/* this file contains the code to put up misc messages: xe_msg().
X * everything goes to stdout and into the message dialog.
X * the latter can be toggled on/off from the main menu.
X */
X
X#include <stdio.h>
X#include <ctype.h>
X#include <math.h>
X#if defined(__STDC__)
X#include <stdlib.h>
X#endif
X#include <X11/Xlib.h>
X#include <Xm/Xm.h>
X#include <Xm/Form.h>
X#include <Xm/PushB.h>
X#include <Xm/Text.h>
X#include <Xm/MessageB.h>
X
X#if defined(__STDC__) || defined(__cplusplus)
X#define P_(s) s
X#else
X#define P_(s) ()
X#endif
X
Xextern void prompt_map_cb P_((Widget w, XtPointer client, XtPointer call));
Xextern void set_something P_((Widget w, char *resource, char *value));
Xextern void set_xmstring P_((Widget w, char *resource, char *txt));
X
Xvoid msg_manage P_((void));
Xvoid xe_msg P_((char *msg, int app_modal));
Xvoid msg_cursor P_((Cursor c));
Xstatic void msg_on_top P_((void));
Xstatic void msg_create_w P_((void));
Xstatic void msg_erase_cb P_((Widget w, XtPointer client, XtPointer call));
Xstatic void msg_close_cb P_((Widget w, XtPointer client, XtPointer call));
Xstatic void msg_add P_((char *msg));
Xstatic void msg_scroll_down P_((void));
Xstatic void msg_stdout P_((char *msg));
Xstatic void xe_msg_modal P_((char *p));
X
X#undef P_
X
Xextern Widget toplevel_w;
X
Xstatic Widget msg_w;        /* main form dialog widget */
Xstatic Widget txt_w;        /* scrolled text widget */
Xstatic int txtl;        /* current length of text in txt_w */
X
X/* called to force the scrolling message window to be up.
X * bring it down if it already is.
X */
Xvoid
Xmsg_manage()
X{
X    if (!msg_w)
X        msg_create_w();
X
X    if (!XtIsManaged(msg_w))
X        XtManageChild (msg_w);
X    else
X        XtUnmanageChild (msg_w);
X}
X
X/* if app_modal != 0
X *   make an application modal info box,
X * else
X *   create dialog if not already made.
X *   add message to the scrolling text window.
X *   also write message to stdout.
X */
Xvoid
Xxe_msg (msg, app_modal)
Xchar *msg;
Xint app_modal;
X{
X    if (app_modal) {
X        xe_msg_modal (msg);
X        return;
X    }
X
X    if (!msg_w)
X        msg_create_w();
X    msg_add (msg);
X
X    /* if we are managed we bully right under the cursor. */
X    if (XtIsManaged(msg_w))
X        msg_on_top();
X
X    msg_stdout (msg);
X}
X
X/* called to put up or remove the watch cursor.  */
Xvoid
Xmsg_cursor (c)
XCursor c;
X{
X    Window win;
X
X    if (msg_w && (win = XtWindow(msg_w))) {
X        Display *dsp = XtDisplay(msg_w);
X        if (c)
X        XDefineCursor (dsp, win, c);
X        else
X        XUndefineCursor (dsp, win);
X    }
X}
X
X/* do everything we can to make sure the message dialog is visible.
X */
Xstatic void
Xmsg_on_top()
X{
X    Window win;
X
X    /* bring to the top right under the cursor
X     * N.B. we can't if this is the first time but don't need to in that
X     *   case anyway because that means nothing is on top of us yet.
X     */
X    win = XtWindow(XtParent(msg_w));
X    if (win) {
X        XRaiseWindow(XtDisplay(XtParent(msg_w)), win);
X        XtCallCallbacks (msg_w, XmNmapCallback, NULL);
X    }
X
X    XFlush (XtDisplay(msg_w));
X}
X
X/* create the message dialog */
Xstatic void
Xmsg_create_w()
X{
X    static struct {
X        char *name;
X        void (*cb)();
X    } cb[] = {
X        {"Erase", msg_erase_cb},
X        {"Close", msg_close_cb},
X    };
X    Widget w;
X    Arg args[20];
X    int i, n;
X
X    /* make the help shell form-dialog widget */
X
X    n = 0;
X    XtSetArg (args[n], XmNdefaultPosition, False); n++;
X    XtSetArg (args[n], XmNautoUnmanage, False); n++;
X    XtSetArg (args[n], XmNfractionBase, 9); n++;
X    msg_w = XmCreateFormDialog (toplevel_w, "Message", args, n);
X    XtAddCallback (msg_w, XmNmapCallback, prompt_map_cb, NULL);
X    set_something (XtParent(msg_w), XmNtitle, "xephem Message");
X
X    /* make the control buttons */
X
X    for (i = 0; i < XtNumber(cb); i++) {
X        n = 0;
X        XtSetArg (args[n], XmNbottomAttachment, XmATTACH_FORM); n++;
X        XtSetArg (args[n], XmNleftAttachment, XmATTACH_POSITION); n++;
X        XtSetArg (args[n], XmNleftPosition, 1 + 5*i); n++;
X        XtSetArg (args[n], XmNrightAttachment, XmATTACH_POSITION); n++;
X        XtSetArg (args[n], XmNrightPosition, 3 + 5*i); n++;
X        w = XmCreatePushButton (msg_w, cb[i].name, args, n);
X        XtAddCallback (w, XmNactivateCallback, cb[i].cb, NULL);
X        XtManageChild (w);
X    }
X
X    /* make the scrolled text area to hold the messages */
X
X    n = 0;
X    XtSetArg (args[n], XmNtopAttachment, XmATTACH_FORM); n++;
X    XtSetArg (args[n], XmNbottomAttachment, XmATTACH_WIDGET); n++;
X    XtSetArg (args[n], XmNbottomWidget, w); n++;
X    XtSetArg (args[n], XmNleftAttachment, XmATTACH_FORM); n++;
X    XtSetArg (args[n], XmNrightAttachment, XmATTACH_FORM); n++;
X    XtSetArg (args[n], XmNeditMode, XmMULTI_LINE_EDIT); n++;
X    XtSetArg (args[n], XmNeditable, False); n++;
X    XtSetArg (args[n], XmNcursorPositionVisible, False); n++;
X    txt_w = XmCreateScrolledText (msg_w, "ScrolledText", args, n);
X    XtManageChild (txt_w);
X}
X
X/* callback from the erase pushbutton */
X/* ARGSUSED */
Xstatic void
Xmsg_erase_cb (w, client, call)
XWidget w;
XXtPointer client;
XXtPointer call;
X{
X    XmTextReplace (txt_w, 0, txtl, "");
X    txtl = 0;
X    XFlush (XtDisplay(msg_w));
X}
X
X/* callback from the close pushbutton */
X/* ARGSUSED */
Xstatic void
Xmsg_close_cb (w, client, call)
XWidget w;
XXtPointer client;
XXtPointer call;
X{
X    XtUnmanageChild (msg_w);
X}
X
X/* add msg to the txt_w widget.
X * Always set the vertical scroll bar to the extreme bottom.
X */
Xstatic void
Xmsg_add (msg)
Xchar *msg;
X{
X    int l;
X
X    l = strlen(msg);
X    if (l == 0)
X        return;
X
X    XmTextReplace (txt_w, txtl, txtl, msg);
X    txtl += l;
X
X    if (msg[l-1] != '\n')
X        msg_add ("\n");
X    else
X        msg_scroll_down();
X}
X
X/* make sure the text is scrolled to the bottom */
Xstatic void
Xmsg_scroll_down()
X{
X    XmTextSetInsertionPosition (txt_w, txtl);
X}
X
X/* write msg to stdout.
X * add newline if one is not already at the end of msg[].
X */
Xstatic void
Xmsg_stdout (msg)
Xchar *msg;
X{
X    int l = strlen(msg);
X
X    printf ("%s", msg);
X    if (l > 0 && msg[l-1] != '\n')
X        (void) putchar ('\n');
X    (void) fflush (stdout);
X}
X
X/* print a message, p, in an app-modal dialog. */
Xstatic void
Xxe_msg_modal (p)
Xchar *p;
X{
X    static Widget apmsg_w;
X    Arg args[20];
X    int n;
X
X    if (!apmsg_w) {
X        XmString button_string;
X        XmString title_string;
X        Widget w;
X
X        button_string = XmStringCreate("Ok", XmSTRING_DEFAULT_CHARSET);
X        title_string = XmStringCreate ("xephem Modal Message",
X                        XmSTRING_DEFAULT_CHARSET);
X
X        /* Create MessageBox dialog. */
X        n = 0;
X        XtSetArg (args[n], XmNdefaultPosition, False);  n++;
X        XtSetArg (args[n], XmNtitle, "xephem Modal Message");  n++;
X        XtSetArg (args[n], XmNokLabelString, button_string);  n++;
X        XtSetArg (args[n], XmNdialogStyle, XmDIALOG_APPLICATION_MODAL); n++;
X        apmsg_w = XmCreateWarningDialog (toplevel_w, "message", args, n);
X        XtAddCallback (apmsg_w, XmNmapCallback, prompt_map_cb, NULL);
X
X        XmStringFree (title_string);
X        XmStringFree (button_string);
X
X        w = XmMessageBoxGetChild (apmsg_w, XmDIALOG_CANCEL_BUTTON);
X        XtUnmanageChild (w);
X        w = XmMessageBoxGetChild (apmsg_w, XmDIALOG_HELP_BUTTON);
X        XtUnmanageChild (w);
X    }
X
X    set_xmstring (apmsg_w, XmNmessageString, p);
X
X    /* Display help window. rely on autoUnmanage to bring back down. */
X    XtManageChild (apmsg_w);
X}
END_OF_FILE
  if test 7013 -ne `wc -c <'msgmenu.c'`; then
    echo shar: \"'msgmenu.c'\" unpacked with wrong size!
  fi
  # end of 'msgmenu.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'\" \(2135 characters\)
  sed "s/^X//" >'nutation.c' <<'END_OF_FILE'
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X#include "preferences.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 */
Xvoid
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 (pref_get(PREF_ALGO) == PREF_FAST) {
X        *deps = 0.0;
X        *dpsi = 0.0;
X        return;
X    }
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*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 2135 -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'\" \(426 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 */
Xvoid
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 426 -ne `wc -c <'obliq.c'`; then
    echo shar: \"'obliq.c'\" unpacked with wrong size!
  fi
  # end of 'obliq.c'
fi
if test -f 'patchlevel.h' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'patchlevel.h'\"
else
  echo shar: Extracting \"'patchlevel.h'\" \(266 characters\)
  sed "s/^X//" >'patchlevel.h' <<'END_OF_FILE'
X/* this string constant defines the current version of xephem.
X * it is used by versionmenu.c to display the version in the Intro banner.
X * it is also updated by any subsequent patches that are issued for xephem.
X */
X
X#define    PATCHLEVEL    "Version 2.4b May 10, 1993"
END_OF_FILE
  if test 266 -ne `wc -c <'patchlevel.h'`; then
    echo shar: \"'patchlevel.h'\" unpacked with wrong size!
  fi
  # end of 'patchlevel.h'
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'\" \(5005 characters\)
  sed "s/^X//" >'pelement.c' <<'END_OF_FILE'
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X
X#if defined(__STDC__) || defined(__cplusplus)
X#define P_(s) s
X#else
X#define P_(s) ()
X#endif
X
Xextern void range P_((double *v, double r));
X
Xvoid pelement P_((double mjd, double plan[8][9]));
X
X#undef P_
X
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 */
Xvoid
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 5005 -ne `wc -c <'pelement.c'`; then
    echo shar: \"'pelement.c'\" unpacked with wrong size!
  fi
  # end of 'pelement.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'\" \(1974 characters\)
  sed "s/^X//" >'reduce.c' <<'END_OF_FILE'
X#include <math.h>
X#include "astro.h"
X
X
X#if defined(__STDC__) || defined(__cplusplus)
X#define P_(s) s
X#else
X#define P_(s) ()
X#endif
X
Xextern void range P_((double *v, double r));
X
Xvoid reduce_elements P_((double mjd0, double mjd, double inc0, double ap0, double om0, double *inc, double *ap, double *om));
X
X#undef P_
X
X/* convert those orbital elements that change from epoch mjd0 to epoch mjd.
X */
Xvoid
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 1974 -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'\" \(1867 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 */
Xvoid
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 */
Xvoid
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 1867 -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'\" \(4887 characters\)
  sed "s/^X//" >'riset.c' <<'END_OF_FILE'
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X
X#if defined(__STDC__) || defined(__cplusplus)
X#define P_(s) s
X#else
X#define P_(s) ()
X#endif
X
Xextern void range P_((double *v, double r));
X
Xvoid riset P_((double ra, double dec, double lat, double dis, double *lstr, double *lsts, double *azr, double *azs, int *status));
X
X#undef P_
X
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 lower 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 */
Xvoid
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)) != 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 4887 -ne `wc -c <'riset.c'`; then
    echo shar: \"'riset.c'\" unpacked with wrong size!
  fi
  # end of 'riset.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'\" \(2228 characters\)
  sed "s/^X//" >'sun.c' <<'END_OF_FILE'
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X
X#if defined(__STDC__) || defined(__cplusplus)
X#define P_(s) s
X#else
X#define P_(s) ()
X#endif
X
Xextern void anomaly P_((double ma, double s, double *nu, double *ea));
Xextern void range P_((double *v, double r));
X
Xvoid sunpos P_((double mjd, double *lsn, double *rsn));
X
X#undef P_
X
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 */
Xvoid
Xsunpos (mjd, lsn, rsn)
Xdouble mjd;
Xdouble *lsn, *rsn;
X{
X    static double last_mjd = -3691, last_lsn, last_rsn;
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    if (mjd == last_mjd) {
X        *lsn = last_lsn;
X        *rsn = last_rsn;
X        return;
X    }
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    range (lsn, 2*PI);
X    last_lsn = *lsn;
X    *rsn = last_rsn = 1.0000002*(1-s*cos(ea))+dr;
X    last_mjd = mjd;
X}
END_OF_FILE
  if test 2228 -ne `wc -c <'sun.c'`; then
    echo shar: \"'sun.c'\" unpacked with wrong size!
  fi
  # end of 'sun.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'\" \(1600 characters\)
  sed "s/^X//" >'utc_gst.c' <<'END_OF_FILE'
X#include "astro.h"
X
X#if defined(__STDC__) || defined(__cplusplus)
X#define P_(s) s
X#else
X#define P_(s) ()
X#endif
X
Xextern void cal_mjd P_((int mn, double dy, int yr, double *Mjd));
Xextern void mjd_cal P_((double Mjd, int *mn, double *dy, int *yr));
Xextern void range P_((double *v, double r));
X
Xvoid utc_gst P_((double mjd, double utc, double *gst));
Xvoid gst_utc P_((double mjd, double gst, double *utc));
Xstatic double tnaught P_((double mjd));
X
X#undef P_
X
X
X
X/* given a modified julian date, mjd, and a universally coordinated time, utc,
X * return greenwich mean siderial time, *gst.
X */
Xvoid
Xutc_gst (mjd, utc, gst)
Xdouble mjd;
Xdouble utc;
Xdouble *gst;
X{
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 */
Xvoid
Xgst_utc (mjd, gst, utc)
Xdouble mjd;
Xdouble gst;
Xdouble *utc;
X{
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 1600 -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 21 \(of 21\).
cp /dev/null ark21isdone
MISSING=""
for I in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 ; do
    if test ! -f ark${I}isdone ; then
    MISSING="${MISSING} ${I}"
    fi
done
if test "${MISSING}" = "" ; then
    echo You have unpacked all 21 archives.
    rm -f ark[1-9]isdone ark[1-9][0-9]isdone
    echo Building ephem.db
    cat > ephem.db.Z.uu ephem.db.Z.uu.?
    uudecode ephem.db.Z.uu
    rm ephem.db.Z.uu ephem.db.Z.uu.?
    uncompress ephem.db.Z
    echo Building skyviewmenu.c
    cat > skyviewmenu.c skyviewmenu.c.?
    rm skyviewmenu.c.?
    echo Building smallfm.xbm
    cat > smallfm.xbm smallfm.xbm.?
    rm smallfm.xbm.?
else
    echo You still must unpack the following archives:
    echo "        " ${MISSING}
fi
exit 0
exit 0 # Just in case...
--
  // [e-mail]            | Send comp.sources.x submissions to:
\X/  Amiga - The only way to fly!      |
 "It's intuitively obvious to the most |    [e-mail]
  casual observer..."                  |