home *** CD-ROM | disk | FTP | other *** search
/ rtsi.com / 2014.01.www.rtsi.com.tar / www.rtsi.com / OS9 / OSK / EFFO / forum16.lzh / SOFTWARE / C / LUNISOLAR / lunisolar.c < prev    next >
Text File  |  1991-01-06  |  14KB  |  467 lines

  1. /* Prints the phase of the moon and generates LaTeX commands */
  2. /* that produce lunisolar calendars. */
  3. /* Usage: lunisolar
  4.    gives the phase of the moon,
  5.    and: lunisolar <year> <time_zone>
  6.    generates a lunisolar calendar for LaTeX.
  7.    Time zone may be one of:
  8.    GMT    NST    AST    EST    CST    
  9.    MST    PST    YST    HST    BST    
  10.    JST
  11. */
  12. /* Construct with the command "cc -O -o lunisolar lunisolar.c -lm". */
  13. /* John D. Ramsdell - May 1990 */
  14. static char copyright[] = "Copyright 1990 by The MITRE Corporation.";
  15. /*
  16.  * Permission to use, copy, modify, and distribute this
  17.  * software and its documentation for any purpose and without
  18.  * fee is hereby granted, provided that the above copyright
  19.  * notice appear in all copies.  The MITRE Corporation
  20.  * makes no representations about the suitability of this
  21.  * software for any purpose.  It is provided "as is" without
  22.  * express or implied warranty.
  23.  */
  24.  
  25. #include <stdio.h>
  26. #include <math.h>
  27. #include <time.h>
  28. #include <string.h>
  29.  
  30. #define HALF_PI (PI / 2.0)
  31. #define TWO_PI (2.0 * PI)
  32. #define RADIANS_PER_DEGREE (PI / 180.0)
  33. #define MINUTES_PER_DAY (60 * 24)
  34. #define SECONDS_PER_DAY (60 * MINUTES_PER_DAY)
  35. #define DAYS_PER_JULEAN_CENTURY 36525
  36.  
  37. typedef 
  38.   struct {
  39.     char *name;            /* Name of time zone. */
  40.     int offset;            /* Offset in minutes. */
  41.   } tz_t; 
  42.  
  43. tz_t tz_map[] =
  44. {
  45.   { "GMT",   0*60 },        /* Greenwich Mean Time */
  46.   { "NST",   7*30 },        /* Newfoundland is 3.5 hours */
  47.                 /* different from GMT. */
  48.   { "AST",   4*60 },        /* Alantic Standard Time. */
  49.   { "EST",   5*60 },        /* Eastern Standard Time. */
  50.   { "CST",   6*60 },        /* Central Standard Time. */
  51.   { "MST",   7*60 },        /* Mountain Standard Time. */
  52.   { "PST",   8*60 },        /* Pacific Standard Time. */
  53.   { "YST",   9*60 },        /* Yukon Standard Time. */
  54.   { "HST",  10*60 },        /* Hawaiian Standard Time. */
  55.   { "BST",  11*60 },        /* Bering Standard Time. */
  56.   { "JST",  -9*60 },        /* Japan Standard Time. */
  57.   { NULL ,   0    }        /* Mark end of list with NULL. */
  58. };
  59.  
  60. char *time_zone_name;        /* Selected time zone name. */
  61. double time_zone_offset;    /* Selected offset in minutes. */  
  62.  
  63. int select_time_zone (name)
  64. char *name;
  65. {
  66.   tz_t *t;
  67.   for (t = tz_map; t->name != NULL; t++)
  68.     if (strcmp (name, t->name) == 0) {
  69.       time_zone_name = name;
  70.       time_zone_offset = (double) t->offset;
  71.       return 0;            /* Found match. */
  72.     }
  73.   return 1;            /* No match found. */
  74. }
  75.  
  76. int leap_year (year)    /* True if year is a leap_year. */
  77. int year;
  78. {
  79.   return year % 4 == 0 && year % 100 != 0 || year % 400 == 0;
  80. }
  81.  
  82. /* Time is most often represented as a double precision number */
  83. /* in units of days.  Angles are in radians. */
  84.  
  85. /* J2000 is the date January 1, 2000; 12:00:00 GMT */
  86. /* This date is really called J2000.0. */
  87. struct tm J2000 = {0, 0, 12, 1, 0, 100, 6, 0, 0};
  88.  
  89. double days_after_J2000 ()    /* Returns the current time, */
  90. {                /* in units of days, after J2000.0. */
  91.   time_t current_seconds, J2000_seconds;
  92.   struct tm *current;
  93.   
  94.   if (time(¤t_seconds) == -1) { 
  95.     fprintf(stderr, "The current time not available.\n");
  96.     exit(1);
  97.   }
  98.   current = gmtime(¤t_seconds);
  99.   if (current == NULL) {
  100.     fprintf(stderr, "Coordinated Universal Time is not available.\n");
  101.     exit(1);
  102.   }
  103.   current_seconds = mktime (current);
  104.   if (current_seconds == -1) { 
  105.     fprintf(stderr, "Internal error -- cannot represent current time.\n");
  106.     exit (1);
  107.   }
  108.   J2000_seconds = mktime (&J2000);
  109.   if (J2000_seconds == -1) { 
  110.     fprintf(stderr, "Internal error -- cannot represent J2000.0.\n");
  111.     exit (1);
  112.   }
  113.   return difftime(current_seconds, J2000_seconds) / SECONDS_PER_DAY;
  114. }
  115.  
  116. double normalize_angle (angle)    /* Returns the angle between */
  117. double angle;
  118. {                    /* -PI < angle <= PI. */
  119.   if (angle > PI)
  120.     do angle -= TWO_PI; while (angle > PI);
  121.   else
  122.     while (angle <= -PI) angle += TWO_PI;
  123.   return angle;
  124. }
  125.  
  126. /*******************************************************************/
  127.  
  128. /* Astronomical almanac */
  129.  
  130. /*
  131.  * All formulas are from:
  132.  * The Astronomical Almanac for the Year 1984, 
  133.  * US Naval Observatory and Royal Greenwich Observatory,
  134.  * US Government Printing Office, Washington DC, 1984.
  135.  */
  136.  
  137. /* Angular position of the sun to a */
  138. /* precision of 0.01 degrees. (Page C24). */
  139.  
  140. #define SUN0 (RADIANS_PER_DEGREE * 280.460)
  141. #define SUN1 (RADIANS_PER_DEGREE *   0.9856474)
  142. #define SUN2 (RADIANS_PER_DEGREE * 357.528)
  143. #define SUN3 (RADIANS_PER_DEGREE *   0.9856003)
  144. #define SUN4 (RADIANS_PER_DEGREE *   1.915)
  145. #define SUN5 (RADIANS_PER_DEGREE *   0.020)
  146.  
  147. double sun_position (days)
  148. double days;
  149. {
  150.   double mean_longitude_of_sun, mean_anomaly, ecliptic_longitude;
  151.   mean_longitude_of_sun =
  152.     normalize_angle (SUN0 + SUN1 * days);
  153.   mean_anomaly =
  154.     normalize_angle (SUN2 + SUN3 * days);
  155.   ecliptic_longitude =
  156.     normalize_angle (mean_longitude_of_sun
  157.              + SUN4 * sin (mean_anomaly)
  158.              + SUN5 * sin (2.0 * mean_anomaly));
  159.   return ecliptic_longitude;
  160. }
  161.  
  162. /* Angular velocity of the sun.  Derivative of sun_position. */
  163.  
  164. double sun_velocity (days)
  165. double days;
  166. {
  167.   double mean_anomaly =
  168.     normalize_angle (SUN2 + SUN3 * days);
  169.   return SUN1 + SUN4 * SUN3 * cos (mean_anomaly)
  170.               + SUN5 * 2.0 * SUN3 * cos (2.0 * mean_anomaly);
  171. }
  172.     
  173. /* Angular position of the moon to a */
  174. /* precision of 0.3 degrees. (Page D46). */
  175.  
  176. #define RADIAN_CENTURY (RADIANS_PER_DEGREE / DAYS_PER_JULEAN_CENTURY)
  177.  
  178. #define MOON0  (RADIANS_PER_DEGREE * 218.32)
  179. #define MOON1  (RADIAN_CENTURY * 481267.883)
  180. #define MOON2A (RADIANS_PER_DEGREE *   6.29)
  181. #define MOON2B (RADIANS_PER_DEGREE * 134.9)
  182. #define MOON2C (RADIAN_CENTURY * 477198.85)
  183. #define MOON3A (RADIANS_PER_DEGREE *  -1.27)
  184. #define MOON3B (RADIANS_PER_DEGREE * 259.2)
  185. #define MOON3C (RADIAN_CENTURY * -413335.38)
  186. #define MOON4A (RADIANS_PER_DEGREE *   0.66)
  187. #define MOON4B (RADIANS_PER_DEGREE * 235.7)
  188. #define MOON4C (RADIAN_CENTURY * 890534.23)
  189. #define MOON5A (RADIANS_PER_DEGREE *   0.21)
  190. #define MOON5B (RADIANS_PER_DEGREE * 269.9)
  191. #define MOON5C (RADIAN_CENTURY * 954397.70)
  192. #define MOON6A (RADIANS_PER_DEGREE *  -0.19)
  193. #define MOON6B (RADIANS_PER_DEGREE * 357.5)
  194. #define MOON6C (RADIAN_CENTURY * 35999.05)
  195. #define MOON7A (RADIANS_PER_DEGREE *  -0.11)
  196. #define MOON7B (RADIANS_PER_DEGREE * 186.6)
  197. #define MOON7C (RADIAN_CENTURY * 966404.05)
  198.  
  199. double moon_position (days)
  200. double days;
  201. {
  202.  
  203. return normalize_angle (MOON0
  204.               + MOON1 * days
  205.               + MOON2A * sin (MOON2B + MOON2C * days)
  206.               + MOON3A * sin (MOON3B + MOON3C * days)
  207.               + MOON4A * sin (MOON4B + MOON4C * days)
  208.               + MOON5A * sin (MOON5B + MOON5C * days)
  209.               + MOON6A * sin (MOON6B + MOON6C * days)
  210.               + MOON7A * sin (MOON7B + MOON7C * days));
  211.  
  212. }
  213.  
  214. /****************************************************************/
  215.  
  216. /* Prints an English sentence giving the current phase of the moon. */
  217. #define PHASE_LIMIT MOON1
  218. int print_moon ()
  219. {
  220.   double days, phase;        /* Computes the moon's phase by */
  221.   int percent;            /* computing the difference between */
  222.   days = days_after_J2000 ();    /* the sun and moon's ecliptic longitude. */
  223.   phase = sun_position (days);
  224.   phase = normalize_angle (moon_position (days) - phase);
  225.   percent = 50.0 * (1.0 - cos (phase)) + 0.5; /* Visable fraction. */
  226.   printf("The moon is ");
  227.   if (fabs (phase) < PHASE_LIMIT)
  228.     printf ("new");
  229.   else if (fabs (normalize_angle (phase + PI)) < PHASE_LIMIT)
  230.     printf ("full");
  231.   else if (fabs (phase - HALF_PI) < PHASE_LIMIT)
  232.     printf ("first quarter (%d%% of full)", percent);
  233.   else if (fabs (phase + HALF_PI) < PHASE_LIMIT)
  234.     printf ("last quarter (%d%% of full)", percent);
  235.   else if (phase > HALF_PI)
  236.    printf ("waxing and gibbous (%d%% of full)", percent);
  237.   else if (phase > 0.0)
  238.    printf ("a waxing crescent (%d%% of full)", percent);
  239.   else if (phase > -HALF_PI)
  240.    printf ("a waning crescent (%d%% of full)", percent);
  241.   else
  242.    printf ("waning and gibbous (%d%% of full)", percent);
  243.   printf (".\n");
  244.   return 0;
  245. }
  246.  
  247. /**********************************************************/
  248.  
  249. /* lunisolar calendar routines. */
  250.  
  251. int first_day_of_year (year)/* Returns the integer number of days */
  252. int year;
  253. {                /* between the start of year and */    
  254.   int days = 365 * (year - 2000);/* J2000.0. */ 
  255.   if (year > 2000)
  256.     do { 
  257.       year--;
  258.       if (leap_year (year)) days++;
  259.     } while (year > 2000);
  260.   else
  261.     for (; year < 2000; year++)
  262.       if (leap_year (year)) days--;
  263.   return days;
  264. }
  265.  
  266. /* Routines that find the seasons. */
  267.  
  268. #define DIGITS 15
  269. /* Root finder using Newton's method. */
  270. int zero (x, f, fp)
  271. int x;
  272. double (*f)(),(*fp)();
  273. {
  274.   int i;
  275.   double y, midnite, noon;
  276.   y = x;
  277.   for (i = 0; i < DIGITS; i++)
  278.     y = y - f(y)/fp(y);
  279.   noon = 0.5 + time_zone_offset / MINUTES_PER_DAY;
  280.   midnite = floor (y - noon) + noon;
  281.   if (f (midnite) * f (midnite + 1.0) <= 0.0)
  282.     return midnite;
  283.   x = midnite;
  284.   printf ("%%Not sure about the season change for day %d.\n", x);
  285.   return x;
  286. }
  287.   
  288. double phase;            /* sun_zero has a root at the */
  289. double sun_zero (days)    /* desired day.  Used with zero */
  290. double days;
  291. {                /* to find the seasons. */
  292.   return normalize_angle (sun_position (days) - phase);
  293. }
  294.  
  295. void find_seasons (first_day, seasons)
  296. int first_day, *seasons;
  297. {                /* Remember Spring is the */
  298.   int i;            /* first season of a year. */
  299.   phase = -HALF_PI;        /* Find start of winters. */
  300.   seasons[0] = zero (first_day - 11, sun_zero, sun_velocity);
  301.   seasons[4] = zero (seasons[0] + 365, sun_zero, sun_velocity);
  302.   phase = 0.0;            /* Find start of other seasons. */
  303.   for (i = 1; i < 4; i++, phase += HALF_PI)
  304.     seasons[i] = zero (seasons[i-1] + 91, sun_zero, sun_velocity);
  305.   printf ("%% Seasons relative to January 1:");
  306.   for (i = 0; i < 5; i++)
  307.     printf (" %d", seasons[i] - first_day);
  308.   printf (".\n");
  309. }
  310.  
  311. /* Computes the position of the moon for each day at noon local time. */
  312. void make_moon_table (seasons, moon)
  313. int *seasons;
  314. float *moon;
  315. {
  316.   int i, day;
  317.   for (i = 0, day = seasons[0]; day < seasons[4]; i++, day++) {
  318.     double dday = day + time_zone_offset / MINUTES_PER_DAY;
  319.     moon[i] = normalize_angle (moon_position (dday) - sun_position (dday));
  320.   }
  321. }
  322.  
  323. /* Routines that output LaTeX commands. */
  324.  
  325. /* Dates spiral inward by an amount DELTA_RADIUS. */
  326. #define START_RADIUS 1.0
  327. #define DELTA_RADIUS 0.005
  328. float radius;
  329.  
  330. int month, day, moon_index;
  331. int days_per_month[12] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
  332.  
  333. /* Makes LaTeX statements that place the dates. */
  334. void mark_dates (year, from, to, moon) 
  335. int year, from, to;
  336. float *moon;
  337. {
  338.   radius = START_RADIUS;
  339.   for (; from < to; moon_index++, from++) {
  340.     printf ("\\put(%1.6f,%1.6f){\\makebox(0,0){%d/%d}}\n",
  341.         radius * sin (moon[moon_index]),
  342.         radius * cos (moon[moon_index]),
  343.         month, day);
  344.     radius -= DELTA_RADIUS;
  345.     day++;
  346.     if (day > days_per_month[month-1])
  347.       if (month == 2 && leap_year (year) && day == 29);
  348.       else {
  349.     day = 1;
  350.     month++;
  351.     if (month > 12) month = 1;
  352.       }
  353.   }
  354. }
  355.  
  356. void header (season, year)/* Start of each season. */
  357. char *season;
  358. int year;
  359. {
  360.   printf ("\\begin{figure}\n");
  361.   printf ("\\begin{center}\n");
  362.   printf ("\\begin{picture}(2.0,2.0)(-1.0,-1.0)\n");
  363.   printf ("\\tiny\n");
  364.   printf ("\\put(0,0){\\makebox(0,0){\\Huge %s %d}}\n",
  365.       season, year);
  366. }
  367.  
  368. void trailer ()        /* End of each season. */
  369. {
  370.   printf ("\\put(-1.0,0.0){\\line(1,0){0.5}}\n");
  371.   printf ("\\put(0.5,0.0){\\line(1,0){0.5}}\n");
  372.   printf ("\\put(0.0,-1.0){\\line(0,1){0.5}}\n");
  373.   printf ("\\put(0.0,-0.4){\\circle{0.1}}\n");
  374.   printf ("\\put(0.0,-0.3){\\makebox(0,0)[b]{\\large Full Moon}}\n");
  375.   printf ("\\put(0.0,0.5){\\line(0,1){0.5}}\n");
  376.   printf ("\\put(0.0,0.4){\\circle*{0.1}}\n");
  377.   printf ("\\put(0.0,0.3){\\makebox(0,0)[t]{\\large New Moon}}\n");
  378.   printf ("\\end{picture}\n");
  379.   printf ("\\\\ {\\Large Lunisolar Calendar}\n");
  380.   printf ("\\\\ {\\large Dates mark the lunar phase at noon %s.}\n",
  381.       time_zone_name);
  382.   printf ("\\end{center}\n");
  383.   printf ("\\end{figure}\n");
  384. }
  385.  
  386. char *season_titles[4] =
  387. { "Winter", "Spring", "Summer", "Fall"};
  388.  
  389. void LaTeXize_tables (year, first_day, seasons, moon)
  390. int year, first_day, *seasons;
  391. float *moon;
  392. {
  393.   int a_season;
  394.   printf ("\\documentstyle{article}\n");
  395.   printf ("\\pagestyle{empty}\n");
  396.   printf ("\\begin{document}\n");
  397.   printf ("\\Large\n");
  398.   printf ("\\setlength{\\unitlength}{60mm}\n");
  399.   month = 12;            /* December */
  400.   day = 32 - first_day + seasons[0];
  401.   moon_index = 0;
  402.   for (a_season = 0; a_season < 4; a_season++) {
  403.     header (season_titles[a_season], a_season == 0 ? year - 1 : year);
  404.     mark_dates (year, seasons[a_season], seasons[a_season+1], moon);
  405.     trailer ();
  406.   }
  407.   printf ("\\end{document}\n");
  408. }
  409.     
  410. /* Lunisolar master routine. */
  411.   
  412. int seasons[5];            /* Stores days that mark season changes. */
  413. float moon[370];        /* Stores moon phases for each day. */
  414.  
  415. /* Constructs a LaTeX file that generates a lunisolar calendar */
  416. /* for the year year and time zone tz. */
  417. int lunisolar (y, tz)
  418. char *y, *tz;
  419. {                
  420.   int year;
  421.   if (sscanf (y, "%d", &year) != 1) return 1;
  422.   if (year < 1950 || year > 2050) {
  423.     printf ("Program useful between the years 1950 and 2050.\n");
  424.     return 1;            /* error return. */
  425.   } 
  426.   else if (select_time_zone (tz) != 0)
  427.     return 1;
  428.   else {
  429.     int day_of_Jan1 = first_day_of_year (year);
  430.     printf ("%% Lunisolar calendar for %d.\n", year);
  431.     printf ("%% Constructed for %s, %1.2f hours %s of Greenwich.\n",
  432.         time_zone_name, fabs (time_zone_offset) / 60.0,
  433.         (time_zone_offset >= 0.0 ? "west" : "east"));
  434.     find_seasons (day_of_Jan1, seasons);
  435.     make_moon_table (seasons, moon);
  436.     LaTeXize_tables (year, day_of_Jan1, seasons, moon);
  437.     return 0;
  438.   }
  439. }
  440.  
  441. int main (argc, argv)/* Invokes print_moon with */
  442. int argc;
  443. char **argv;
  444. {                /* no arguments, and */
  445.   int i;            /* lunisolar with two. */
  446.   if (argc == 1) return print_moon ();
  447.   else if (argc == 3 && lunisolar (argv[1], argv[2]) == 0) return 0;
  448.   else {            /* print bad use error message. */
  449.     fprintf (stderr, "Bad args:");
  450.     for (i = 0; i < argc; i++)
  451.       fprintf (stderr, " %s", argv[i]);
  452.     fprintf (stderr,
  453.          "\nUsage: %s\ngives the phase of the moon,\n",
  454.          argv[0]);
  455.     fprintf (stderr, "and: %s <year> <time_zone>\n", argv[0]);
  456.     fprintf (stderr, "generates a lunisolar calendar for LaTeX.\n");
  457.     fprintf (stderr, "Time zone may be one of:");
  458.     for (i = 0; tz_map[i].name != NULL; i++) {
  459.       if (i % 5 == 0) fprintf (stderr, "\n");
  460.       fprintf (stderr, "%s\t", tz_map[i].name);
  461.     }
  462.     fprintf (stderr, "\n");
  463.     return 1;
  464.   }
  465. }
  466.  
  467.