Date: Thu, 5 Nov 87 09:45:56 EST
From: ramsdell%linus@mitre-bedford.ARPA
To: texhax@score.stanford.edu, rs@uunet.UU.NET
Subject: Lunisolar calendar generator for LaTeX

Here is a program that prints the current phase of the moon.  It also
generates the LaTeX commands that produce lunisolar calendars.  It is
a handy program for astronomers and sailors, as a lunisolar calendar
tells when the moon will be bright.  
John


#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create the files:
#	lunisolar.c
# This archive created: Thu Nov  5 08:53:57 1987
export PATH; PATH=/bin:$PATH
if test -f 'lunisolar.c'
then
	echo shar: will not over-write existing file "'lunisolar.c'"
else
cat << \SHAR_EOF > 'lunisolar.c'
/* Prints the phase of the moon and generates LaTeX commands */
/* that produce lunisolar calendars. */
/* Usage: lunisolar
   gives the phase of the moon,
   and: lunisolar <year> <time_zone>
   generates a lunisolar calendar for LaTeX.
   Time zone may be one of:
   GMT	NST	AST	EST	CST	
   MST	PST	YST	HST	BST	
   JST
*/
/* Construct with the command "cc -O -o lunisolar lunisolar.c -lm". */
/* John D. Ramsdell - November 1987 */
static char copyright[] =
  "Copyright 1987 by The MITRE Corporation.  All rights reserved.";
/*
 * Permission to use, copy, modify, and distribute this
 * software and its documentation for any purpose and without
 * fee is hereby granted, provided that the above copyright
 * notice appear in all copies.  The MITRE Corporation
 * makes no representations about the suitability of this
 * software for any purpose.  It is provided "as is" without
 * express or implied warranty.
 */

#include <stdio.h>
#include <math.h>

#define PI 3.141592653589793
#define NEGATIVE_PI (-1.0 * PI)
#define TWO_PI (2.0 * PI)
#define RADIANS_PER_DEGREE (PI / 180.0)

struct {
  char *name;			/* Name of time zone. */
  int offset;			/* Offset in minutes. */
} tz_map[] =
{
  { "GMT",   0*60 },		/* Greenwich Mean Time */
  { "NST",   7*30 },		/* Newfoundland is 3.5 hours */
				/* different from GMT. */
  { "AST",   4*60 },		/* Alantic Standard Time. */
  { "EST",   5*60 },		/* Eastern Standard Time. */
  { "CST",   6*60 },		/* Central Standard Time. */
  { "MST",   7*60 },		/* Mountain Standard Time. */
  { "PST",   8*60 },		/* Pacific Standard Time. */
  { "YST",   9*60 },		/* Yukon Standard Time. */
  { "HST",  10*60 },		/* Hawaiian Standard Time. */
  { "BST",  11*60 },		/* Bering Standard Time. */
  { "JST",  -9*60 },		/* Japan Standard Time. */
  { ""   ,   0    }		/* Mark end of list with "". */
};

char *time_zone_name;		/* Selected time zone name. */
double time_zone_offset;	/* Selected offset in minutes. */  

int select_time_zone (name)
     char *name;
{
  int i;
  if (strcmp (name, "") == 0) return 1;	/* Return error. */
  for (i = 0; strcmp (tz_map[i].name, "") != 0; i++)
    if (strcmp (name, tz_map[i].name) == 0) {
      time_zone_name = name;
      time_zone_offset = (double) tz_map[i].offset;
      return 0;			/* Found match. */
    }
  return 1;			/* No match found. */
}

int leap_year (year)		/* True if year is a leap_year. */
     int year;
{
  return year % 4 == 0 && year % 100 != 0 || year % 400 == 0;
}

/* Time is most often represented as a double precision number */
/* in units of days.  Angles are in radians. */

/* J2000 is the number of days between January 1, 2000; 12:00:00 GMT */
/* and the local origin of time. */
double J2000;			/* This date is really called J2000.0. */

void make_J2000 ()		/* Assumes time origin of */
{				/* January 1, 1970; 00:00:00 GMT. */
  int year;
  J2000 = (2000 - 1970) * 365 + 0.5;
  for (year = 1970; year < 2000; year++) /* Account for leap years. */
    if (leap_year (year)) J2000 += 1.0;
}

double days_after_J2000 ()	/* Returns the current time, */
{				/* in units of days, after J2000.0. */
  long seconds = time((long *) 0);
  double seconds_per_day = 24 * 60 * 60;
  return seconds / seconds_per_day - J2000;
}

