home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / misc / volume17 / calentool / part02 / datelib.c next >
Encoding:
C/C++ Source or Header  |  1991-04-04  |  43.8 KB  |  1,813 lines

  1. /*
  2.  * $Header: datelib.c,v 2.5 91/01/30 12:45:28 billr Exp $
  3.  *
  4.  * datelib.c - Calendar (date) computation library
  5.  *
  6.  * R. P. C. Rodgers, UCSF, November 1988
  7.  * (with thanks to Greg Couch, Conrad Huang, and Dave Yee for their helpful
  8.  *  suggestions)
  9.  *
  10.  *        Copyright 1988, 1989 R. P. C. Rodgers,  All Rights Reserved
  11.  * Permission is granted by the author for use of this code under the following
  12.  * conditions: 1) improvements and corrections to this code are shared with the
  13.  * author; 2) the code is made freely available to anyone requesting it; 3) the
  14.  * code is not used for any profit-making venture without the express
  15.  * permission of the author; and, 3) all copies of the code will preserve the
  16.  * above authorship and copyright information and this notice.
  17.  *
  18.  * (requires math library)
  19.  * 
  20.  * Source of formulae:
  21.  * 
  22.  * 1) Schocken, Wolfgang Alexander
  23.  *    The calculated confusion of calendars; puzzles in Christian, Jewish and
  24.  *       Moslem calendars.
  25.  *    1st Ed.
  26.  *    Vantage Press
  27.  *    New York
  28.  *    1976
  29.  * 
  30.  * 2) Meeus, Jean
  31.  *    Astronomical Formulae for Calculators
  32.  *    Monografieen over Astronomie en Astrofysica
  33.  *    Volkssterrenwacht Urania V.Z.W.
  34.  *    Mattheessensstraat 62, B 2540 Hove, Belgium
  35.  *    Vereniging voor Sterrenkunde V.Z.W.
  36.  *    Ringlaan 3, B 1180 Brussel, Belgium
  37.  *    Vol. 4
  38.  *    Derde Druk
  39.  *    October 1980
  40.  *    (3rd edition of 1985 available from Willmann-Bell, address below)
  41.  *    (formulae for Easter, Julian and Gregorian date formation, and
  42.  *     solstices/equinoxes)
  43.  * 
  44.  * 3) Assorted contributions to the Hewlett-Packard HP-41C Software Library
  45.  * 
  46.  * In his 1987 Sun program "moontool," John Walker (Autodesk, Sausalito, CA)
  47.  * mentions several other potentially useful references:
  48.  * 
  49.  * 4) Pierre Bretagnon and Jean-Louis Simon
  50.  *    Planetary Programs and Tables from -4000 to +2800
  51.  *    Willmann-Bell
  52.  *    1986
  53.  *    (for utmost accuracy in planetary computations)
  54.  * 
  55.  * 5) Eric Burgess
  56.  *    Celestial BASIC
  57.  *    Revised Edition
  58.  *    Sybex
  59.  *    1985
  60.  *    (cookbook oriented, many algorithms hard to dig out of turgid BASIC)
  61.  * 
  62.  * Many of these references can be obtained from Willmann-Bell, P.O. Box
  63.  * 35025, Richmond, VA 23235, USA.  Phone: (804) 320-7016.  In addition
  64.  * to their own publications, they stock most of the standard references
  65.  * for mathematical and positional astronomy.
  66.  * 
  67.  * NOTES:
  68.  * check ranges on days (0-365)
  69.  * islamic new year and jewish dates not thoroughly tested
  70.  * 
  71.  */
  72.  
  73. #include "ct.h"        /* for the NO_HOLIDAYS #define */
  74.  
  75. #include    <stdio.h>    /* for NULL */
  76. #include    <math.h>
  77.  
  78. double    julian_day();
  79.  
  80. #ifndef NO_HOLIDAYS
  81. extern    char    *malloc();
  82. double  easter_offset(), passover_offset();
  83.  
  84. static char    *monthname[] = { 
  85.         "January", "February", "March", "April", "May",
  86.         "June", "July", "August", "September",
  87.         "October", "November", "December"    };
  88. static char    *dayname[] = { 
  89.         "Sunday", "Monday", "Tuesday", "Wednesday",
  90.         "Thursday", "Friday", "Saturday"    };
  91. static char    timebuf[16];
  92. static double    passoverJD, easterJD;
  93. static int    passoverJY;
  94.  
  95. /*
  96.  * date_string:
  97.  * Return date string of form: "Monday 12 December 1988"
  98.  * given day of week (0-6), day (1-31), month (1-12), year
  99.  */
  100. char *
  101. date_string(day_of_week, day, month, year)
  102.     double    day;
  103.     int    day_of_week, month, year;
  104. {
  105.     char *date;
  106.  
  107.     date = (char *) malloc(81 * sizeof(char));
  108.     if (date == NULL)
  109.         err_rpt("out of memory", FATAL);
  110.     (void) sprintf(date, "%s %d %s %d",
  111.         dayname[day_of_week], (int) day, monthname[--month], year);
  112.     return date;
  113. }
  114.  
  115. /*
  116.  * date_string_2:
  117.  * Return date string of form: "Monday 12 December 1988 (NYEAR)"
  118.  * given day of week (0-6), day (1-31), month (1-12), year, and alternate
  119.  * (that is, non-Gregorian) year, NYEAR
  120.  */
  121. char *
  122. date_string_2(day_of_week, day, month, year, nyear)
  123.     double    day;
  124.     int    day_of_week, month, nyear, year;
  125. {
  126.     char *date;
  127.  
  128.     date = (char *) malloc(81 * sizeof(char));
  129.     if (date == NULL)
  130.         err_rpt("out of memory", FATAL);
  131.     (void) sprintf(date, "%s %d %s %d (%d)",
  132.         dayname[day_of_week], (int) day, monthname[--month],
  133.             year, nyear);
  134.     return date;
  135. }
  136.  
  137. /*
  138.  * date_time_string:
  139.  * Return date/time string of form: "10:42 Monday 12 December 1988"
  140.  * given hour (0-24), minute (0-59), day of week (0-6), day (1-31),
  141.  * month (1-12), year
  142.  */
  143. char *
  144. date_time_string(hour, min, day_of_week, day, month, year)
  145.     double    day;
  146.     int    day_of_week, hour, min, month, year;
  147. {
  148.     char *date;
  149.  
  150.     date = (char *) malloc(81 * sizeof(char));
  151.     if (date == NULL)
  152.         err_rpt("out of memory", FATAL);
  153.     (void) sprintf(date, "%02d:%02d %s %d %s %d",
  154.         hour, min, dayname[day_of_week], (int) day,
  155.         monthname[--month], year);
  156.     return date;
  157. }
  158.  
  159. /*
  160.  * election_day:
  161.  * Compute date of Election Day given year (1584-9999)
  162.  * Election day is defined as the first Tuesday following the first
  163.  * Monday in November.
  164.  */
  165. double
  166. election_day(year)
  167.     int    year;
  168. {
  169.     double day, nth_mday_of_month();
  170.  
  171.     /* find first Tuesday in Nov. */
  172.     day = nth_mday_of_month(1, 2, 11, year);
  173.     /* if it's the first day of the month then Election day is next week */
  174.     if (day == 1.)
  175.         day = 8.;
  176.     return julian_day(day, 11, year);
  177. }
  178.  
  179. /*
  180.  * dday2:
  181.  * Compute dday2 for various holiday routines given x value, year (>1583)
  182.  */
  183. dday2(x, year)
  184.     int    x, year;
  185. {
  186.     int    atmp, btmp;
  187.  
  188.     atmp = 1.25 * year;
  189.     btmp = year / 100;
  190.     btmp = 0.75 * (1 + btmp);
  191.     return (x + atmp - btmp) % 7;
  192. }
  193. #endif    /* NO_HOLIDAYS */
  194.  
  195. /*
  196.  * days_remaining_in_year:
  197.  * Compute remaining days of year (0-365)
  198.  * given day (1-31), month (1-12), year (1901-2009)
  199.  */
  200. days_remaining_in_year(day, month, year)
  201.     double    day;
  202.     int    month, year;
  203. {
  204.     int    length_of_year(), day_of_year();
  205.  
  206.     return length_of_year(year) - day_of_year(day, month, year);
  207. }
  208.  
  209. /*
  210.  * length_of_year:
  211.  * Compute days in year (365 or 366)
  212.  * given year (1901-2009)
  213.  */
  214. length_of_year(year)
  215.     int    year;
  216. {
  217.     int    ylength;
  218.  
  219.     if ((year % 400) == 0)
  220.         ylength = 366;
  221.     else if ((year % 100) == 0)
  222.         ylength = 365;
  223.     else if ((year % 4) == 0)
  224.         ylength = 366;
  225.     else
  226.         ylength = 365;
  227.     return ylength;
  228. }
  229.  
  230. /*
  231.  * get_day_of_week:
  232.  * Compute day of week (0-6 for Sunday-Saturday)
  233.  * given day (1-31), month (1-12), year (1901-2009)
  234.  */
  235. get_day_of_week(day, month, year)
  236.     double    day;
  237.     int    month, year;
  238. {
  239.     int    atmp;
  240.  
  241.     atmp = julian_day(day, month, year) + 1.5;
  242.     return atmp % 7;
  243. }
  244.  
  245. /*
  246.  * day_of_year:
  247.  * Compute day of year (1-366)
  248.  * given day (1-31), month (1-12), year (1901-2009)
  249.  */
  250. day_of_year(day, month, year)
  251.     double    day;
  252.     int    month, year;
  253. {
  254.     return (int) julian_day(day, month, year)
  255.         - (int) julian_day(0.0, 1, year);
  256. }
  257.  
  258. /*
  259.  * nth_mday_of_month:
  260.  * Compute nth m-day of the month (1-31)
  261.  * given n (1-5), day of week (0-6, 0 for Sunday), month (1-12),
  262.  * year (1583-9999)
  263.  */
  264. double
  265. nth_mday_of_month(n, day_of_week, month, year)
  266.     int    day_of_week, month, n, year;
  267. {
  268.     int    atmp, btmp, ctmp, dtmp, etmp, tmonth, tyear;
  269.  
  270.     if (month > 2) {
  271.         tmonth = month + 1;
  272.         tyear = year;
  273.     }
  274.     else {
  275.         tmonth = month + 13;
  276.         tyear = year - 1;
  277.     }
  278.     atmp = 2.6 * tmonth;
  279.     btmp = 1.25 * tyear;
  280.     ctmp = (tyear / 100) - 7;
  281.     dtmp = 0.75 * ctmp;
  282.     etmp = (day_of_week - atmp - btmp + dtmp) % 7;
  283.     return (double) (etmp + (n * 7));
  284. }
  285.  
  286. /*
  287.  * julian_day:
  288.  * Compute Julian day (>=1)
  289.  * given day (1-31), month (1-12), year (1901-2009)
  290.  *
  291.  * Notes: The Gregorian reform is taken into account (that is, the day
  292.  * following 4 October 1582 is 15 October 1582; dates on this latter day or
  293.  * later are said to be in the Gregorian calendar).  B.C. years are counted
  294.  * astronomically (the year prior to year 1 is year 0).  The Julian day
  295.  * begins at GMT noon (the day is expressed in floating point).
  296.  * Example: to obtain Julian day for 4 Jul 1979: julian_day(4.0, 7, 1979)
  297.  */
  298. double
  299. julian_day(day, month, year)
  300.     double    day;
  301.     int    month, year;
  302. {
  303.     int    atmp, monthp, yearp;
  304.     double    ctmp = 1720994.5 + day;
  305.  
  306.     if (month > 2) {
  307.         monthp = month + 1;
  308.         yearp = year;
  309.     }
  310.     else {
  311.         monthp = month + 13;
  312.         yearp = year - 1;
  313.     }
  314.     if ((year > 1582) || (year == 1582 && month >= 10)
  315.         || (year == 1582 && month == 10 && day >= 15)) {
  316.         atmp = year / 100;
  317.         ctmp += 2 - atmp + (int)(atmp / 4);
  318.     }
  319.     ctmp += (int)(365.25 * yearp) + (int)(30.6001 * monthp);
  320.     return ctmp;
  321. }
  322.  
  323. #ifndef NO_HOLIDAYS
  324. /*
  325.  * datelib_int:
  326.  * Calculate often used quantities (e.g. Easter, Passover) as an
  327.  * optimization.
  328.  */
  329. datelib_init(year)
  330.     int    year;
  331. {
  332.     void passover_init(), easter_init();
  333.  
  334.     easter_init(year);
  335.     passover_init(year);
  336. }
  337.  
  338. /*
  339.  * corrected_julian_day:
  340.  * Correct Julian day (>=1) for conversion from JULIAN CALENDAR
  341.  * to GREGORIAN CALENDAR.
  342.  */
  343. double
  344. corrected_julian_day(jday)
  345.     double    jday;
  346. {
  347.     double    day;
  348.     int    atmp, btmp, month, year;
  349.  
  350.     gregorian_date(&day, &month, &year, jday);
  351.     atmp = year / 100;
  352.     btmp = ((3 * atmp) - 5) / 4;
  353.     return julian_day(day, month, year) + btmp;
  354. }
  355.  
  356. /*
  357.  * gregorian_date:
  358.  * Return pointers to day (1-31), month (1-12), year (1901-2009),
  359.  * given Julian day (>=0)
  360.  *
  361.  * Notes: The Gregorian reform is taken into account (that is, the day
  362.  * following 4 October 1582 is 15 October 1582; dates on this latter day or
  363.  * later are said to be in the Gregorian calendar).  B.C. years are counted
  364.  * astronomically (the year prior to year 1 is year 0).  The Julian day
  365.  * begins at GMT noon.  The Julian day can be expressed in
  366.  * floating point below to get a day result with decimal places.  This method
  367.  * is valid only for positive Julian day numbers.
  368.  */
  369. gregorian_date(p_day, p_month, p_year, jday)
  370.     int    *p_month, *p_year;
  371.     double    *p_day, jday;
  372. {
  373.     int    atmp, btmp, ctmp, dtmp, etmp, gtmp, ztmp;
  374.     double    ftmp;
  375.  
  376.     ztmp = jday + 0.5;
  377.     ftmp = (jday + 0.5) - ztmp;
  378.     if (ztmp >= 2299161) {
  379.         gtmp = (ztmp - 1867216.25) / 36524.25;
  380.         ctmp = gtmp / 4;
  381.         atmp = ztmp + 1 + gtmp - ctmp;
  382.     }
  383.     else
  384.         atmp = ztmp;
  385.     btmp = atmp + 1524;
  386.     ctmp = (btmp - 122.1) / 365.25;
  387.     dtmp = 365.25 * ctmp;
  388.     etmp = ((btmp - dtmp) / 30.6001);
  389.     ztmp = 30.6001 * etmp;
  390.     *p_day = btmp - dtmp - ztmp + ftmp;
  391.     if (etmp > 13.5)
  392.         *p_month = etmp - 13;
  393.     else
  394.         *p_month = etmp - 1;
  395.     if (*p_month > 2.5)
  396.         *p_year = ctmp - 4716;
  397.     else
  398.         *p_year = ctmp - 4715;
  399. }
  400.  
  401. /*
  402.  * mdays_between_dates:
  403.  * Compute number of mdays between two dates (exclusive: don't count the
  404.  * days themselves) given day of week (0-6, O for Sunday),
  405.  * two sets of day (1-31), month (1-12), year (1583-9999)
  406.  */
  407. mdays_between_dates(day_of_week, day1, month1, year1, day2, month2, year2)
  408.     int    day_of_week, month1, year1, month2, year2;
  409.     double    day1, day2;
  410. {
  411.     return wday(day_of_week, day2, month2, year2)
  412.         - wday(day_of_week, day1, month1, year1);
  413. }
  414.  
  415. /*
  416.  * years_date_is_mday:
  417.  * Compute year(s) for which a given date is an m-day
  418.  * given starting year, ending year,
  419.  * day of week (0-6, 0 for Sunday), day, and month
  420.  * algorithm said to be valid for 1 Mar 1900 to 28 Feb 2100.
  421.  */
  422. years_date_is_mday(day_of_week, day, month, start_year, end_year, yearlist,
  423.     number_of_years)
  424.     int    day_of_week, end_year, month, *number_of_years, yearlist[],
  425.         start_year;
  426.     double    day;
  427. {
  428.     int        diff, index, year, tdow;
  429.     static int    augment[] = {6, 11, 6, 5}; 
  430.  
  431.     *number_of_years = 0;
  432.     if (start_year > end_year) return -1;
  433.     for (year = start_year; year <= end_year; year++ ) {
  434.         tdow = get_day_of_week(day, month, year);
  435.         if (tdow == day_of_week) {
  436.             yearlist[(*number_of_years)++] = year;
  437.         }
  438.         if (*number_of_years == 2 || *number_of_years == 3) {
  439.             diff = yearlist[*number_of_years]
  440.                 - yearlist[*number_of_years - 1];
  441.             if (diff == 5) {
  442.                 year++;
  443.                 index = 0;
  444.                 break;
  445.             }
  446.             else if (diff == 11) {
  447.                 year++;
  448.                 index = 2;
  449.                 break;
  450.             }
  451.         }
  452.     }
  453.     for ( ; year <= end_year; year++ ) {
  454.         yearlist[(*number_of_years + 1)] =
  455.             yearlist[*number_of_years] + augment[index++ % 4];
  456.         (*number_of_years)++;
  457.     }
  458.     return 0;
  459. }
  460.  
  461. /*
  462.  * weekdays_between_dates:
  463.  * Compute weekdays between any two dates
  464.  * given two sets of day (1-31), month (1-12), year (1901-2009)
  465.  */
  466. weekdays_between_dates(day1, month1, year1, day2, month2, year2)
  467.     int    month1, month2, year1, year2;
  468.     double    day1, day2;
  469. {
  470.     return wday2(day2, month2, year2) - wday2(day1, month1, year1);
  471. }
  472.  
  473. /*
  474.  * wday:
  475.  * Compute wday for mdays_between_dates routine
  476.  * given day of week, day, month, year
  477.  */
  478. wday(day_of_week, day, month, year)
  479.     int    day_of_week, month, year;
  480.     double    day;
  481. {
  482.     int    atmp, btmp, ctmp, dtmp;
  483.  
  484.     atmp = dday(day, month, year) - day_of_week;
  485.     btmp = atmp / 7;
  486.     ctmp = atmp % 7;
  487.     dtmp = 0.11 * ctmp + 0.9;
  488.     return btmp + (0.5 * dtmp);
  489. }
  490.  
  491. /*
  492.  * wday2:
  493.  * Compute wday2 for weekdays_between_dates routine
  494.  * given day, month, year
  495.  */
  496. wday2(day, month, year)
  497.     int    month, year;
  498.     double    day;
  499. {
  500.     int    atmp, btmp, ctmp, dtmp;
  501.  
  502.     atmp = dday(day, month, year);
  503.     btmp = atmp / 7;
  504.     ctmp = atmp % 7;
  505.     dtmp = 1.801 * ctmp;
  506.     return (5 * btmp) + (0.5 * dtmp);
  507. }
  508.  
  509. /*
  510.  * dday:
  511.  * Compute dday for other routines
  512.  * given day (1-31), month (1-12), year (1901-2009)
  513.  */
  514. dday(day, month, year)
  515.     int    month, year;
  516.     double    day;
  517. {
  518.     int    atmp, btmp, ctmp, dtmp, tmonth, tyear;
  519.  
  520.     if (month > 2) {
  521.         tmonth = month + 1;
  522.         tyear = year;
  523.     }
  524.     else {
  525.         tmonth = month + 13;
  526.         tyear = year - 1;
  527.     }
  528.     atmp = tyear / 100;
  529.     btmp = 0.75 * (atmp - 7);
  530.     ctmp = 365.25 * tyear;
  531.     dtmp = 30.6001 * tmonth;
  532.     return (int) day - btmp + ctmp + dtmp;
  533. }
  534.  
  535. /*
  536.  * vernal_equinox:
  537.  * Compute Julian day for Vernal (March) Equinox given year (1901-2009)
  538.  */
  539. double
  540. vernal_equinox(year)
  541.     int    year;
  542. {
  543.     double    solstice_equinox_exact();
  544.  
  545.     return solstice_equinox_exact(0, year);
  546. }
  547.  
  548. /*
  549.  * summer_solstice:
  550.  * Compute Julian day for Summer Solstice (June) given year (1901-2009)
  551.  */
  552. double
  553. summer_solstice(year)
  554.     int    year;
  555. {
  556.     double    solstice_equinox_exact();
  557.  
  558.     return solstice_equinox_exact(1, year);
  559. }
  560.  
  561. /*
  562.  * autumn_equinox:
  563.  * Compute Julian day for Autumnal (September) Equinox given year (1901-2009)
  564.  */
  565. double
  566. autumn_equinox(year)
  567.     int    year;
  568. {
  569.     double    solstice_equinox_exact();
  570.  
  571.     return solstice_equinox_exact(2, year);
  572. }
  573.  
  574. /*
  575.  * winter_solstice:
  576.  * Compute Julian day for Winter (December) Solstice given year (1901-2009)
  577.  */
  578. double
  579. winter_solstice(year)
  580.     int    year;
  581. {
  582.     double    day, jday;
  583.     int    month, nyear;
  584.     double    solstice_equinox_exact();
  585.  
  586.     jday = solstice_equinox_exact(3, year);
  587.     gregorian_date(&day, &month, &nyear, jday);
  588.     if (nyear == year)
  589.         return jday;
  590.     else
  591.         return solstice_equinox_exact(3, (year + 1));
  592. }
  593.  
  594. /*
  595.  * solstice_equinox__exact:
  596.  * Compute more precise date for Solstice/Equinox
  597.  * given code (0-3, 0 for Vernal Equinox), year (1901-2009)
  598.  */
  599. double
  600. solstice_equinox_exact(code, year)
  601.     int    code, year;
  602. {
  603.     int    count;
  604.  
  605.     double    corr, jday_app;
  606.     double    solstice_equinox_approx(), solstice_equinox_correction();
  607.  
  608.     /* compute first approximation to Julian day */
  609.     jday_app = solstice_equinox_approx(code-1, year);
  610.     /* iteratively recompute corrected Julian day */
  611.     for(count = 0; count < 15; count++) {
  612.         corr = solstice_equinox_correction(jday_app, code);
  613.         jday_app += corr;
  614.         if (fabs(corr) < 0.00001)
  615.             break;
  616.     }
  617.     return jday_app;
  618. }
  619.  
  620. /*
  621.  * solstice_equinox_correction:
  622.  * Compute correction for Solstice/Equinox Julian day
  623.  * approximate Julian day, code (0-3, 0 for Vernal Equinox)
  624.  */
  625. double
  626. solstice_equinox_correction(jday, code)
  627.     int    code;
  628.     double    jday;
  629. {
  630.     double    apparent_longitude(), sin_degrees();
  631.  
  632.     return(58.0 * sin_degrees((code * 90.0) - apparent_longitude(jday)));
  633. }
  634.  
  635. /*
  636.  * apparent_longitude:
  637.  * Compute apparent longitude (true equinox) of Sun for
  638.  * Solstice/Equinox given approximate Julian day
  639.  */
  640. double
  641. apparent_longitude(jday)
  642.     double    jday;
  643. {
  644.     double    btmp, ctmp, dtmp, etmp, ftmp;
  645.     double    sin_degrees();
  646.  
  647.     btmp = (jday - 2415020.0) / 36525.0;    /* time in Julian centuries: */
  648.                         /* epoch 0.5 January 1900 */
  649.     ctmp = 279.69668 + (36000.76892 * btmp) /* geometric mean longitude */
  650.         + (0.0003025 * btmp * btmp);    /* (mean equinox of the date) */
  651.     dtmp = 358.47583 + (35999.04975 * btmp)    /* mean anomaly of Sun */
  652.         + (0.00015 * btmp * btmp)
  653.         + (0.0000033 * btmp * btmp * btmp);
  654.     etmp = (1.919460 - (0.004789 * btmp)    /* Sun's eqn of the center */
  655.         - (0.000014 * btmp)) * sin_degrees(dtmp)
  656.         + (0.020094 - (0.0001 * btmp)) * sin_degrees(2 * dtmp)
  657.         + 0.000293 * sin_degrees(3 * dtmp);
  658.     ftmp = ctmp + etmp;            /* Sun's true longitude */
  659.     return ftmp - 0.00569            /* app. long. of Sun */
  660.         - 0.00479 * sin_degrees(259.18    /* (true equinox of the date) */
  661.         - (1934.142 * btmp));
  662. }
  663.  
  664. /*
  665.  * sin_degrees:
  666.  * Compute sin for argument in degrees
  667.  */
  668. double
  669. sin_degrees(degrees)
  670. #define    PI_CORR    0.01745329252                /* pi / 180 */
  671.     double    degrees;
  672. {
  673.     return sin(degrees * PI_CORR);
  674. }
  675.  
  676. /*
  677.  * solstice_equinox_approx:
  678.  * Compute approximate date for Solstice/Equinox
  679.  * given code (0-3, 0 for Vernal Equinox), year (1901-2009)
  680.  */
  681. double
  682. solstice_equinox_approx(code, year)
  683.     int    code, year;
  684. {
  685.     return (365.2422 * (year + (code / 4))) + 1721141.3;
  686. }
  687.  
  688. /*
  689.  * easter:
  690.  * Return Julian date for Easter, given year (1901-2009)
  691.  *
  692.  * Offered in the book of Meeus, which cites:
  693.  *
  694.  * 1) Spencer Jones
  695.  *    General Astronomy
  696.  *    1922
  697.  *    pg. 73-74
  698.  *
  699.  * 2) Journal of the British Astronomical Association
  700.  *    Vol. 8
  701.  *    Pg. 91
  702.  *    Dec. 1977
  703.  *    (which attributes method to Butcher's Ecclesiastical Calendar of 1876)
  704.  *
  705.  * Method valid for all dates in the Gregorian calendar
  706.  * (from 15 October 1583 on)
  707.  */
  708. void
  709. easter_init(year)
  710.     int    year;
  711. {
  712.     double    day;
  713.     int    atmp, btmp, ctmp, dtmp, etmp, ftmp,
  714.         gtmp, htmp, itmp, ktmp, ltmp, mtmp;
  715.     int    month;
  716.  
  717.     atmp = year % 19;
  718.     btmp = year / 100;
  719.     ctmp = year % 100;
  720.     dtmp = btmp / 4;
  721.     etmp = btmp % 4;
  722.     ftmp = (btmp + 8) / 25;
  723.     gtmp = (btmp - ftmp + 1) / 3;
  724.     htmp = ((19 * atmp) + btmp - dtmp - gtmp + 15) % 30;
  725.     itmp = ctmp / 4;
  726.     ktmp = ctmp % 4;
  727.     ltmp = (32 + (2 * etmp) + (2 * itmp) - htmp - ktmp) % 7;
  728.     mtmp = (atmp + (11 * htmp) + (22 * ltmp)) / 451;
  729.     month = (htmp + ltmp - (7 * mtmp) + 114) / 31;
  730.     day = ((htmp + ltmp - (7 * mtmp) + 114) % 31) + 1;
  731.     easterJD = julian_day(day, month, year);
  732. }
  733.  
  734. double
  735. easter(year)
  736.     int    year;
  737. {
  738.     return easterJD;
  739. }
  740.  
  741. /*
  742.  * first_sunday_advent:
  743.  * Christian holidays: compute Julian day for First Sunday in Advent
  744.  * (closest Sunday to St. Andrew's day, the last day in November)
  745.  * given year (>1583)
  746.  */
  747. double
  748. first_sunday_advent(year)
  749.     int    year;
  750. {
  751.     int    atmp;
  752.  
  753.     atmp = get_day_of_week(30.0, 11, year);
  754.     if (atmp <= 3) {
  755.         return julian_day((30.0 - (double) atmp), 11, year);
  756.     }
  757.     else {
  758.         return julian_day(30.0, 11, year) + (7 - atmp);
  759.     }
  760. }
  761.  
  762. /*
  763.  * easter_offset:
  764.  * Christian holidays: compute Julian day as offset from Easter
  765.  * given year (>1583)
  766.  */
  767. double
  768. easter_offset(offset, year)
  769.     double    offset;
  770.     int    year;
  771. {
  772.     return easterJD + offset;
  773. }
  774.  
  775. /*
  776.  * septuagesima:
  777.  * Christian holidays: compute Julian day for Septuagesima Sunday
  778.  * (Third Sunday before Lent)
  779.  * given year (>1583)
  780.  */
  781. double
  782. septuagesima(year)
  783.     int    year;
  784. {
  785.     return easter_offset(-63.0, year);
  786. }
  787.  
  788. /*
  789.  * sexagesima:
  790.  * Christian holidays: compute Julian day for Sexagesima Sunday
  791.  * (Second Sunday before Lent)
  792.  * given year (>1583)
  793.  */
  794. double
  795. sexagesima(year)
  796.     int    year;
  797. {
  798.     return easter_offset(-56.0, year);
  799. }
  800.  
  801. /*
  802.  * quinquagesima:
  803.  * Christian holidays: compute Julian day for Quinquagesima (Shrove) Sunday
  804.  * (Sunday before Lent begins on Ash Wednesday)
  805.  * given year (>1583)
  806.  */
  807. double
  808. quinquagesima(year)
  809.     int    year;
  810. {
  811.     return easter_offset(-49.0, year);
  812. }
  813.  
  814. /*
  815.  * shrove_monday:
  816.  * Christian holidays: compute Julian day for Shrove Monday
  817.  * (Two days before Lent begins on Ash Wednesday)
  818.  * given year (>1583)
  819.  */
  820. double
  821. shrove_monday(year)
  822.     int    year;
  823. {
  824.     return easter_offset(-48.0, year);
  825. }
  826.  
  827. /*
  828.  * shrove_tuesday:
  829.  * Christian holidays: compute Julian day for Shrove Tuesday
  830.  * (Day before Lent begins on Ash Wednesday; Mardi Gras)
  831.  * given year (>1583)
  832.  */
  833. double
  834. shrove_tuesday(year)
  835.     int    year;
  836. {
  837.     return easter_offset(-47.0, year);
  838. }
  839.  
  840. /*
  841.  * ash_wednesday:
  842.  * Christian holidays: compute Julian day for Ash Wednesday
  843.  * given year (>1583)
  844.  */
  845. double
  846. ash_wednesday(year)
  847.     int    year;
  848. {
  849.     return easter_offset(-46.0, year);
  850. }
  851.  
  852. /*
  853.  * first_sunday_lent:
  854.  * Christian holidays: compute Julian day for First Sunday in Lent
  855.  * given year (>1583)
  856.  */
  857. double
  858. first_sunday_lent(year)
  859.     int    year;
  860. {
  861.     return easter_offset(-42.0, year);
  862. }
  863.  
  864. /*
  865.  * second_sunday_lent:
  866.  * Christian holidays: compute Julian day for Second Sunday in Lent
  867.  * given year (>1583)
  868.  */
  869. double
  870. second_sunday_lent(year)
  871.     int    year;
  872. {
  873.     return easter_offset(-35.0, year);
  874. }
  875.  
  876. /*
  877.  * third_sunday_lent:
  878.  * Christian holidays: compute Julian day for Third Sunday in Lent
  879.  * given year (>1583)
  880.  */
  881. double
  882. third_sunday_lent(year)
  883.     int    year;
  884. {
  885.     return easter_offset(-28.0, year);
  886. }
  887.  
  888. /*
  889.  * fourth_sunday_lent:
  890.  * Christian holidays: compute Julian day for Fourth Sunday in Lent
  891.  * given year (>1583)
  892.  */
  893. double
  894. fourth_sunday_lent(year)
  895.     int    year;
  896. {
  897.     return easter_offset(-21.0, year);
  898. }
  899.  
  900. /*
  901.  * passion_sunday:
  902.  * Christian holidays: compute Julian day for Passion Sunday
  903.  * (Second Sunday before Easter)
  904.  * given year (>1583)
  905.  */
  906. double
  907. passion_sunday(year)
  908.     int    year;
  909. {
  910.     return easter_offset(-14.0, year);
  911. }
  912.  
  913. /*
  914.  * palm_sunday:
  915.  * Christian holidays: compute Julian day for Palm Sunday
  916.  * (Sunday before Easter)
  917.  * given year (>1583)
  918.  */
  919. double
  920. palm_sunday(year)
  921.     int    year;
  922. {
  923.     return easter_offset(-7.0, year);
  924. }
  925.  
  926. /*
  927.  * maundy_thursday:
  928.  * Christian holidays: compute Julian day for Maundy Thursday
  929.  * (Three days prior to Easter)
  930.  * given year (>1583)
  931.  */
  932. double
  933. maundy_thursday(year)
  934.     int    year;
  935. {
  936.     return easter_offset(-3.0, year);
  937. }
  938.  
  939. /*
  940.  * good_friday:
  941.  * Christian holidays: compute Julian day for Good Friday
  942.  * (Two days prior to Easter)
  943.  * given year (>1583)
  944.  */
  945. double
  946. good_friday(year)
  947.     int    year;
  948. {
  949.     return easter_offset(-2.0, year);
  950. }
  951.  
  952. /*
  953.  * easter_monday:
  954.  * Christian holidays: compute Julian day for Easter Monday (Canada)
  955.  * given year (>1583)
  956.  */
  957. double
  958. easter_monday(year)
  959.     int    year;
  960. {
  961.     return easter_offset(1.0, year);
  962. }
  963.  
  964. /*
  965.  * rogation_sunday:
  966.  * Christian holidays: compute Julian day for Rogation Sunday
  967.  * (Rogation is the period of three days prior to Ascension Day; strictly
  968.  * speaking, this is the period Mon-Wed)
  969.  * given year (>1583)
  970.  */
  971. double
  972. rogation_sunday(year)
  973.     int    year;
  974. {
  975.     return easter_offset(35.0, year);
  976. }
  977.  
  978. /*
  979.  * ascension_day:
  980.  * Christian holidays: compute Julian day for Ascension Day
  981.  * (Holy Thursday; fortieth day after Easter, inclusive)
  982.  * given year (>1583)
  983.  */
  984. double
  985. ascension_day(year)
  986.     int    year;
  987. {
  988.     return easter_offset(39.0, year);
  989. }
  990.  
  991. /*
  992.  * whitsunday:
  993.  * Christian holidays: compute Julian day for Whitsunday (Pentecost)
  994.  * (Seventh Sunday after Easter)
  995.  * given year (>1583)
  996.  */
  997. double
  998. whitsunday(year)
  999.     int    year;
  1000. {
  1001.     return easter_offset(49.0, year);
  1002. }
  1003.  
  1004. /*
  1005.  * trinity_sunday:
  1006.  * Christian holidays: compute Julian day for Trinity Sunday
  1007.  * (Eighth Sunday after Easter)
  1008.  * given year (>1583)
  1009.  */
  1010. double
  1011. trinity_sunday(year)
  1012.     int    year;
  1013. {
  1014.     return easter_offset(56.0, year);
  1015. }
  1016.  
  1017. /*
  1018.  * corpus_christi:
  1019.  * Christian holidays: compute Julian day for Corpus Christi
  1020.  * (First Thursday after Trinity Sunday)
  1021.  * given year (>1583)
  1022.  */
  1023. double
  1024. corpus_christi(year)
  1025.     int    year;
  1026. {
  1027.     return easter_offset(60.0, year);
  1028. }
  1029.  
  1030. /*
  1031.  * passover:
  1032.  * Jewish holidays: compute Julian day of Pesach (first day of Passover)
  1033.  * and establish the Gregorian day of month and month, and Jewish year,
  1034.  * given Julian year (>1583)
  1035.  * This formula, due to Karl Friedrich Gauss (1777-1855) is described in
  1036.  * excruciating detail in the book by Schocken P. 51-61).  Note that it
  1037.  * computes the Julian date, which has to be corrected to a Gregorian date,
  1038.  * and that it exhibits the eccentricity of always computing the day
  1039.  * relative to March, so that April 2 appears as March 33, which is also
  1040.  * corrected below.
  1041.  *
  1042.  * Floating point implementation by R.P.C. Rodgers; integer implementation
  1043.  * (for faster calculation) by Amos Shapir (amos@nsc.com).
  1044.  */
  1045. void
  1046. passover_init(year)
  1047.     int    year;
  1048. {
  1049.     int    etmp, p_day;
  1050.     int    atmp, btmp, ctmp, day_of_week, dtmp, ftmp, gtmp;
  1051.     int    p_month;
  1052.  
  1053.     atmp = year + 3760;
  1054.     passoverJY = atmp;
  1055.     btmp = (12 * atmp + 17) % 19;
  1056.     ctmp = atmp % 4;
  1057.     etmp = (765433 * btmp) - (1565 * atmp)
  1058.         + (ctmp * 123120) + 15781075;
  1059.     dtmp = etmp / 492480;
  1060.     etmp %= 492480;
  1061.         /* day_of_week is not to be confused with the
  1062.          value returned by the day_of_week routine; here, Sunday = 1 */
  1063.     day_of_week = ((3 * atmp) + (5 * ctmp) + dtmp + 5) % 7;
  1064.     if (day_of_week == 0 && btmp > 11 && etmp >= 442111)
  1065.         p_day = dtmp + 1;
  1066.     else if (day_of_week == 1 && btmp > 6 && etmp >= 311676)
  1067.         p_day = dtmp + 2;
  1068.     else if (day_of_week == 2 || day_of_week == 4 || day_of_week ==6)
  1069.         p_day = dtmp + 1;
  1070.     else {
  1071.         p_day = dtmp;
  1072.     }
  1073.     ftmp = year / 100;        /* Correct to Gregorian date */
  1074.     gtmp = ((3 * ftmp) - 5) / 4;
  1075.     p_day += gtmp;
  1076.     if (p_day > 31) {        /* Correct for March days > 31 */
  1077.         p_day -= 31;
  1078.         p_month = 4;
  1079.         }
  1080.     else
  1081.         p_month = 3;
  1082.     passoverJD = julian_day((double)p_day, p_month, year);
  1083. }
  1084.  
  1085. double
  1086. passover(year, jyear)
  1087.     int year, *jyear;
  1088. {
  1089.     *jyear = passoverJY;
  1090.     return passoverJD;
  1091. }
  1092.  
  1093. /*
  1094.  * passover_offset:
  1095.  * Jewish holidays: compute Julian day as offset from Passover
  1096.  * given year (>1583)
  1097.  */
  1098. double
  1099. passover_offset(offset, year, jyear)
  1100.     double    offset;
  1101.     int    *jyear, year;
  1102. {
  1103.     *jyear = passoverJY;
  1104.     return passoverJD + offset;
  1105. }
  1106.  
  1107. /*
  1108.  * purim:
  1109.  * Jewish holidays: compute Julian day and Jewish year for Purim (Feast of Lots)
  1110.  * given year (>1583)
  1111.  */
  1112. double
  1113. purim(year, jyear)
  1114.     int    *jyear, year;
  1115. {
  1116.     return passover_offset(-30.0, year, jyear);
  1117. }
  1118.  
  1119. /*
  1120.  * shavuot:
  1121.  * Jewish holidays: compute Julian day and Jewish year for First day of Shavuot
  1122.  * given year (>1583)
  1123.  */
  1124. double
  1125. shavuot(year, jyear)
  1126.     int    *jyear, year;
  1127. {
  1128.     return passover_offset(50.0, year, jyear);
  1129. }
  1130.  
  1131. /*
  1132.  * rosh_hashanah:
  1133.  * Jewish holidays: compute Julian day and Jewish year for first day of
  1134.  * Rosh Hashanah (New Year) given year (>1583)
  1135.  */
  1136. double
  1137. rosh_hashanah(year, jyear)
  1138.     int    *jyear, year;
  1139. {
  1140.     double    atmp;
  1141.     atmp = passover_offset(163.0, year, jyear);
  1142.     (*jyear)++;
  1143.     return atmp;
  1144. }
  1145.  
  1146. /*
  1147.  * yom_kippur:
  1148.  * Jewish holidays: compute Julian day and Jewish year for Yom Kippur
  1149.  * given year (>1583)
  1150.  */
  1151. double
  1152. yom_kippur(year, jyear)
  1153.     int    *jyear, year;
  1154. {
  1155.     double    atmp;
  1156.     atmp = passover_offset(172.0, year, jyear);
  1157.     (*jyear)++;
  1158.     return atmp;
  1159. }
  1160.  
  1161. /*
  1162.  * sukkot:
  1163.  * Jewish holidays: compute Julian day and Jewish year for
  1164.  * First day of Sukkot (9 days) given year (>1583)
  1165.  */
  1166. double
  1167. sukkot(year, jyear)
  1168.     int    *jyear, year;
  1169. {
  1170.     double    atmp;
  1171.     atmp = passover_offset(177.0, year, jyear);
  1172.     (*jyear)++;
  1173.     return atmp;
  1174. }
  1175.  
  1176. /*
  1177.  * simchat_torah:
  1178.  * Jewish holidays: compute Julian day and Jewish year for Simchat Torah
  1179.  * (in Diapsora) given year (>1583)
  1180.  */
  1181. double
  1182. simchat_torah(year, jyear)
  1183.     int    *jyear, year;
  1184. {
  1185.     double    atmp;
  1186.     atmp = passover_offset(185.0, year, jyear);
  1187.     (*jyear)++;
  1188.     return atmp;
  1189. }
  1190.  
  1191. /*
  1192.  * chanukah:
  1193.  * Jewish holidays: compute Julian day and Jewish year for
  1194.  * first day of Chanukah (8 days) given year (>1583)
  1195.  */
  1196. double
  1197. chanukah(year, jyear)
  1198.     int    *jyear, year;
  1199. {
  1200.     double    atmp, ptmp;
  1201.     int    btmp, ytmp;
  1202.  
  1203.     atmp = passover(year, jyear);
  1204.     /* we need top compute passover for next year, so
  1205.      * save current info and restore when done
  1206.      */
  1207.     ptmp = passoverJD;
  1208.     ytmp = passoverJY;
  1209.     passover_init(year + 1);
  1210.     btmp = passoverJD - atmp;
  1211.     passoverJD = ptmp;
  1212.     passoverJY = ytmp;
  1213.     (*jyear)++;
  1214.     if (btmp == 355 || btmp == 385)
  1215.         return atmp + 247.0;
  1216.     else
  1217.         return atmp + 246.0;
  1218. }
  1219.  
  1220. /*
  1221.  * islamic_date:
  1222.  * Islamic holidays: compute Gregorian date(s) and Islamic year for any
  1223.  * Islamic day, given ISLAMIC day (1-30), ISLAMIC month(1-12),
  1224.  * and GREGORIAN year (>1583)
  1225.  * This is complicated by the fact that a given Islamic date can appear as
  1226.  * often as twice within the same Gregorian year.
  1227.  */
  1228. islamic_date(
  1229.     mday, mmonth, number_of_days, day1, month1, day2, month2, year, myear1,
  1230.     myear2)
  1231.     double    *day1, *day2, mday;
  1232.     int    mmonth, *month1, *month2, *myear1, *myear2,
  1233.         *number_of_days, year;
  1234. {
  1235.     double    day;
  1236.     int    count, month, nyear;
  1237.     double    islamic_to_julian();
  1238.  
  1239.     *myear1 = year - 621;        /* approx. >= Muslim year */
  1240.     nyear = year - 1;
  1241.     for (count = 0; nyear != year && count <= 100; count++) {
  1242.         gregorian_date(&day, &month, &nyear,
  1243.             islamic_to_julian(mday, mmonth, *myear1));
  1244.         *myear1 = *myear1 + (year - nyear);
  1245.     }
  1246.     if (nyear == year) {
  1247.         *day1 = day;
  1248.         *month1 = month;
  1249.         /*
  1250.          * See if there is a second occurrence in same Gregorian year
  1251.          */
  1252.         gregorian_date(&day, &month, &nyear,
  1253.             islamic_to_julian(mday, mmonth, (*myear1 + 1)));
  1254.         if (nyear == year) {
  1255.             *day2 = day;
  1256.             *month2 = month;
  1257.             *number_of_days = 2;
  1258.             *myear2 = *myear1 + 1;
  1259.         }
  1260.         else {
  1261.             *number_of_days = 1;
  1262.             gregorian_date(&day, &month, &nyear,
  1263.                 islamic_to_julian(mday, mmonth, (*myear1 - 1)));
  1264.             if (nyear == year) {
  1265.                 *day2 = day;
  1266.                 *month2 = month;
  1267.                 *number_of_days = 2;
  1268.                 *myear2 = *myear1 + 1;
  1269.             }
  1270.         }
  1271.     }
  1272. /*    else return -1; */
  1273.     return 0;
  1274. }
  1275.  
  1276. /*
  1277.  * islamic_to_julian:
  1278.  * Islamic holidays: compute Julian day for any Islamic day,
  1279.  * given ISLAMIC day (1-30), ISLAMIC month(1-12), and ISLAMIC year (>962)
  1280.  * Formula from Schocken (p. 66)
  1281.  */
  1282. double
  1283. islamic_to_julian(mday, mmonth, myear)
  1284.     double    mday;
  1285.     int    mmonth, myear;
  1286. {
  1287.     double    etmp, ftmp, jday;
  1288.     int    atmp, btmp, ctmp, dtmp, nyear;
  1289.     double    corrected_julian_day();
  1290.  
  1291.     nyear = myear + 621;        /* approx. Julian year */
  1292.     atmp = ((19 * myear) - 4) % 30;
  1293.     btmp = nyear % 4;
  1294.     ctmp = (mmonth - 1) / 2;
  1295.     dtmp = (mmonth - 1) % 2;
  1296.     etmp = mday + (59.0 * ctmp) + (30.0 * dtmp) + (atmp / 30.0)
  1297.          + (btmp / 4.0) - (10.8833333 * myear) + 146.8833333;
  1298.     if (etmp < 0.0) {
  1299.         ftmp = (4.0 * fabs(etmp)) / 1461.0;
  1300.         atmp = ftmp;
  1301.         dtmp = (4.0 * fabs(etmp));
  1302.         dtmp = dtmp % 1461;
  1303.         nyear -= (atmp + 1);
  1304.         ctmp = nyear % 4;
  1305.         jday = julian_day(1.0, 3, nyear)
  1306.             + (((1461.0 - dtmp) + ctmp - btmp) / 4.0)
  1307.             - 1.0;
  1308.         return corrected_julian_day(jday);
  1309.     }
  1310.     jday = julian_day(1.0, 3, nyear) + etmp - 1.0;
  1311.     return corrected_julian_day(jday);
  1312. }
  1313.  
  1314. /*
  1315.  * islamic_new_year:
  1316.  * Islamic holidays: compute Julian date(s) and Islamic year(s)
  1317.  * for Islamic New Year given JULIAN year (>1583)
  1318.  * (note: due to the length of the Islamic year (355 d), there can be portions
  1319.  * of as many as two Islamic New Years within a given Gregorian year; for
  1320.  * example, according to the algorithm below, the Islamic New Year occured
  1321.  * twice in 1975: Wednesday 1 January, Sunday 21 December)
  1322.  *
  1323.  * Algorithm: Schocken, page 66, which agrees with dates I have obtained for
  1324.  * the (Gregorian) years 1962-2000.  Schocken outlines a Muslim calendar
  1325.  * consisting of 12 months which are alternately 30 and 29 days in length
  1326.  * (the first month, Muharram, is 30 days in length).  In a 30-year cycle,
  1327.  * there are 11 leap years in which a 30th day is added to the final month
  1328.  * of the year.  See full details below.
  1329.  *
  1330.  * According to Dr. Omar Afzal of Cornell University ((607)255-5118,
  1331.  * 277-6707; Chairman of the Committee for Crescent Observation;
  1332.  * irfan@mvax.ee.cornell.edu), this is an antiquated system which has been
  1333.  * in use for some 400 years, but today the Muslim calendar follows a more
  1334.  * strictly lunar basis.  I wish to acknowledge the Pakistani Consular
  1335.  * Office of San Francisco (415)788-0677, who referred me to Dr. Muzammil
  1336.  * Siddiqui of the Orange County Islamic Center, (714)531-1722, who in turn
  1337.  * referred me to to Dr. Afzal.  I thank Dr. Afzal for his patient and lucid
  1338.  * explanations of the Islamic calendar, and for providing printed matter and
  1339.  * tables of dates to assist me.  Among the sources he provided:
  1340.  *
  1341.  * %A U. V. Tsybulsky
  1342.  * %B Calendar of Middle Eastern Countries
  1343.  * %I Nauka Publishing House
  1344.  * %C Moscow
  1345.  * %L English
  1346.  * %D 1979
  1347.  *
  1348.  * which provides a scholarly discussion of calendrical conventions in various
  1349.  * Muslim countries, as well as simple formulas for conversion between
  1350.  * Gregorian and Islamic dates.  It also includes translations of tables which
  1351.  * originally appeared in:
  1352.  *
  1353.  * %A F. R. Unat
  1354.  * %B Hicri Tarihleri
  1355.  * %I Turktarih Kurumu Basimevi
  1356.  * %C Ankara
  1357.  * %L Turkish
  1358.  * %D 1959
  1359.  *
  1360.  * Additional tabular material is available in:
  1361.  *
  1362.  * %A G. S. P. Freeman-Grenville
  1363.  * %B The Muslim and Christian Calendars
  1364.  * %I Oxford University Press
  1365.  * %C London
  1366.  * %D 1963
  1367.  *
  1368.  * The Islamic calendar is also known as the Hijri calendar, and dates often
  1369.  * have "A.H." or "a.h." appended to them to indicate "Anno Hijri."  A
  1370.  * listing of the months and their lengths (in the old scheme):
  1371.  *
  1372.  *    Month Name     Days         Days to add to 1 Muharram to get 1st of month
  1373.  * 1  Muharram        30          0
  1374.  * 2  Safar           29         30
  1375.  * 3  Rabi' al-Awwal  30         59
  1376.  * 4  Rabi' ath-Thani 29         89
  1377.  * 5  Jumada al-Ula   30        118
  1378.  * 6  Jumada al-Akhir 29        148
  1379.  * 7  Rajab           30        177
  1380.  * 8  Sha'ban         29        207
  1381.  * 9  Ramadan         30        236
  1382.  * 10 Shawwal         29        266
  1383.  * 11 Dhul-Qa'da      30        295
  1384.  * 12 Dhul-Hijja      29 (30)   325
  1385.  *
  1386.  * In the old scheme, used here for simplicity, adherence to this set of rules
  1387.  * led to a gradual lag of the month behind the actual crescent moon.  Thus,
  1388.  * a leap day was appended to the final month in each of 11 years out of every
  1389.  * 30 year cycle: (2, 5, 7, 10, 13, 16, 18, 21, 24, 26, 29).  If the Hijra year
  1390.  * is divided by 30, and the remainder is equal to one of these numbers, it is
  1391.  * a leap year.
  1392.  *
  1393.  * The Islamic day begins after sunset, and the beginning of the month is tied
  1394.  * to the first VISIBLE appearance of the crescent moon, which often lags behind
  1395.  * the astronomical new moon.  Also, the crescent was sometimes computed
  1396.  * according to Mecca time.  There are numerous methods for defining the date
  1397.  * of the crescent moon (and hence calendar dates) among the various Muslim
  1398.  * regions of the world.  Given these uncertainties, the dates computed here are
  1399.  * likely to be accurate to only +/- 2 days.
  1400.  *
  1401.  * Any errors in this code are my own fault, and
  1402.  * not intended to offend any members of the Islamic faith.
  1403.  */
  1404. islamic_new_year(year, number_of_dates, date1, date2, myear1, myear2)
  1405.     double    *date1, *date2;
  1406.     int    *number_of_dates, *myear1, *myear2, year;
  1407. {
  1408.     double    day1, day2;
  1409.     int    month1, month2;
  1410.  
  1411.     if (islamic_date(1.0, 1, number_of_dates, &day1, &month1, &day2,
  1412.         &month2, year, myear1, myear2) < 0) return -1;
  1413.     *date1 = julian_day(day1, month1, year);
  1414.     if (*number_of_dates == 2) *date2 = julian_day(day2, month2, year);
  1415.     return 0;
  1416. }
  1417.  
  1418. /*
  1419.  * islamic_offset:
  1420.  * Islamic holidays: compute Julian day(s) for an Islamic date as an offset
  1421.  * from the Islamic New Year, given (Gregorian) year
  1422.  */
  1423. islamic_offset(offset, year, number_of_dates, date1, date2, myear1, myear2)
  1424.     double    *date1, *date2, offset;
  1425.     int    *number_of_dates, *myear1, *myear2, year;
  1426. {
  1427.     double    day, tdate1, tdate2;
  1428.     int    month, tyear1, tyear2;
  1429.  
  1430.     (void) islamic_new_year(year, number_of_dates, date1, date2,
  1431.         myear1, myear2);
  1432.     if (*number_of_dates == 2) {
  1433.         tdate2 = *date2 + offset;
  1434.         gregorian_date(&day, &month, &tyear2, tdate2);
  1435.         if (tyear2 > year) {
  1436.             tdate2 = *date1 + offset;
  1437.             (void) islamic_new_year((year - 1), number_of_dates,
  1438.                 date1, date2, myear1, myear2);
  1439.             tdate1 = *date1 + offset;
  1440.             gregorian_date(&day, &month, &tyear1, tdate1);
  1441.             if (tyear1 == year) {
  1442.                 *date1 = tdate1;
  1443.                 *date2 = tdate2;
  1444.                 *number_of_dates = 2;
  1445.                 *myear2 = *myear1 + 1;
  1446.             }
  1447.             else {
  1448.                 *date1 = tdate2;
  1449.                 *number_of_dates = 1;
  1450.                 *myear1 = *myear1 + 1;
  1451.             }
  1452.         }
  1453.         else {
  1454.             *date1 += offset;
  1455.             *date2 = tdate2;
  1456.         }
  1457.     }
  1458.     else {
  1459.         tdate2 = *date1 + offset;
  1460.         gregorian_date(&day, &month, &tyear2, tdate2);
  1461.         if (tyear2 > year) {
  1462.             (void) islamic_new_year((year - 1), number_of_dates,
  1463.                 date1, date2, myear1, myear2);
  1464.             if (*number_of_dates == 2) {
  1465.                 *date1 = *date2 + offset;
  1466.                 *number_of_dates = 1;
  1467.                 *myear1 = *myear2;
  1468.             }
  1469.             else {
  1470.                 *date1 = *date1 + offset;
  1471.             }
  1472.         }
  1473.         else {
  1474.             (void) islamic_new_year((year - 1), number_of_dates,
  1475.                 date1, date2, myear1, myear2);
  1476.             if (*number_of_dates == 2) {
  1477.                 tdate1 = *date2 + offset;
  1478.                 *myear1 = *myear2;
  1479.             }
  1480.             else {
  1481.                 tdate1 = *date1 + offset;
  1482.             }
  1483.             gregorian_date(&day, &month, &tyear1, tdate1);
  1484.             if (tyear1 == year) {
  1485.                 *date1 = tdate1;
  1486.                 *date2 = tdate2;
  1487.                 *number_of_dates = 2;
  1488.                 *myear2 = *myear1 + 1;
  1489.             }
  1490.             else {
  1491.                 *date1 = tdate2;
  1492.                 *myear1 = *myear1 + 1;
  1493.                 *number_of_dates = 1;
  1494.             }
  1495.         }
  1496.     }
  1497. }
  1498.  
  1499. /*
  1500.  * muharram_9:
  1501.  * Islamic holidays: compute Julian day(s) and Islamic year(s) for
  1502.  * Muharram 9 (Day of fasting) given year (>1583)
  1503.  */
  1504. muharram_9(year, number_of_dates, date1, date2, myear1, myear2)
  1505.     double    *date1, *date2;
  1506.     int    *number_of_dates, *myear1, *myear2, year;
  1507. {
  1508.     islamic_offset(
  1509.         8.0, year, number_of_dates, date1, date2, myear1, myear2);
  1510. }
  1511.  
  1512. /*
  1513.  * muharram_10:
  1514.  * Islamic holidays: compute Julian day(s) and Islamic year(s) for
  1515.  * Muharram 10 (Day of deliverance of Moses from the Pharoah; for Shia Islam,
  1516.  * martyrdom of Husain) given year (>1583)
  1517.  */
  1518. muharram_10(year, number_of_dates, date1, date2, myear1, myear2)
  1519.     double    *date1, *date2;
  1520.     int    *number_of_dates, *myear1, *myear2, year;
  1521. {
  1522.     islamic_offset(
  1523.         9.0, year, number_of_dates, date1, date2, myear1, myear2);
  1524. }
  1525.  
  1526. /*
  1527.  * muharram_16:
  1528.  * Islamic holidays: compute Julian day(s) and Islamic year(s) for
  1529.  * Muharram 16 (Imamat Day; Ismaili Khoja) given year (>1583)
  1530.  */
  1531. muharram_16(year, number_of_dates, date1, date2, myear1, myear2)
  1532.     double    *date1, *date2;
  1533.     int    *number_of_dates, *myear1, *myear2, year;
  1534. {
  1535.     islamic_offset(
  1536.         15.0, year, number_of_dates, date1, date2, myear1, myear2);
  1537. }
  1538.  
  1539. /*
  1540.  * eid_i_milad_un_nabi:
  1541.  * Islamic holidays: compute Julian day(s) and Islamic year(s) for
  1542.  * Rabi I 12 (Eid-i-Milad-un-Nabi: The Prophet's Birthday) given year (>1583)
  1543.  */
  1544. eid_i_milad_un_nabi(year, number_of_dates, date1, date2, myear1, myear2)
  1545.     double    *date1, *date2;
  1546.     int    *number_of_dates, *myear1, *myear2, year;
  1547. {
  1548.     islamic_offset(
  1549.         70.0, year, number_of_dates, date1, date2, myear1, myear2);
  1550. }
  1551.  
  1552. /*
  1553.  * jumada_al_akhir_23:
  1554.  * Islamic holidays: compute Julian day(s) and Islamic year(s) for
  1555.  * Jumada al-Akhir 23 (Birth of Agha Khan IV, Ismaili) given year (>1583)
  1556.  */
  1557. jumada_al_akhir_23(year, number_of_dates, date1, date2, myear1, myear2)
  1558.     double    *date1, *date2;
  1559.     int    *number_of_dates, *myear1, *myear2, year;
  1560. {
  1561.     islamic_offset(
  1562.         170.0, year, number_of_dates, date1, date2, myear1, myear2);
  1563. }
  1564.  
  1565. /*
  1566.  * shab_e_miraj:
  1567.  * Islamic holidays: compute Julian day(s) and Islamic year(s) for
  1568.  * Rajab 27 (Shab-e-Mi'raj: The Prophet's Ascension) given year (>1583)
  1569.  */
  1570. shab_e_miraj(year, number_of_dates, date1, date2, myear1, myear2)
  1571.     double    *date1, *date2;
  1572.     int    *number_of_dates, *myear1, *myear2, year;
  1573. {
  1574.     islamic_offset(
  1575.         203.0, year, number_of_dates, date1, date2, myear1, myear2);
  1576. }
  1577.  
  1578. /*
  1579.  * shab_e_barat:
  1580.  * Islamic holidays: compute Julian day(s) and Islamic year(s) for
  1581.  * Shaban 15 (Shab-e-Bara't: Night, followed by day of fasting) given year (>1583)
  1582.  */
  1583. shab_e_barat(year, number_of_dates, date1, date2, myear1, myear2)
  1584.     double    *date1, *date2;
  1585.     int    *number_of_dates, *myear1, *myear2, year;
  1586. {
  1587.     islamic_offset(
  1588.         221.0, year, number_of_dates, date1, date2, myear1, myear2);
  1589. }
  1590.  
  1591. /*
  1592.  * ramadan:
  1593.  * Islamic holidays: compute Julian day(s) and Islamic year(s) for
  1594.  * Ramadan 1 (Fasting month begins) given year (>1583)
  1595.  */
  1596. ramadan(year, number_of_dates, date1, date2, myear1, myear2)
  1597.     double    *date1, *date2;
  1598.     int    *number_of_dates, *myear1, *myear2, year;
  1599. {
  1600.     islamic_offset(
  1601.         236.0, year, number_of_dates, date1, date2, myear1, myear2);
  1602. }
  1603.  
  1604. /*
  1605.  * shab_e_qadr:
  1606.  * Islamic holidays: compute Julian day(s) and Islamic year(s) for
  1607.  * Ramadan 27 (Shab-e Qadr: Night vigil) given year (>1583)
  1608.  */
  1609. shab_e_qadr(year, number_of_dates, date1, date2, myear1, myear2)
  1610.     double    *date1, *date2;
  1611.     int    *number_of_dates, *myear1, *myear2, year;
  1612. {
  1613.     islamic_offset(
  1614.         262.0, year, number_of_dates, date1, date2, myear1, myear2);
  1615. }
  1616.  
  1617. /*
  1618.  * eid_al_fitr:
  1619.  * Islamic holidays: compute Julian day(s) and Islamic year(s) for
  1620.  * Shawwal 1 (Eid-al-Fitr: Day of Feast) given year (>1583)
  1621.  */
  1622. eid_al_fitr(year, number_of_dates, date1, date2, myear1, myear2)
  1623.     double    *date1, *date2;
  1624.     int    *number_of_dates, *myear1, *myear2, year;
  1625. {
  1626.     islamic_offset(
  1627.         266.0, year, number_of_dates, date1, date2, myear1, myear2);
  1628. }
  1629.  
  1630. /*
  1631.  * dhul_hijja_9:
  1632.  * Islamic holidays: compute Julian day(s) and Islamic year(s) for
  1633.  * Dhul-Hijj 9 (Day of Pilgrimage at Arafat, Mecca) given year (>1583)
  1634.  */
  1635. dhul_hijja_9(year, number_of_dates, date1, date2, myear1, myear2)
  1636.     double    *date1, *date2;
  1637.     int    *number_of_dates, *myear1, *myear2, year;
  1638. {
  1639.     islamic_offset(
  1640.         333.0, year, number_of_dates, date1, date2, myear1, myear2);
  1641. }
  1642.  
  1643. /*
  1644.  * eid_al_adha:
  1645.  * Islamic holidays: compute Julian day(s) and Islamic year(s) for
  1646.  * Dhul-Hijj 10 (Eid-al-Adha: Day of Abraham's Sacrifice) given year (>1583)
  1647.  */
  1648. eid_al_adha(year, number_of_dates, date1, date2, myear1, myear2)
  1649.     double    *date1, *date2;
  1650.     int    *number_of_dates, *myear1, *myear2, year;
  1651. {
  1652.     islamic_offset(
  1653.         334.0, year, number_of_dates, date1, date2, myear1, myear2);
  1654. }
  1655.  
  1656. /*
  1657.  * ghadir:
  1658.  * Islamic holidays: compute Julian day(s) and Islamic year(s) for
  1659.  * Dhul-Hijj 18 (Ghadir: Ali's Nomination) given year (>1583)
  1660.  */
  1661. ghadir(year, number_of_dates, date1, date2, myear1, myear2)
  1662.     double    *date1, *date2;
  1663.     int    *number_of_dates, *myear1, *myear2, year;
  1664. {
  1665.     islamic_offset(
  1666.         342.0, year, number_of_dates, date1, date2, myear1, myear2);
  1667. }
  1668.  
  1669. /*
  1670.  * print_islamic_string:
  1671.  * Print formatted date string(s) of form:
  1672.  *    Monday 1 August 1988 (Islamic year 1409)
  1673.  * given string labelling date, (Gregorian) year of event, and function which
  1674.  * takes (Gregorian) year as an argument and produces: the (integer) number of
  1675.  * Gregorian dates (1 or 2), the first and second (Julian) days, and the first
  1676.  * and second Islamic years.
  1677.  */
  1678. print_islamic_string(string, year, func)
  1679.     char    *string;
  1680.     int    year;
  1681.     int    (*func) ();
  1682. {
  1683.  
  1684.     char    *date;
  1685.     double    date1, date2, day, holiday;
  1686.     int    month, myear1, myear2, number_of_dates;
  1687.     char    *date_string();
  1688.  
  1689.     holiday = (*func) (
  1690.         year, &number_of_dates, &date1, &date2, &myear1, &myear2);
  1691.     if (holiday < 0)
  1692.         (void) printf("Can not determine requested date\n");
  1693.     else {
  1694.         (void) printf("%s (%d) Julian day: %f\n", string, year, date1);
  1695.         gregorian_date(&day, &month, &year, date1);
  1696.         date = date_string(get_day_of_week(day, month, year),
  1697.             day, month, year);
  1698.         (void) printf("%s (%d): %s", string, year, date);
  1699.         (void) printf(" (Islamic year %d)\n", myear1);
  1700.         if (number_of_dates == 2) {
  1701.             (void) printf(
  1702.                 "%s (%d) Julian day: %f\n", string, year,
  1703.                 date2);
  1704.             gregorian_date(&day, &month, &year, date2);
  1705.             date = date_string(get_day_of_week(day, month, year),
  1706.                 day, month, year);
  1707.             (void) printf("%s (%d): %s", string, year, date);
  1708.             (void) printf(" (Islamic year %d)\n",
  1709.                 myear2);
  1710.         }
  1711.     }
  1712. }
  1713.  
  1714. /*
  1715.  * print_date_and_time_string:
  1716.  * Print formatted date/time string of form: 10:42 Thursday 8 December 1988
  1717.  * given string labelling date, year of event, and function which takes year
  1718.  * as an argument and produces Julian day, 
  1719.  */
  1720. print_date_and_time_string(string, year, func)
  1721.     char    *string;
  1722.     int    year;
  1723.     double    (*func)();
  1724. {
  1725.     char    *date;
  1726.     double    btmp, ctmp, day, holiday;
  1727.     int    atmp, hour, min, month;
  1728.     char    *date_time_string();
  1729.  
  1730.     holiday = (*func) (year);
  1731.     (void) printf("%s (%d) Julian day: %f\n", string, year, holiday);
  1732.     gregorian_date(&day, &month, &year, holiday);
  1733.     atmp = day;
  1734.     btmp = day - atmp;
  1735.     hour = btmp * 24.0;
  1736.     ctmp = (btmp - (hour / 24.0)) * 1440.0;
  1737.     min = ctmp;
  1738.     date = date_time_string(hour, min, get_day_of_week(day, month, year),
  1739.         day, month, year);
  1740.     (void) printf("%s (%d): %s\n", string, year, date);
  1741. }
  1742.  
  1743. /*
  1744.  * print_date_string:
  1745.  * Print formatted date string of form: Thursday 8 December 1988
  1746.  * given string labelling date, year of event, and function which takes year
  1747.  * as an argument and produces Julian day. 
  1748.  */
  1749. print_date_string(string, year, func)
  1750.     char    *string;
  1751.     int    year;
  1752.     double    (*func) ();
  1753. {
  1754.     char    *date;
  1755.     double    day, holiday;
  1756.     int    month;
  1757.     char    *date_string();
  1758.  
  1759.     holiday = (*func) (year);
  1760.     (void) printf("%s (%d) Julian day: %f\n", string, year, holiday);
  1761.     gregorian_date(&day, &month, &year, holiday);
  1762.     date = date_string(get_day_of_week(day, month, year), day, month, year);
  1763.     (void) printf("%s (%d): %s\n", string, year, date);
  1764. }
  1765.  
  1766. /*
  1767.  * print_jewish_string:
  1768.  * Print formatted date string of form: Thursday 8 December 1988 (NYEAR)
  1769.  * given string labelling date, year of event, and function which takes year
  1770.  * as an argument, sets the value of a non-Gregorian year (NYEAR), and produces
  1771.  * the Julian day.
  1772.  */
  1773. print_jewish_string(string, year, func)
  1774.     char    *string;
  1775.     int    year;
  1776.     double    (*func) ();
  1777. {
  1778.     char    *date;
  1779.     double    day, holiday;
  1780.     int    month, nyear;
  1781.     char    *date_string();
  1782.  
  1783.     holiday = (*func) (year, &nyear);
  1784.     (void) printf("%s (%d) Julian day: %f\n", string, year, holiday);
  1785.     gregorian_date(&day, &month, &year, holiday);
  1786.     date = date_string_2(get_day_of_week(day, month, year), day, month,
  1787.         year, nyear);
  1788.     (void) printf("%s (%d): %s\n", string, year, date);
  1789. }
  1790.  
  1791. /*
  1792.  * extract time of day from Julian day
  1793.  */
  1794. char *
  1795. julian_time(jday)
  1796. double jday;
  1797. {
  1798.     double    h1, floor();
  1799.     int    hours, minutes;
  1800.  
  1801.     jday += 0.5;  /* Julian day starts at noon GMT */
  1802.     h1 = (jday - floor(jday)) * 24.; /* number of hours & frac of hours */
  1803.     hours = (int) floor(h1);
  1804.     minutes = (int) ((h1 - (double) hours) * 60.);
  1805.     sprintf(timebuf, " at %d:%02d", hours, minutes);
  1806.     return(timebuf);
  1807. }
  1808.  
  1809. /*
  1810.  * End of code
  1811.  */
  1812. #endif    /* NO_HOLIDAYS */
  1813.