ramsdell@linus.mitre.org (05/26/90)
Posting-number: Volume 13, Issue 4
Submitted-by: ramsdell@linus.mitre.org
Archive-name: lunisolar/part01
Enclosed is an ANSI C program for printing the phase of the moon and
generating LaTeX source which typesets a lunisolar calendar.
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: Wed May 23 08:44:37 1990
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 - May 1990 */
static const char copyright[] = "Copyright 1990 by The MITRE Corporation.";
/*
* 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.
*/
#define _ANSI_C_SOURCE
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
#define PI 3.141592653589793
#define HALF_PI (PI / 2.0)
#define TWO_PI (2.0 * PI)
#define RADIANS_PER_DEGREE (PI / 180.0)
#define MINUTES_PER_DAY (60 * 24)
#define SECONDS_PER_DAY (60 * MINUTES_PER_DAY)
#define DAYS_PER_JULEAN_CENTURY 36525
typedef
struct {
char *name; /* Name of time zone. */
int offset; /* Offset in minutes. */
} tz_t;
tz_t 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. */
{ NULL , 0 } /* Mark end of list with NULL. */
};
char *time_zone_name; /* Selected time zone name. */
double time_zone_offset; /* Selected offset in minutes. */
int select_time_zone (char *name)
{
tz_t *t;
for (t = tz_map; t->name != NULL; t++)
if (strcmp (name, t->name) == 0) {
time_zone_name = name;
time_zone_offset = (double) t->offset;
return 0; /* Found match. */
}
return 1; /* No match found. */
}
int leap_year (int year) /* True if year is a leap_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 date January 1, 2000; 12:00:00 GMT */
/* This date is really called J2000.0. */
struct tm J2000 = {0, 0, 12, 1, 0, 100, 6, 0, 0};
double days_after_J2000 (void) /* Returns the current time, */
{ /* in units of days, after J2000.0. */
time_t current_seconds, J2000_seconds;
struct tm *current;
if (time(¤t_seconds) == -1) {
fprintf(stderr, "The current time not available.\n");
exit(1);
}
current = gmtime(¤t_seconds);
if (current == NULL) {
fprintf(stderr, "Coordinated Universal Time is not available.\n");
exit(1);
}
current_seconds = mktime (current);
if (current_seconds == -1) {
fprintf(stderr, "Internal error -- cannot represent current time.\n");
exit (1);
}
J2000_seconds = mktime (&J2000);
if (J2000_seconds == -1) {
fprintf(stderr, "Internal error -- cannot represent J2000.0.\n");
exit (1);
}
return difftime(current_seconds, J2000_seconds) / SECONDS_PER_DAY;
}
double normalize_angle (double angle) /* Returns the angle between */
{ /* -PI < angle <= PI. */
if (angle > PI)
do angle -= TWO_PI; while (angle > PI);
else
while (angle <= -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 (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 (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 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 (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
int print_moon (void)
{
double days, phase; /* Computes the moon's phase by */
int percent; /* computing the difference between */
days = days_after_J2000 (); /* the sun and moon's 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 - HALF_PI) < PHASE_LIMIT)
printf ("first quarter (%d%% of full)", percent);
else if (fabs (phase + HALF_PI) < PHASE_LIMIT)
printf ("last quarter (%d%% of full)", percent);
else if (phase > HALF_PI)
printf ("waxing and gibbous (%d%% of full)", percent);
else if (phase > 0.0)
printf ("a waxing crescent (%d%% of full)", percent);
else if (phase > -HALF_PI)
printf ("a waning crescent (%d%% of full)", percent);
else
printf ("waning and gibbous (%d%% of full)", percent);
printf (".\n");
return 0;
}
/**********************************************************/
/* lunisolar calendar routines. */
int first_day_of_year (int year)/* Returns the integer number of days */
{ /* between the start of year and */
int days = 365 * (year - 2000);/* J2000.0. */
if (year > 2000)
do {
year--;
if (leap_year (year)) days++;
} while (year > 2000);
else
for (; year < 2000; year++)
if (leap_year (year)) days--;
return days;
}
/* Routines that find the seasons. */
#define DIGITS 15
/* Root finder using Newton's method. */
int zero (int x, double (*f)(double), double (*fp)(double))
{
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 / MINUTES_PER_DAY;
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 (double days) /* desired day. Used with zero */
{ /* to find the seasons. */
return normalize_angle (sun_position (days) - phase);
}
void find_seasons (int first_day, int *seasons)
{ /* Remember Spring is the */
int i; /* first season of a year. */
phase = -HALF_PI; /* 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 += HALF_PI)
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");
}
/* Computes the position of the moon for each day at noon local time. */
void make_moon_table (int *seasons, float *moon)
{
int i, day;
for (i = 0, day = seasons[0]; day < seasons[4]; i++, day++) {
double dday = day + time_zone_offset / MINUTES_PER_DAY;
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};
/* Makes LaTeX statements that place the dates. */
void mark_dates (int year, int from, int to, 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 (char *season, int year)/* Start of each season. */
{
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 (void) /* 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 (int year, int first_day, int *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. */
/* Constructs a LaTeX file that generates a lunisolar calendar */
/* for the year year and time zone tz. */
int lunisolar (char *y, char *tz)
{
int year;
if (sscanf (y, "%d", &year) != 1) return 1;
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;
}
}
int main (int argc, char **argv)/* Invokes print_moon with */
{ /* no arguments, and */
int i; /* lunisolar with two. */
if (argc == 1) return print_moon ();
else if (argc == 3 && lunisolar (argv[1], argv[2]) == 0) return 0;
else { /* print bad use error message. */
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; tz_map[i].name != NULL; i++) {
if (i % 5 == 0) fprintf (stderr, "\n");
fprintf (stderr, "%s\t", tz_map[i].name);
}
fprintf (stderr, "\n");
return 1;
}
}
SHAR_EOF
fi # end of overwriting check
# End of shell archive
exit 0