double normalize_angle (angle)	/* Returns the angle between */
     double angle;		/* -PI < angle <= PI. */
{
  if (angle > PI)
    do angle -= TWO_PI; while (angle > PI);
  else
    while (angle <= NEGATIVE_PI) angle += TWO_PI;
  return angle;
}

/*******************************************************************/

/* Astronomical almanac */

/*
 * All formulas are from:
 * The Astronomical Almanac for the Year 1984, 
 * US Naval Observatory and Royal Greenwich Observatory,
 * US Government Printing Office, Washington DC, 1984.
 */

/* Angular position of the sun to a */
/* precision of 0.01 degrees. (Page C24). */

#define SUN0 (RADIANS_PER_DEGREE * 280.460)
#define SUN1 (RADIANS_PER_DEGREE *   0.9856474)
#define SUN2 (RADIANS_PER_DEGREE * 357.528)
#define SUN3 (RADIANS_PER_DEGREE *   0.9856003)
#define SUN4 (RADIANS_PER_DEGREE *   1.915)
#define SUN5 (RADIANS_PER_DEGREE *   0.020)

double sun_position (days)
     double days;
{
  double mean_longitude_of_sun, mean_anomaly, ecliptic_longitude;
  mean_longitude_of_sun =
    normalize_angle (SUN0 + SUN1 * days);
  mean_anomaly =
    normalize_angle (SUN2 + SUN3 * days);
  ecliptic_longitude =
    normalize_angle (mean_longitude_of_sun
		     + SUN4 * sin (mean_anomaly)
		     + SUN5 * sin (2.0 * mean_anomaly));
  return ecliptic_longitude;
}

/* Angular velocity of the sun.  Derivative of sun_position. */

double sun_velocity (days)
     double days;
{
  double mean_anomaly =
    normalize_angle (SUN2 + SUN3 * days);
  return SUN1 + SUN4 * SUN3 * cos (mean_anomaly)
              + SUN5 * 2.0 * SUN3 * cos (2.0 * mean_anomaly);
}
    
/* Angular position of the moon to a */
/* precision of 0.3 degrees. (Page D46). */

#define DAYS_PER_JULEAN_CENTURY 36525.0
#define RADIAN_CENTURY (RADIANS_PER_DEGREE / DAYS_PER_JULEAN_CENTURY)

#define MOON0  (RADIANS_PER_DEGREE * 218.32)
#define MOON1  (RADIAN_CENTURY * 481267.883)
#define MOON2A (RADIANS_PER_DEGREE *   6.29)
#define MOON2B (RADIANS_PER_DEGREE * 134.9)
#define MOON2C (RADIAN_CENTURY * 477198.85)
#define MOON3A (RADIANS_PER_DEGREE *  -1.27)
#define MOON3B (RADIANS_PER_DEGREE * 259.2)
#define MOON3C (RADIAN_CENTURY * -413335.38)
#define MOON4A (RADIANS_PER_DEGREE *   0.66)
#define MOON4B (RADIANS_PER_DEGREE * 235.7)
#define MOON4C (RADIAN_CENTURY * 890534.23)
#define MOON5A (RADIANS_PER_DEGREE *   0.21)
#define MOON5B (RADIANS_PER_DEGREE * 269.9)
#define MOON5C (RADIAN_CENTURY * 954397.70)
#define MOON6A (RADIANS_PER_DEGREE *  -0.19)
#define MOON6B (RADIANS_PER_DEGREE * 357.5)
#define MOON6C (RADIAN_CENTURY * 035999.05)
#define MOON7A (RADIANS_PER_DEGREE *  -0.11)
#define MOON7B (RADIANS_PER_DEGREE * 186.6)
#define MOON7C (RADIAN_CENTURY * 966404.05)

double moon_position (days)
     double days;
{
  return normalize_angle (MOON0
			  + MOON1 * days
			  + MOON2A * sin (MOON2B + MOON2C * days)
			  + MOON3A * sin (MOON3B + MOON3C * days)
			  + MOON4A * sin (MOON4B + MOON4C * days)
			  + MOON5A * sin (MOON5B + MOON5C * days)
			  + MOON6A * sin (MOON6B + MOON6C * days)
			  + MOON7A * sin (MOON7B + MOON7C * days));
}

/****************************************************************/

