Source Code 1994 March/Source_Code_CD-ROM_Walnut_Creek_March_1994.iso/compsrcs/misc/volume39/remind/patch09d

Text File | 1993-10-04 | 24.3 KB | 686 lines

Newsgroups: comp.sources.misc
From: [e-mail] (David F. Skoll)
Subject: v39i118:  remind - A replacement for calendar, Patch09d/4
Message-ID: <[e-mail]>
X-Md4-Signature: bc2071672c4bbe2241298bea3dc0f5c3
Sender: [e-mail] (Kent Landfield)
Organization: Dept. of Electronics, Carleton University
Date: Mon, 4 Oct 1993 20:40:13 GMT
Approved: [e-mail]

Submitted-by: [e-mail] (David F. Skoll)
Posting-number: Volume 39, Issue 118
Archive-name: remind/patch09d
Environment: UNIX, MS-DOS, OS/2
Patch-To: remind: Volume 33, Issue 58-69

#! /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".
# Contents:  moon.c
# Wrapped by kent@sparky on Mon Oct  4 15:03:31 1993
PATH=/bin:/usr/bin:/usr/ucb:/usr/local/bin:/usr/lbin ; export PATH
echo If this archive is complete, you will see the following message:
echo '          "shar: End of archive 4 (of 4)."'
if test -f 'moon.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'moon.c'\"
else
  echo shar: Extracting \"'moon.c'\" \(22461 characters\)
  sed "s/^X//" >'moon.c' <<'END_OF_FILE'
