home *** CD-ROM | disk | FTP | other *** search
/ Frozen Fish 2: PC / frozenfish_august_1995.bin / bbs / d09xx / d0975.lha / PCal / moonphas.c < prev    next >
C/C++ Source or Header  |  1992-02-19  |  18KB  |  649 lines
  1. /*
  2.  *  moonphas.c - encapsulates routines used by Pcal for moon phase calculation
  3.  *
  4.  *  The following suite of routines are to calculate the phase of the moon
  5.  *  for a given month, day, year.  They compute the phase of the moon for
  6.  *  noon (UT) on the day requested, the start of the Julian day.
  7.  *
  8.  *  Contents:
  9.  *
  10.  *        calc_phase
  11.  *        find_moonfile
  12.  *        find_phase
  13.  *        read_moonfile
  14.  *
  15.  *  Revision history:
  16.  *
  17.  *    4.3a    SAG    12/27/91    Apply timezone offset to JD calculation
  18.  *
  19.  *    4.3    AWR    12/05/91    Search for moon file in Pcal's path
  20.  *                    as well as in calendar file's path
  21.  *
  22.  *            10/25/91    Many changes for support of moon
  23.  *                    phase wildcards and -Z flag
  24.  *
  25.  *    4.11    AWR    08/23/91    Revise is_quarter() to eliminate
  26.  *                    occasional missing or duplicate
  27.  *                    quarter-moons in "-m" mode; add
  28.  *                    gen_phases()
  29.  *
  30.  *    4.1    AWR    08/02/91    Fix bug in julday() where result is
  31.  *                    calculated incorrectly if floor() has
  32.  *                    no prototype
  33.  *
  34.  *      4.01    RLD     03/19/91        Upgraded calc_phase() to accurately
  35.  *                                      calculate the lunar phase for any
  36.  *                                      day without needing to resort to
  37.  *                                      a moon file.
  38.  *
  39.  *    4.0    AWR    03/07/91    Add find_moonfile()
  40.  *
  41.  *            01/15/91    Author: translated PostScript
  42.  *                    routines to C and added moon
  43.  *                    file routines
  44.  */
  45.  
  46. /*
  47.  *  Offset between local timezone and UT.  Units are hours west of Greenwich.
  48.  *    Example values are 5 for EST in the USA, 8 for PST in the USA and
  49.  *    -10 for EST in Australia.
  50.  *
  51.  */
  52. #define UTOFFSET 6        /* default is CST in the USA */
  53. /*
  54.  * Standard headers:
  55.  */
  56.  
  57. #include <stdio.h>
  58. #include <string.h>
  59. #include <ctype.h>
  60. #include <math.h>
  61.  
  62. /*
  63.  * Pcal-specific definitions:
  64.  */
  65.  
  66. #include "pcaldefs.h"
  67. #include "pcalglob.h"
  68. #include "pcallang.h"
  69.  
  70. /*
  71.  * Macros:
  72.  */
  73.  
  74. /*  Astronomical constants. */
  75.  
  76. #define epoch        2444238.5       /* 1980 January 0.0 */
  77.  
  78. /*  Constants defining the Sun's apparent orbit. */
  79.  
  80. #define elonge        278.833540       /* ecliptic longitude of the Sun
  81.                         at epoch 1980.0 */
  82. #define elongp        282.596403       /* ecliptic longitude of the Sun at
  83.                         perigee */
  84. #define eccent      0.016718       /* eccentricity of Earth's orbit */
  85.  
  86. /*  Elements of the Moon's orbit, epoch 1980.0. */
  87.  
  88. #define mmlong      64.975464      /* moon's mean lonigitude at the epoch */
  89. #define mmlongp     349.383063       /* mean longitude of the perigee at the
  90.                                         epoch */
  91. #define mlnode        151.950429       /* mean longitude of the node at the
  92.                         epoch */
  93. #define synmonth    29.53058868    /* synodic month (new Moon to new Moon) */
  94.  
  95. #define PI 3.14159265358979323846  /* assume not near black hole nor in
  96.                         Tennessee ;) */
  97.  
  98.  
  99. /*  Handy mathematical functions. */
  100.  
  101. #define sgn(x) (((x) < 0) ? -1 : ((x) > 0 ? 1 : 0))      /* extract sign */
  102. #ifndef abs
  103. #define abs(x) ((x) < 0 ? (-(x)) : (x))           /* absolute val */
  104. #endif
  105. #define fixangle(a) ((a) - 360.0 * (floor((a) / 360.0)))  /* fix angle      */
  106. #define torad(d) ((d) * (PI / 180.0))              /* deg->rad      */
  107. #define todeg(d) ((d) * (180.0 / PI))              /* rad->deg      */
  108. #define dsin(x) (sin(torad((x))))              /* sin from deg */
  109. #define dcos(x) (cos(torad((x))))              /* cos from deg */
  110. #define FNITG(x) (sgn (x) * floor (abs (x)))
  111.  
  112.  
  113. /* Misc. stuff, to disappear someday along with moon file routines */
  114.  
  115. #define HOUR        12        /* hour of day when phase calculated */
  116. #define PERIOD        synmonth    /* for backward compatibility */
  117.  
  118. /* interpolate phase for day "d" from moon_info array elements "n1" and "n2" */
  119. #define CALC_PHASE(d, n1, n2) \
  120.      moon_info[n1].phase + ((d) - moon_info[n1].doy) * \
  121.     ((moon_info[n2].phase - moon_info[n1].phase) / \
  122.      (moon_info[n2].doy - moon_info[n1].doy))
  123.  
  124. /* generate error message, close file, and quit */
  125. #define ERR_EXIT(msg)    \
  126.     if (1) { ERR(msg); fclose(fp); return FALSE; } else
  127.  
  128. /* day and phase sequence error conditions - cf. read_moonfile() */
  129. #define DAY_TOO_SOON    (nrec > 1 && doy < prevdoy + 6)
  130. #define DAY_TOO_LATE    (doy > prevdoy + 9)
  131. #define WRONG_PHASE    (nrec > 1 && ph != (prevph + 1) % 4)
  132.  
  133.  
  134. /*
  135.  * Globals:
  136.  */
  137.  
  138. typedef struct {
  139.     int    doy;     /* day of year (1..366) */
  140.     int    quarter; /* quarter (MOON_NM, MOON_1Q, etc.) */
  141.     double    phase;     /* moon phase (cycles since new moon prior to 1/1) */
  142. } MOON_INFO;
  143.  
  144. static MOON_INFO moon_info[60];        /* quarter moons for year + dummies */
  145.  
  146.  
  147. /*
  148.  *  Routines to accurately calculate the phase of the moon
  149.  *
  150.  *  Originally adapted from "moontool.c" by John Walker, Release 2.0.
  151.  *
  152.  *      This routine (calc_phase) and its support routines were adapted from
  153.  *      phase.c (v 1.2 88/08/26 22:29:42 jef) in the program "xphoon"        
  154.  *      (v 1.9 88/08/26 22:29:47 jef) by Jef Poskanzer and Craig Leres.
  155.  *      The necessary notice follows...
  156.  *
  157.  *      Copyright (C) 1988 by Jef Poskanzer and Craig Leres.
  158.  *
  159.  *      Permission to use, copy, modify, and distribute this software and its
  160.  *      documentation for any purpose and without fee is hereby granted,
  161.  *      provided that the above copyright notice appear in all copies and that
  162.  *      both that copyright notice and this permission notice appear in
  163.  *      supporting documentation.  This software is provided "as is" without
  164.  *      express or implied warranty.
  165.  *
  166.  *      These were added to "pcal" by RLD on 19-MAR-1991
  167.  */
  168.  
  169.  
  170. /*
  171.  *  julday -- calculate the julian date from input month, day, year
  172.  *      N.B. - The Julian date is computed for noon UT.
  173.  *
  174.  *      Adopted from Peter Duffett-Smith's book `Astronomy With Your
  175.  *      Personal Computer' by Rick Dyson 18-MAR-1991
  176.  */
  177. #ifdef PROTOS
  178. static double julday(int month,
  179.              int day,
  180.              int year)
  181. #else
  182. static double julday(month, day, year)
  183.     int month, day, year;
  184. #endif
  185. {
  186.         int mn1, yr1;
  187.         double a, b, c, d, djd;
  188.  
  189.         mn1 = month;
  190.         yr1 = year;
  191.         if ( yr1 < 0 ) yr1 = yr1 + 1;
  192.         if ( month < 3 )
  193.             {
  194.                 mn1 = month + 12;
  195.                 yr1 = yr1 - 1;
  196.             }
  197.         if (( year < 1582 ) ||
  198.             ( year == 1582  && month < 10 ) ||
  199.                 ( year == 1582  && month == 10 && day < 15.0 ))
  200.                     b = 0;
  201.                 else
  202.                 {
  203.                     a = floor (yr1 / 100.0);
  204.                     b = 2 - a + floor (a / 4);
  205.                 }
  206.         if ( yr1 >= 0 )
  207.             c = floor (365.25 * yr1) - 694025;
  208.         else
  209.             c = FNITG ((365.25 * yr1) - 0.75) - 694025;
  210.         d = floor (30.6001 * (mn1 + 1));
  211.         djd = b + c + d + day + 2415020.0;
  212.         return djd;
  213. }
  214.  
  215.  
  216. /*
  217.  *  kepler - solve the equation of Kepler
  218.  */
  219. #ifdef PROTOS
  220. static double kepler(double m,
  221.                      double ecc)
  222. #else
  223. static double kepler(m, ecc)
  224.     double m, ecc;
  225. #endif
  226. {
  227.     double e, delta;
  228. #define EPSILON 1E-6
  229.  
  230.     e = m = torad(m);
  231.     do {
  232.        delta = e - ecc * sin(e) - m;
  233.        e -= delta / (1 - ecc * cos(e));
  234.     } while (abs(delta) > EPSILON);
  235.     return e;
  236. }
  237.  
  238.  
  239. /*
  240.  *  calc_phase - calculate phase of moon as a fraction:
  241.  *
  242.  *    The argument is the time for which the phase is requested,
  243.  *    expressed as the month, day and year.  It returns the phase of
  244.  *    the moon (0.0 -> 0.99) with the ordering as New Moon, First Quarter,
  245.  *      Full Moon, and Last Quarter.
  246.  *
  247.  *    Converted from the subroutine phase.c used by "xphoon.c" (see
  248.  *      above disclaimer) into calc_phase() for use in "moonphas.c"
  249.  *    by Rick Dyson 18-MAR-1991
  250.  */
  251. #ifdef PROTOS
  252. double calc_phase(int month,
  253.           int inday,
  254.           int year)
  255. #else
  256. double calc_phase(month, inday, year)
  257.     int month, inday, year;
  258. #endif
  259. {
  260.     double Day, N, M, Ec, Lambdasun, ml, MM, MN, Ev, Ae, A3, MmP,
  261.            mEc, A4, lP, V, lPP, NP, y, x, MoonAge, pdate;
  262.  
  263.     /*  need to convert month, day, year into a Julian pdate */
  264.     pdate = julday (month, inday, year);
  265.     /*  Account for timezone (want JD for local noon not UT noon) */
  266.     pdate += UTOFFSET/24.0;
  267.  
  268.         /*  Calculation of the Sun's position. */
  269.  
  270.     Day = pdate - epoch;            /* date within epoch */
  271.     N = fixangle((360 / 365.2422) * Day);    /* mean anomaly of the Sun */
  272.     M = fixangle(N + elonge - elongp);      /* convert from perigee
  273.                              co-ordinates to epoch 1980.0 */
  274.     Ec = kepler(M, eccent);            /* solve equation of Kepler */
  275.     Ec = sqrt((1 + eccent) / (1 - eccent)) * tan(Ec / 2);
  276.     Ec = 2 * todeg(atan(Ec));        /* true anomaly */
  277.         Lambdasun = fixangle(Ec + elongp);    /* Sun's geocentric ecliptic
  278.                              longitude */
  279.  
  280.         /*  Calculation of the Moon's position. */
  281.  
  282.         /*  Moon's mean longitude. */
  283.     ml = fixangle(13.1763966 * Day + mmlong);
  284.  
  285.         /*  Moon's mean anomaly. */
  286.     MM = fixangle(ml - 0.1114041 * Day - mmlongp);
  287.  
  288.         /*  Moon's ascending node mean longitude. */
  289.     MN = fixangle(mlnode - 0.0529539 * Day);
  290.  
  291.     /*  Evection. */
  292.     Ev = 1.2739 * sin(torad(2 * (ml - Lambdasun) - MM));
  293.  
  294.     /*  Annual equation. */
  295.     Ae = 0.1858 * sin(torad(M));
  296.  
  297.     /*  Correction term. */
  298.     A3 = 0.37 * sin(torad(M));
  299.  
  300.     /*  Corrected anomaly. */
  301.     MmP = MM + Ev - Ae - A3;
  302.  
  303.     /*  Correction for the equation of the centre. */
  304.     mEc = 6.2886 * sin(torad(MmP));
  305.  
  306.     /*  Another correction term. */
  307.     A4 = 0.214 * sin(torad(2 * MmP));
  308.  
  309.     /*  Corrected longitude. */
  310.     lP = ml + Ev + mEc - Ae + A4;
  311.  
  312.     /*  Variation. */
  313.     V = 0.6583 * sin(torad(2 * (lP - Lambdasun)));
  314.  
  315.     /*  True longitude. */
  316.     lPP = lP + V;
  317.  
  318.     /*  Calculation of the phase of the Moon. */
  319.  
  320.     /*  Age of the Moon in degrees. */
  321.     MoonAge = lPP - Lambdasun;
  322.  
  323.         return (fixangle (MoonAge) / 360.0);
  324. }
  325.  
  326.  
  327. /*
  328.  *  is_quarter - if current phase of moon coincides with quarter moon, return
  329.  *  MOON_NM, MOON_1Q, etc.; otherwise return MOON_OTHER;
  330.  *  
  331.  */
  332. #ifdef PROTOS
  333. static int is_quarter(double prev,
  334.               double curr,
  335.               double next)
  336. #else
  337. static int is_quarter(prev, curr, next)
  338.     double prev, curr, next;
  339. #endif
  340. {
  341.     int quarter;
  342.     double phase, diff;
  343.  
  344.     if (curr < prev)    /* adjust phases for 1 -> 0 wraparound */
  345.         curr++;
  346.     if (next < prev)
  347.         next++;
  348.  
  349.     /* if a quarter moon has occurred between "prev" and "next", return
  350.      * TRUE if "curr" is closer to it than "prev" or "next" is    
  351.      */
  352.     for (quarter = 1; quarter <= 4; quarter++) 
  353.         if (prev < (phase = quarter/4.0) && next > phase &&
  354.             (diff = abs(curr - phase)) < phase - prev &&
  355.             diff < next - phase)
  356.             return quarter % 4;        /* MOON_NM == 0 */
  357.  
  358.     return MOON_OTHER;
  359. }
  360.  
  361.  
  362. /*
  363.  * Routines to read moon file and calculate moon phase from data within
  364.  */
  365.  
  366.  
  367. /*
  368.  * make_moonfile - create the base name for the moon file from the supplied
  369.  * template and year; fill it in and return pointer to it
  370.  */
  371. #ifdef PROTOS
  372. static char *make_moonfile(char *filename,
  373.                char *template,
  374.                int year)
  375. #else
  376. static char *make_moonfile(filename, template, year)
  377.     char *filename;        /* base file name (output) */
  378.     char *template;        /* file name template */
  379.     int year;        /* year to substitute for "%y" in template */
  380. #endif
  381. {
  382.     char *p;
  383.  
  384.     /* copy the template to the file name; replace "%y" (if present)
  385.      * with last two digits of year (as per other year substitutions)
  386.      */
  387.     strcpy(filename, template);
  388.     if ((p = strchr(filename, '%')) && p[1] == 'y') {
  389.         *p++ = '0' + (year / 10) % 10;
  390.         *p   = '0' + year % 10;
  391.     }
  392.  
  393.     return filename;
  394. }
  395.  
  396.  
  397. /*
  398.  * find_moonfile - look for moon file for specified year.  If it exists
  399.  * and is readable, return its full path name; else return NULL.  (There
  400.  * are admittedly ways to do this without attempting to open the file -
  401.  * cf. find_executable() in pcalutil.c - but they may not be portable.)
  402.  */
  403. #ifdef PROTOS
  404. char *find_moonfile(int year)
  405. #else
  406. char *find_moonfile(year)
  407.     int year;
  408. #endif
  409. {
  410.     static char filename[STRSIZ];
  411.     char path[STRSIZ], *pathlist[3], moonfile[STRSIZ];
  412.     FILE *fp;
  413.  
  414.     /* create list of paths for alt_fopen() to search */
  415.     pathlist[0] = mk_path(path, datefile);    /* datefile path */
  416.     pathlist[1] = progpath;            /* program path */
  417.     pathlist[2] = NULL;            /* terminate list */
  418.  
  419.     /* attempt to open the moon file */
  420.     fp = alt_fopen(filename, make_moonfile(moonfile, MOONFILE, year),
  421.                pathlist, "r");
  422.  
  423. #ifdef ALT_MOONFILE
  424.     if (!fp)             /* try again with alternate name */
  425.         fp = alt_fopen(filename, make_moonfile(moonfile, ALT_MOONFILE,
  426.                    year), pathlist, "r");
  427. #endif
  428.  
  429.     return fp ? (fclose(fp), filename) : NULL;
  430. }
  431.  
  432.  
  433. /*
  434.  * read_moonfile - looks for moon data file (in same directory as .calendar);
  435.  * if found, reads file, fills in moon_info[] and returns TRUE; if not found
  436.  * (or error encountered), returns FALSE
  437.  */
  438. #ifdef PROTOS
  439. int read_moonfile(int year)
  440. #else
  441. int read_moonfile(year)
  442.     int year;
  443. #endif
  444. {
  445.     char *filename;
  446.     int line, nrec, month, day, hh, mm;
  447.     int ph, prevph = MOON_OTHER, doy, prevdoy, n, quarter;
  448.     double phase;
  449.     FILE *fp;
  450.     static char *words[MAXWORD];    /* avoid conflicts with globals */
  451.     static char lbuf[LINSIZ];
  452.  
  453.     if (! *datefile)            /* skip if no datefile */
  454.         return FALSE;
  455.  
  456.     /* get name of moon file and attempt to open it */
  457.  
  458.     if ((filename = find_moonfile(year)) == NULL ||
  459.         (fp = fopen(filename, "r")) == NULL) {
  460.         if (DEBUG(DEBUG_MOON) || DEBUG(DEBUG_PATHS))
  461.             FPR(stderr, "No moon file for %d\n", year);
  462.         return FALSE;
  463.     }
  464.  
  465.     if (DEBUG(DEBUG_MOON) || DEBUG(DEBUG_PATHS))
  466.         FPR(stderr, "Using moon file %s\n", filename);
  467.  
  468.     /*
  469.      * Moon file entries are of the form <phase> <date> {<time>}; each
  470.      * is converted below to a moon_info[] record (note that only the
  471.      * initial phase of the moon is directly calculated from <phase>;
  472.      * it is subsequently used only for error checking).  Dummy entries
  473.      * in moon_info[] precede and follow the information read from the
  474.      * moon file; these are used for subsequent interpolation of dates
  475.      * before the first / after the last quarter of the year.
  476.      */
  477.  
  478.     prevdoy = 0;
  479.     line = 0;
  480.  
  481.     for (nrec = 1; getline(fp, lbuf, &line); nrec++) {
  482.  
  483.         /* ensure line is recognizable */
  484.         if ((n = loadwords(words, lbuf)) < 2 ||
  485.             (ph = get_phase(words[0])) == MOON_OTHER)
  486.             ERR_EXIT(E_INV_LINE);
  487.  
  488.         if (nrec == 1)            /* phase at initial quarter */
  489.             quarter = ph == MOON_NM ? 4 : ph;
  490.  
  491.         /* extract the month and day fields (in appropriate order) */
  492.  
  493.         (void) split_date(words[1],
  494.                   date_style == USA_DATES ? &month : &day,
  495.                   date_style == USA_DATES ? &day : &month,
  496.                   NULL);
  497.  
  498.         /* validate the date and phase */
  499.  
  500.         if (!is_valid(month, day, year))    /* date OK? */
  501.             ERR_EXIT(E_INV_DATE);
  502.  
  503.         doy = DAY_OF_YEAR(month, day, year);    /* in sequence? */
  504.         if (DAY_TOO_SOON || DAY_TOO_LATE || WRONG_PHASE)
  505.             ERR_EXIT(E_DATE_SEQ);
  506.  
  507.         prevdoy = doy;            /* save for sequence check */
  508.         prevph = ph;
  509.  
  510.         /* calculate moon phase, factoring in time (if present) */
  511.  
  512.         phase = 0.25 * quarter++;
  513.         if (n > 2) {            /* extract hour and minute */
  514.             (void) split_date(words[2], &hh, &mm, NULL);
  515.             phase += (HOUR - (hh + (mm / 60.0))) / (24 * PERIOD); 
  516.         }
  517.         moon_info[nrec].doy = doy;    /* enter day and phase */
  518.         moon_info[nrec].quarter = ph % 4;
  519.         moon_info[nrec].phase = phase;
  520.  
  521.         if (DEBUG(DEBUG_MOON))
  522.             FPR(stderr, "moon_info[%2d]: %3d %d %6.3f\n", nrec,
  523.                 moon_info[nrec].doy, moon_info[nrec].quarter,
  524.                 moon_info[nrec].phase);
  525.     }
  526.  
  527.     /* check to see that file is all there */
  528.  
  529.     doy = YEAR_LEN(year) + 1;    /* day after end of year */
  530.     if (DAY_TOO_LATE)
  531.         ERR_EXIT(E_PREM_EOF);
  532.  
  533.     /* extrapolate dummy entries from nearest lunar month */
  534.  
  535.     moon_info[nrec].doy = doy;    /* day after end of year */
  536.     moon_info[nrec].quarter = MOON_OTHER;
  537.     moon_info[nrec].phase = CALC_PHASE(doy, nrec-5, nrec-1);
  538.     if (DEBUG(DEBUG_MOON))
  539.         FPR(stderr, "moon_info[%2d]: %3d %d %6.3f\n", nrec,
  540.             moon_info[nrec].doy, moon_info[nrec].quarter,
  541.             moon_info[nrec].phase);
  542.  
  543.     moon_info[0].doy = 0;        /* day before start of year */
  544.     moon_info[0].quarter = MOON_OTHER;
  545.     moon_info[0].phase = CALC_PHASE(0, 1, 5);
  546.     if (DEBUG(DEBUG_MOON))
  547.         FPR(stderr, "moon_info[%2d]: %3d %d %6.3f\n", 0,
  548.             moon_info[0].doy, moon_info[0].quarter,
  549.             moon_info[0].phase);
  550.     
  551.     fclose(fp);
  552.     return TRUE;
  553. }
  554.  
  555.  
  556. /*
  557.  * gen_phases - fill array with moon phases for all days in month/year (plus
  558.  * extra entries at beginning and end for last day of previous month and
  559.  * first day of next month, respectively)
  560.  */
  561. #ifdef PROTOS
  562. static void gen_phases(double phase[],
  563.                int month,
  564.                int year)
  565. #else
  566. static void gen_phases(phase, month, year)
  567.     double phase[];
  568.     int month, year;
  569. #endif
  570. {
  571.     int day, len;
  572.     DATE date;
  573.  
  574.     /* start with moon phase for last day of previous month */
  575.     MAKE_DATE(date, month, 0, year);
  576.     normalize(&date);
  577.     phase[0] = calc_phase(date.mm, date.dd, date.yy);
  578.  
  579.     /* add the moon phases for all days in the current month */
  580.     for (day = 1, len = LENGTH_OF(month, year); day <= len; day++)
  581.         phase[day] = calc_phase(month, day, year);
  582.  
  583.     /* append the moon phase for the first day of next month */
  584.     MAKE_DATE(date, month, len + 1, year);
  585.     normalize(&date);
  586.     phase[len + 1] = calc_phase(date.mm, date.dd, date.yy);
  587. }
  588.  
  589.  
  590. /*
  591.  * find_phase - look up phase of moon in moon phase file (if possible);
  592.  * otherwise calculate it using calc_phase() above.  Sets *pquarter to
  593.  * MOON_NM, MOON_1Q, etc. if quarter moon, MOON_OTHER if not
  594.  */
  595. #ifdef PROTOS
  596. double find_phase(int month,
  597.           int day,
  598.           int year,
  599.           int *pquarter)
  600. #else
  601. double find_phase(month, day, year, pquarter)
  602.     int month, day, year;
  603.     int *pquarter;
  604. #endif
  605. {
  606.     static int sv_year = 0, sv_month = 0;
  607.     static int use_file;
  608.     static double mphase[33];    /* 31 days + 2 dummies */
  609.     int i, doy;
  610.     double phase;
  611.  
  612.     if (year != sv_year) {        /* look for file for new year */
  613.         use_file = read_moonfile(year);
  614.     }
  615.  
  616.     if (! use_file) {        /* no file - calculate phase */
  617.         /* new month?  fill mphase[] with moon phases */
  618.         if (month != sv_month || year != sv_year) {
  619.             gen_phases(mphase, month, year);
  620.             sv_month = month;
  621.             sv_year = year;
  622.         }
  623.  
  624.         phase = mphase[day];
  625.         *pquarter = is_quarter(mphase[day-1], phase, mphase[day+1]);
  626.         return phase;
  627.     }
  628.  
  629.     /* moon file found - use the data read by read_moonfile() */
  630.  
  631.     sv_year = year;                /* avoid re-reading same file */
  632.     doy = DAY_OF_YEAR(month, day, year);
  633.     for (i = 1; doy > moon_info[i].doy; i++)    /* find interval */
  634.         ;
  635.  
  636.     /* if day appears in table, return exact value; else interpolate */
  637.  
  638.     if (doy == moon_info[i].doy) {
  639.         *pquarter = moon_info[i].quarter;
  640.         phase = moon_info[i].phase;
  641.     } else {
  642.         *pquarter = MOON_OTHER;
  643.         phase = CALC_PHASE(doy, i-1, i);
  644.     }
  645.  
  646.     return phase - (int)phase;            /* 0.0 <= phase < 1.0 */
  647. }
  648.  
  649.