/* Prints an English sentence giving the current phase of the moon. */
#define PHASE_LIMIT MOON1
void print_moon ()
{
  double days, phase;		/* Computes the moon's phase by */
  int percent;			/* computing the difference between */
  make_J2000 ();		/* the sun and moon's */
  days = days_after_J2000 ();	/* ecliptic longitude. */
  phase = sun_position (days);
  phase = normalize_angle (moon_position (days) - phase);
  percent = 50.0 * (1.0 - cos (phase)) + 0.5; /* Visable fraction. */
  printf("The moon is ");
  if (fabs (phase) < PHASE_LIMIT)
    printf ("new");
  else if (fabs (normalize_angle (phase + PI)) < PHASE_LIMIT)
    printf ("full");
  else if (fabs (phase - PI/2.0) < PHASE_LIMIT)
    printf ("first quarter (%d%% of full)", percent);
  else if (fabs (phase + PI/2.0) < PHASE_LIMIT)
    printf ("last quarter (%d%% of full)", percent);
  else if (phase > PI/2.0)
   printf ("waxing and gibbous (%d%% of full)", percent);
  else if (phase > 0.0)
   printf ("a waxing crescent (%d%% of full)", percent);
  else if (phase > PI/-2.0)
   printf ("a waning crescent (%d%% of full)", percent);
  else
   printf ("waning and gibbous (%d%% of full)", percent);
  printf (".\n");
}

/**********************************************************/

/* lunisolar calendar routines. */

int first_day_of_year (year)	/* Returns the integer number of days */
     int year;			/* between the start of year and */
{				/* J2000.0. */
  int days;
  days = 365 * (year - 2000);
  if (year > 2000)
    for (; year > 2000; year--)
      if (leap_year (year)) days += 1;
      else;			/*  Needed else! */
  else
    for (; year < 2000; year++)
      if (leap_year (year)) days -= 1;
  return days;
}

/* Routines that find the seasons. */

#define DIGITS 15
int zero (x, f, fp)		/* Root finder using */
     int x;			/* Newton's method. */
     double (*f) ();
     double (*fp) ();
{
  int i;
  double y, midnite, noon;
  y = x;
  for (i = 0; i < DIGITS; i++)
    y = y - f(y)/fp(y);
  noon = 0.5 + (time_zone_offset / 24.0 / 60.0);
  midnite = floor (y - noon) + noon;
  if (f (midnite) * f (midnite + 1.0) <= 0.0)
    return midnite;
  x = midnite;
  printf ("%%Not sure about the season change for day %d.\n", x);
  return x;
}
  
double phase;			/* sun_zero has a root at the */
double sun_zero (days)		/* desired day.  Used with zero */
     double days;		/* to find the seasons. */
{
  return normalize_angle (sun_position (days) - phase);
}

void find_seasons (first_day, seasons)
     int first_day, *seasons;
{				/* Remember Spring is the */
  int i;			/* first season of a year. */
  phase = PI/-2.0;		/* Find start of winters. */
  seasons[0] = zero (first_day - 11, sun_zero, sun_velocity);
  seasons[4] = zero (seasons[0] + 365, sun_zero, sun_velocity);
  phase = 0.0;			/* Find start of other seasons. */
  for (i = 1; i < 4; i++, phase += PI/2.0)
    seasons[i] = zero (seasons[i-1] + 91, sun_zero, sun_velocity);
  printf ("%% Seasons relative to January 1:");
  for (i = 0; i < 5; i++)
    printf (" %d", seasons[i] - first_day);
  printf (".\n");
}

void make_moon_table (seasons, moon) /* Computes the position of */
     int *seasons;		/* the moon for each day at */
     float *moon;		/* noon local time. */
{
  int i, day;
  for (i = 0, day = seasons[0]; day < seasons[4]; i++, day++) {
    double dday = day + time_zone_offset / (24.0 * 60.0);
    moon[i] = normalize_angle (moon_position (dday) - sun_position (dday));
  }
}

/* Routines that output LaTeX commands. */

/* Dates spiral inward by an amount DELTA_RADIUS. */
#define START_RADIUS 1.0
#define DELTA_RADIUS 0.005
float radius;

int month, day, moon_index;
int days_per_month[12] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};

void mark_dates (year, from, to, moon) /* Makes LaTeX statements */
     int year, from, to;	/* that place the dates. */
     float *moon;
{
  radius = START_RADIUS;
  for (; from < to; moon_index++, from++) {
    printf ("\\put(%1.6f,%1.6f){\\makebox(0,0){%d/%d}}\n",
	    radius * sin (moon[moon_index]),
	    radius * cos (moon[moon_index]),
	    month, day);
    radius -= DELTA_RADIUS;
    day++;
    if (day > days_per_month[month-1])
      if (month == 2 && leap_year (year) && day == 29);
      else {
	day = 1;
	month++;
	if (month > 12) month = 1;
      }
  }
}