X/***************************************************************/
X/*                                                             */
X/*  MOON.C                                                     */
X/*                                                             */
X/*  Calculations for figuring out moon phases.                 */
X/*                                                             */
X/*  This file is part of REMIND.                               */
X/*  Copyright (C) 1992, 1993 by David F. Skoll.                */
X/*                                                             */
X/***************************************************************/
X
X/* All of these routines were adapted from the program "moontool"
X   by John Walker, February 1988.  Here's the blurb from moontool:
X
X   ... The information is generally accurate to within ten
X   minutes.
X
X   The algorithms used in this program to calculate the positions Sun and
X   Moon as seen from the Earth are given in the book "Practical Astronomy
X   With Your Calculator" by Peter Duffett-Smith, Second Edition,
X   Cambridge University Press, 1981. Ignore the word "Calculator" in the
X   title; this is an essential reference if you're interested in
X   developing software which calculates planetary positions, orbits,
X   eclipses, and the like. If you're interested in pursuing such
X   programming, you should also obtain:
X
X   "Astronomical Formulae for Calculators" by Jean Meeus, Third Edition,
X   Willmann-Bell, 1985. A must-have.
X
X   "Planetary Programs and Tables from -4000 to +2800" by Pierre
X   Bretagnon and Jean-Louis Simon, Willmann-Bell, 1986. If you want the
X   utmost (outside of JPL) accuracy for the planets, it's here.
X
X   "Celestial BASIC" by Eric Burgess, Revised Edition, Sybex, 1985. Very
X   cookbook oriented, and many of the algorithms are hard to dig out of
X   the turgid BASIC code, but you'll probably want it anyway.
X
X   Many of these references can be obtained from Willmann-Bell, P.O. Box
X   35025, Richmond, VA 23235, USA. Phone: [phone]. In addition
X   to their own publications, they stock most of the standard references
X   for mathematical and positional astronomy.
X
X   This program was written by:
X
X      John Walker
X      Autodesk, Inc.
X      2320 Marinship Way
X      Sausalito, CA 94965
X      [phone] Ext. 829
X
X      Usenet: {sun!well}!acad!kelvin
X
X   This program is in the public domain: "Do what thou wilt shall be the
X   whole of the law". I'd appreciate receiving any bug fixes and/or
X   enhancements, which I'll incorporate in future versions of the
X   program. Please leave the original attribution information intact so
X   that credit and blame may be properly apportioned.
X
X*/
X#include "config.h"
X#ifdef HAVE_STDLIB_H
X#include <stdlib.h>
X#endif
X#include <stdio.h>
X#include <math.h>
X#include <time.h>
X#include "types.h"
X#include "protos.h"
X#include "expr.h"
X#include "globals.h"
X#include "err.h"
X
X/* Function prototypes */
XPRIVATE long jdate ARGS((int y, int mon, int day));
XPRIVATE double jtime ARGS((int y, int mon, int day, int hour, int min, int sec));
XPRIVATE void jyear ARGS((double td, int *yy, int *mm, int *dd));
XPRIVATE void jhms ARGS((double j, int *h, int *m, int *s));
XPRIVATE double meanphase ARGS((double sdate, double phase, double *usek));
XPRIVATE double truephase ARGS((double k, double phase));
XPRIVATE double kepler ARGS((double m, double ecc));
XPRIVATE double phase ARGS((double, double *, double *, double *, double *, double *, double *));
X
X
X/*  Astronomical constants  */
X
X#define epoch        2444238.5       /* 1980 January 0.0 */
X
X/*  Constants defining the Sun's apparent orbit  */
X
X#define elonge        278.833540       /* Ecliptic longitude of the Sun
X                      at epoch 1980.0 */
X#define elongp        282.596403       /* Ecliptic longitude of the Sun at
X                      perigee */
X#define eccent      0.016718       /* Eccentricity of Earth's orbit */
X#define sunsmax     1.495985e8     /* Semi-major axis of Earth's orbit, km */
X#define sunangsiz   0.533128       /* Sun's angular size, degrees, at
X                      semi-major axis distance */
X
X/*  Elements of the Moon's orbit, epoch 1980.0  */
X
X#define mmlong      64.975464      /* Moon's mean lonigitude at the epoch */
X#define mmlongp     349.383063       /* Mean longitude of the perigee at the
X                      epoch */
X#define mlnode        151.950429       /* Mean longitude of the node at the
X                      epoch */
X#define minc        5.145396       /* Inclination of the Moon's orbit */
X#define mecc        0.054900       /* Eccentricity of the Moon's orbit */
X#define mangsiz     0.5181         /* Moon's angular size at distance a
X                      from Earth */
X#define msmax       384401.0       /* Semi-major axis of Moon's orbit in km */
X#define mparallax   0.9507       /* Parallax at distance a from Earth */
X#define synmonth    29.53058868    /* Synodic month (new Moon to new Moon) */
X#define lunatbase   2423436.0      /* Base date for E. W. Brown's numbered
X                      series of lunations (1923 January 16) */
X
X/*  Properties of the Earth  */
X
X#define earthrad    6378.16       /* Radius of Earth in kilometres */
X
X#define PI 3.14159265358979323846
X
X/*  Handy mathematical functions  */
X
X#ifdef sgn
X#undef sgn
X#endif
X#define sgn(x) (((x) < 0) ? -1 : ((x) > 0 ? 1 : 0))      /* Extract sign */
X
X#ifdef abs
X#undef abs
X#endif
X#define abs(x) ((x) < 0 ? (-(x)) : (x))           /* Absolute val */
X
X#define fixangle(a) ((a) - 360.0 * (floor((a) / 360.0)))  /* Fix angle      */
X#define torad(d) ((d) * (PI / 180.0))              /* Deg->Rad      */
X#define todeg(d) ((d) * (180.0 / PI))              /* Rad->Deg      */
X#define dsin(x) (sin(torad((x))))              /* Sin from deg */
X#define dcos(x) (cos(torad((x))))              /* Cos from deg */
X
X/***************************************************************/
X/*                                                             */
X/*  jdate                                                      */
X/*                                                             */
X/*  Convert a date and time to Julian day and fraction.        */
X/*                                                             */
X/***************************************************************/
X#ifdef HAVE_PROTOS
XPRIVATE long jdate(int y, int mon, int day)
X#else
Xstatic long jdate(y, mon, day)
Xint y, mon, day;
X#endif
X{
X   long c, m;
X
X   m = mon+1;
X   if (m>2) {
X      m -= 3;
X   } else {
X      m += 9;
X      y--;
X   }
X   c = y/100L;   /* Century */
X   y -= 100L * c;
X   return day + (c*146097L)/4 + (y*1461L)/4 + (m*153L+2)/5 + 1721119L;
X}
X
X/***************************************************************/
X/*                                                             */
X/*  jtime                                                      */
X/*                                                             */
X/*  Convert a GMT date and time to astronomical Julian time,   */
X/*  i.e. Julian date plus day fraction, expressed as a double  */
X/*                                                             */
X/***************************************************************/
X#ifdef HAVE_PROTOS
XPRIVATE double jtime(int y, int mon, int day, int hour, int min, int sec)
X#else
Xstatic double jtime(y, mon, day, hour, min, sec)
Xint y, mon, day, hour, min, sec;
X#endif
X{
X   return (jdate(y, mon, day)-0.5) +
X             (sec + 60L * (long) min + 3600L * (long) hour) / 86400.0;
X}
X
X/***************************************************************/
X/*                                                             */
X/*  jyear                                                      */
X/*                                                             */
X/*  Convert a Julian date to year, month, day.                 */
X/*                                                             */
X/***************************************************************/
X#ifdef HAVE_PROTOS
XPRIVATE void jyear(double td, int *yy, int *mm, int *dd)
X#else
Xstatic void jyear(td, yy, mm, dd)
Xdouble td;
Xint *yy, *mm, *dd;
X#endif
X{
X   double j, d, y, m;
X
X   td += 0.5;         /* Astronomical to civil */
X   j = floor(td);
X   j = j - 1721119.0;
X   y = floor(((4 * j) - 1) / 146097.0);
X   j = (j * 4.0) - (1.0 + (146097.0 * y));
X   d = floor(j / 4.0);
X   j = floor(((4.0 * d) + 3.0) / 1461.0);
X   d = ((4.0 * d) + 3.0) - (1461.0 * j);
X   d = floor((d + 4.0) / 4.0);
X   m = floor(((5.0 * d) - 3) / 153.0);
X   d = (5.0 * d) - (3.0 + (153.0 * m));
X   d = floor((d + 5.0) / 5.0);
X   y = (100.0 * y) + j;
X   if (m < 10.0)
X      m = m + 2;
X   else {
X      m = m - 10;
X      y = y + 1;
X   }
X   *yy = y;
X   *mm = m;
X   *dd = d;
X}
X
X/***************************************************************/
X/*                                                             */
X/*  jhms                                                       */
X/*                                                             */
X/*  Convert a Julian time to hour, minutes and seconds.        */
X/*                                                             */
X/***************************************************************/
X#ifdef HAVE_PROTOS
XPRIVATE void jhms(double j, int *h, int *m, int *s)
X#else
Xstatic void jhms(j, h, m, s)
Xdouble j;
Xint *h, *m, *s;
X#endif
X{
X   long ij;
X
X   j += 0.5;         /* Astronomical to civil */
X   ij = (j - floor(j)) * 86400.0;
X   *h = ij / 3600L;
X   *m = (ij / 60L) % 60L;
X   *s = ij % 60L;
X}
X
X/***************************************************************/
X/*                                                             */
X/*  meanphase                                                  */
X/*                                                             */
X/*  Calculates mean phase of the Moon for a                    */
X/*  given base date and desired phase:                   */
X/*     0.0   New Moon                           */
X/*     0.25  First quarter                       */
X/*     0.5   Full moon                           */
X/*     0.75  Last quarter                       */
X/*  Beware!!!  This routine returns meaningless               */
X/*  results for any other phase arguments.  Don't           */
X/*  attempt to generalise it without understanding           */
X/*  that the motion of the moon is far more complicated           */
X/*  than this calculation reveals.                   */
X/*                                                             */
X/***************************************************************/
X#ifdef HAVE_PROTOS
XPRIVATE double meanphase(double sdate, double phase, double *usek)
X#else
Xstatic double meanphase(sdate, phase, usek)
Xdouble sdate, phase;
Xdouble *usek;
X#endif
X{
X   int yy, mm, dd;
X   double k, t, t2, t3, nt1;
X
X   jyear(sdate, &yy, &mm, &dd);
X
X   k = (yy + ((mm - 1) * (1.0 / 12.0)) - 1900) * 12.3685;
X
X   /* Time in Julian centuries from 1900 January 0.5 */
X   t = (sdate - 2415020.0) / 36525;
X   t2 = t * t;           /* Square for frequent use */
X   t3 = t2 * t;           /* Cube for frequent use */
X
X   *usek = k = floor(k) + phase;
X   nt1 = 2415020.75933 + synmonth * k
X         + 0.0001178 * t2
X         - 0.000000155 * t3
X         + 0.00033 * dsin(166.56 + 132.87 * t - 0.009173 * t2);
X
X   return nt1;
X}
X
X/***************************************************************/
X/*                                                             */
X/*  truephase                                                  */
X/*                                                             */
X/*  Given a K value used to determine the                      */
X/*  mean phase of the new moon, and a phase                    */
X/*  selector (0.0, 0.25, 0.5, 0.75), obtain                    */
X/*  the true, corrected phase time.                            */
X/*                                                             */
X/***************************************************************/
X#ifdef HAVE_PROTOS
XPRIVATE double truephase(double k, double phase)
X#else
Xstatic double truephase(k, phase)
Xdouble k, phase;
X#endif
X{
X   double t, t2, t3, pt, m, mprime, f;
X   int apcor = 0;
X
X   k += phase;           /* Add phase to new moon time */
X   t = k / 1236.85;       /* Time in Julian centuries from
X                     1900 January 0.5 */
X   t2 = t * t;           /* Square for frequent use */
X   t3 = t2 * t;           /* Cube for frequent use */
X   pt = 2415020.75933       /* Mean time of phase */
X        + synmonth * k
X        + 0.0001178 * t2
X        - 0.000000155 * t3
X        + 0.00033 * dsin(166.56 + 132.87 * t - 0.009173 * t2);
X
X        m = 359.2242               /* Sun's mean anomaly */
X       + 29.10535608 * k
X       - 0.0000333 * t2
X       - 0.00000347 * t3;
X        mprime = 306.0253          /* Moon's mean anomaly */
X       + 385.81691806 * k
X       + 0.0107306 * t2
X       + 0.00001236 * t3;
X        f = 21.2964                /* Moon's argument of latitude */
X       + 390.67050646 * k
X       - 0.0016528 * t2
X       - 0.00000239 * t3;
X   if ((phase < 0.01) || (abs(phase - 0.5) < 0.01)) {
X
X      /* Corrections for New and Full Moon */
X
X      pt +=     (0.1734 - 0.000393 * t) * dsin(m)
X           + 0.0021 * dsin(2 * m)
X           - 0.4068 * dsin(mprime)
X           + 0.0161 * dsin(2 * mprime)
X           - 0.0004 * dsin(3 * mprime)
X           + 0.0104 * dsin(2 * f)
X           - 0.0051 * dsin(m + mprime)
X           - 0.0074 * dsin(m - mprime)
X           + 0.0004 * dsin(2 * f + m)
X           - 0.0004 * dsin(2 * f - m)
X           - 0.0006 * dsin(2 * f + mprime)
X           + 0.0010 * dsin(2 * f - mprime)
X           + 0.0005 * dsin(m + 2 * mprime);
X      apcor = 1;
X   } else if ((abs(phase - 0.25) < 0.01 || (abs(phase - 0.75) < 0.01))) {
X      pt +=     (0.1721 - 0.0004 * t) * dsin(m)
X           + 0.0021 * dsin(2 * m)
X           - 0.6280 * dsin(mprime)
X           + 0.0089 * dsin(2 * mprime)
X           - 0.0004 * dsin(3 * mprime)
X           + 0.0079 * dsin(2 * f)
X           - 0.0119 * dsin(m + mprime)
X           - 0.0047 * dsin(m - mprime)
X           + 0.0003 * dsin(2 * f + m)
X           - 0.0004 * dsin(2 * f - m)
X           - 0.0006 * dsin(2 * f + mprime)
X           + 0.0021 * dsin(2 * f - mprime)
X           + 0.0003 * dsin(m + 2 * mprime)
X           + 0.0004 * dsin(m - 2 * mprime)
X           - 0.0003 * dsin(2 * m + mprime);
X      if (phase < 0.5)
X         /* First quarter correction */
X         pt += 0.0028 - 0.0004 * dcos(m) + 0.0003 * dcos(mprime);
X      else
X         /* Last quarter correction */
X         pt += -0.0028 + 0.0004 * dcos(m) - 0.0003 * dcos(mprime);
X      apcor = 1;
X   }
X   if (!apcor) return 0.0;
X   return pt;
X}
X
X/***************************************************************/
X/*                                                             */
X/*  kepler                                                     */
X/*                                                             */
X/*  Solve the equation of Kepler.                              */
X/*                                                             */
X/***************************************************************/
X#ifdef HAVE_PROTOS
XPRIVATE double kepler(double m, double ecc)
X#else
Xstatic double kepler(m, ecc)
Xdouble m, ecc;
X#endif
X{
X   double e, delta;
X#define EPSILON 1E-6
X
X   e = m = torad(m);
X   do {
X      delta = e - ecc * sin(e) - m;
X      e -= delta / (1 - ecc * cos(e));
X   } while (abs(delta) > EPSILON);
X   return e;
X}
X/***************************************************************/
X/*                                                             */
X/*  PHASE  --  Calculate phase of moon as a fraction:          */
X/*                                                             */
X/*   The argument is the time for which the phase is              */
X/*   Requested, expressed as a Julian date and               */
X/*   fraction.  Returns the terminator phase angle as a           */
X/*   percentage of a full circle (i.e., 0 to 1), and           */
X/*   stores into pointer arguments the illuminated           */
X/*   fraction of the Moon's disc, the Moon's age in           */
X/*   days and fraction, the distance of the Moon from           */
X/*   the centre of the Earth, and the angular diameter           */
X/*   subtended by the Moon as seen by an observer at           */
X/*   the centre of the Earth.                       */
X/*                                                             */
X/***************************************************************/
X#ifdef HAVE_PROTOS
XPRIVATE double phase(double pdate,
X                 double *pphase,
X             double *mage,
X             double *dist,
X             double *angdia,
X             double *sudist,
X             double *suangdia)
X#else
Xstatic double phase(pdate, pphase, mage, dist, angdia, sudist, suangdia)
Xdouble pdate;
Xdouble *pphase;            /* Illuminated fraction */
Xdouble *mage;               /* Age of moon in days */
Xdouble *dist;               /* Distance in kilometres */
Xdouble *angdia;            /* Angular diameter in degrees */
Xdouble *sudist;            /* Distance to Sun */
Xdouble *suangdia;                  /* Sun's angular diameter */
X#endif
X{
X
X   double Day, N, M, Ec, Lambdasun, ml, MM, MN, Ev, Ae, A3, MmP,
X          mEc, A4, lP, V, lPP, NP, y, x, Lambdamoon,
X          MoonAge, MoonPhase,
X          MoonDist, MoonDFrac, MoonAng,
X          F, SunDist, SunAng;
X
X        /* Calculation of the Sun's position */
X
X   Day = pdate - epoch;        /* Date within epoch */
X   N = fixangle((360 / 365.2422) * Day); /* Mean anomaly of the Sun */
X   M = fixangle(N + elonge - elongp);    /* Convert from perigee
X                      co-ordinates to epoch 1980.0 */
X   Ec = kepler(M, eccent);     /* Solve equation of Kepler */
X   Ec = sqrt((1 + eccent) / (1 - eccent)) * tan(Ec / 2);
X   Ec = 2 * todeg(atan(Ec));   /* 1 anomaly */
X   Lambdasun = fixangle(Ec + elongp);  /* Sun's geocentric ecliptic
X                          longitude */
X   /* Orbital distance factor */
X   F = ((1 + eccent * cos(torad(Ec))) / (1 - eccent * eccent));
X   SunDist = sunsmax / F;        /* Distance to Sun in km */
X        SunAng = F * sunangsiz;     /* Sun's angular size in degrees */
X
X
X        /* Calculation of the Moon's position */
X
X        /* Moon's mean longitude */
X   ml = fixangle(13.1763966 * Day + mmlong);
X
X        /* Moon's mean anomaly */
X   MM = fixangle(ml - 0.1114041 * Day - mmlongp);
X
X        /* Moon's ascending node mean longitude */
X   MN = fixangle(mlnode - 0.0529539 * Day);
X
X   /* Evection */
X   Ev = 1.2739 * sin(torad(2 * (ml - Lambdasun) - MM));
X
X   /* Annual equation */
X   Ae = 0.1858 * sin(torad(M));
X
X   /* Correction term */
X   A3 = 0.37 * sin(torad(M));
X
X   /* Corrected anomaly */
X   MmP = MM + Ev - Ae - A3;
X
X   /* Correction for the equation of the centre */
X   mEc = 6.2886 * sin(torad(MmP));
X
X   /* Another correction term */
X   A4 = 0.214 * sin(torad(2 * MmP));
X
X   /* Corrected longitude */
X   lP = ml + Ev + mEc - Ae + A4;
X
X   /* Variation */
X   V = 0.6583 * sin(torad(2 * (lP - Lambdasun)));
X
X   /* 1 longitude */
X   lPP = lP + V;
X
X   /* Corrected longitude of the node */
X   NP = MN - 0.16 * sin(torad(M));
X
X   /* Y inclination coordinate */
X   y = sin(torad(lPP - NP)) * cos(torad(minc));
X
X   /* X inclination coordinate */
X   x = cos(torad(lPP - NP));
X
X   /* Ecliptic longitude */
X   Lambdamoon = todeg(atan2(y, x));
X   Lambdamoon += NP;
X
X   /* Calculation of the phase of the Moon */
X
X   /* Age of the Moon in degrees */
X   MoonAge = lPP - Lambdasun;
X
X   /* Phase of the Moon */
X   MoonPhase = (1 - cos(torad(MoonAge))) / 2;
X
X   /* Calculate distance of moon from the centre of the Earth */
X
X   MoonDist = (msmax * (1 - mecc * mecc)) /
X      (1 + mecc * cos(torad(MmP + mEc)));
X
X        /* Calculate Moon's angular diameter */
X
X   MoonDFrac = MoonDist / msmax;
X   MoonAng = mangsiz / MoonDFrac;
X
X   if(pphase)   *pphase = MoonPhase;
X   if(mage)     *mage = synmonth * (fixangle(MoonAge) / 360.0);
X   if(dist)     *dist = MoonDist;
X   if(angdia)   *angdia = MoonAng;
X   if(sudist)   *sudist = SunDist;
X   if(suangdia) *suangdia = SunAng;
X   return fixangle(MoonAge) / 360.0;
X}
X
X/***************************************************************/
X/*                                                             */
X/*  MoonPhase                                                  */
X/*                                                             */
X/*  Interface routine dealing in Remind representations.       */
X/*  Given a local date and time, returns the moon phase at     */
X/*  that date and time as a number from 0 to 360.              */
X/*                                                             */
X/***************************************************************/
X#ifdef HAVE_PROTOS
XPUBLIC int MoonPhase(int date, int time)
X#else
Xint MoonPhase(date, time)
Xint date, time;
X#endif
X{
X   int utcd, utct;
X   int y, m, d;
X   double jd, mp;
X
X   /* Convert from local to UTC */
X   LocalToUTC(date, time, &utcd, &utct);
X
X   /* Convert from Remind representation to year/mon/day */
X   FromJulian(utcd, &y, &m, &d);
X
X   /* Convert to a true Julian date -- sorry for the name clashes! */
X   jd = jtime(y, m, d, (utct / 60), (utct % 60), 0);
X
X   /* Calculate moon phase */
X   mp = 360.0 * phase(jd, NULL, NULL, NULL, NULL, NULL, NULL);
X   return (int) mp;
X}
X
X/***************************************************************/
X/*                                                             */
X/*  HuntPhase                                                  */
X/*                                                             */
X/*  Given a starting date and time and a target phase, find    */
X/*  the first date on or after the starting date and time when */
X/*  the moon hits the specified phase.  Phase must be from     */
X/*  0 to 3 for new, 1stq, full, 3rdq                           */
X/*                                                             */
X/***************************************************************/
X#ifdef HAVE_PROTOS
XPUBLIC void HuntPhase(int startdate, int starttim, int phas, int *date, int *time)
X#else
Xvoid HuntPhase(startdate, starttim, phas, date, time)
Xint startdate, starttim, phas, *date, *time;
X#endif
X{
X   int utcd, utct;
X   int y, m, d;
X   int h, min, s;
X   int d1, t1;
X   double k1, k2, jd, jdorig;
X   double nt1, nt2;
X   /* Convert from local to UTC */
X   LocalToUTC(startdate, starttim, &utcd, &utct);
X
X   /* Convert from Remind representation to year/mon/day */
X   FromJulian(utcd, &y, &m, &d);
X
X   /* Convert to a true Julian date -- sorry for the name clashes! */
X   jdorig = jtime(y, m, d, (utct / 60), (utct % 60), 0);
X   jd = jdorig - 45;
X   nt1 = meanphase(jd, 0.0, &k1);
X   while(1) {
X      jd += synmonth;
X      nt2 = meanphase(jd, 0.0, &k2);
X      if (nt1 <= jdorig && nt2 > jdorig) break;
X      nt1 = nt2;
X      k1 = k2;
X   }
X   jd = truephase(k1, phas/4.0);
X   if (jd < jdorig) jd = truephase(k2, phas/4.0);
X
X   /* Convert back to Remind format */
X   jyear(jd, &y, &m, &d);
X   jhms(jd, &h, &min, &s);
X
X   d1 = Julian(y, m, d);
X   t1 = h*60 + min;
X   UTCToLocal(d1, t1, date, time);
X}
END_OF_FILE
  if test 22461 -ne `wc -c <'moon.c'`; then
    echo shar: \"'moon.c'\" unpacked with wrong size!
  fi
  # end of 'moon.c'
fi
echo shar: End of archive 4 \(of 4\).
cp /dev/null ark4isdone
MISSING=""
for I in 1 2 3 4 ; do
    if test ! -f ark${I}isdone ; then
    MISSING="${MISSING} ${I}"
    fi
done
if test "${MISSING}" = "" ; then
    echo You have unpacked all 4 archives.
    rm -f ark[1-9]isdone
else
    echo You still must unpack the following archives:
    echo "        " ${MISSING}
fi
exit 0
exit 0 # Just in case...