void header (season, year)	/* Start of each season. */
     char *season;
     int year;
{
  printf ("\\begin{figure}\n");
  printf ("\\begin{center}\n");
  printf ("\\begin{picture}(2.0,2.0)(-1.0,-1.0)\n");
  printf ("\\tiny\n");
  printf ("\\put(0,0){\\makebox(0,0){\\Huge %s %d}}\n",
	  season, year);
}

void trailer ()			/* End of each season. */
{
  printf ("\\put(-1.0,0.0){\\line(1,0){0.5}}\n");
  printf ("\\put(0.5,0.0){\\line(1,0){0.5}}\n");
  printf ("\\put(0.0,-1.0){\\line(0,1){0.5}}\n");
  printf ("\\put(0.0,-0.4){\\circle{0.1}}\n");
  printf ("\\put(0.0,-0.3){\\makebox(0,0)[b]{\\large Full Moon}}\n");
  printf ("\\put(0.0,0.5){\\line(0,1){0.5}}\n");
  printf ("\\put(0.0,0.4){\\circle*{0.1}}\n");
  printf ("\\put(0.0,0.3){\\makebox(0,0)[t]{\\large New Moon}}\n");
  printf ("\\end{picture}\n");
  printf ("\\\\ {\\Large Lunisolar Calendar}\n");
  printf ("\\\\ {\\large Dates mark the lunar phase at noon %s.}\n",
	  time_zone_name);
  printf ("\\end{center}\n");
  printf ("\\end{figure}\n");
}

char *season_titles[4] =
{ "Winter", "Spring", "Summer", "Fall"};

void LaTeXize_tables (year, first_day, seasons, moon)
     int year, first_day, *seasons;
     float *moon;
{
  int a_season;
  printf ("\\documentstyle{article}\n");
  printf ("\\pagestyle{empty}\n");
  printf ("\\begin{document}\n");
  printf ("\\Large\n");
  printf ("\\setlength{\\unitlength}{60mm}\n");
  month = 12;			/* December */
  day = 32 - first_day + seasons[0];
  moon_index = 0;
  for (a_season = 0; a_season < 4; a_season++) {
    header (season_titles[a_season], a_season == 0 ? year - 1 : year);
    mark_dates (year, seasons[a_season], seasons[a_season+1], moon);
    trailer ();
  }
  printf ("\\end{document}\n");
}
    
/* Lunisolar master routine. */
  
int seasons[5];			/* Stores days that mark season changes. */
float moon[370];		/* Stores moon phases for each day. */

int lunisolar (year, tz)	/* Constructs a LaTeX file that  */
     int year; char *tz;	/* generates a lunisolar calendar */
{				/* for the year year and time zone tz. */
  if (year < 1950 || year > 2050) {
    printf ("Program useful between the years 1950 and 2050.\n");
    return 1;			/* error return. */
  } 
  else if (select_time_zone (tz) != 0)
    return 1;
  else {
    int day_of_Jan1 = first_day_of_year (year);
    printf ("%% Lunisolar calendar for %d.\n", year);
    printf ("%% Constructed for %s, %1.2f hours %s of Greenwich.\n",
	    time_zone_name, fabs (time_zone_offset) / 60.0,
	    (time_zone_offset >= 0.0 ? "west" : "east"));
    find_seasons (day_of_Jan1, seasons);
    make_moon_table (seasons, moon);
    LaTeXize_tables (year, day_of_Jan1, seasons, moon);
    return 0;
  }
}

main (argc, argv)		/* Invokes print_moon with */
     int argc;			/* no arguments, and */
     char **argv;		/* lunisolar with one. */
{
  int i;
  if (argc == 1) print_moon ();
  else {
    if (argc == 3) {
      int year;
      if (sscanf (argv[1], "%d", &year) == 1 &&
	  lunisolar (year, argv[2]) == 0)
	exit (0);
    }
    fprintf (stderr, "Bad args:");
    for (i = 0; i < argc; i++)
      fprintf (stderr, " %s", argv[i]);
    fprintf (stderr,
	     "\nUsage: %s\ngives the phase of the moon,\n",
	     argv[0]);
    fprintf (stderr, "and: %s <year> <time_zone>\n", argv[0]);
    fprintf (stderr, "generates a lunisolar calendar for LaTeX.\n");
    fprintf (stderr, "Time zone may be one of:");
    for (i = 0; strcmp(tz_map[i].name, "") != 0; i++) {
      if (i % 5 == 0) fprintf (stderr, "\n");
      fprintf (stderr, "%s\t", tz_map[i].name);
    }
    fprintf (stderr, "\n");
    exit (1);
  }
}
SHAR_EOF
fi # end of overwriting check
#	End of shell archive
exit 0
-------
-------