home *** CD-ROM | disk | FTP | other *** search
/ InfoMagic Source Code 1993 July / THE_SOURCE_CODE_CD_ROM.iso / bsd_srcs / games / pom / pom.c < prev    next >
Encoding:
C/C++ Source or Header  |  1991-04-08  |  5.4 KB  |  179 lines

  1. /*
  2.  * Copyright (c) 1989 The Regents of the University of California.
  3.  * All rights reserved.
  4.  *
  5.  * This code is derived from software posted to USENET.
  6.  *
  7.  * Redistribution and use in source and binary forms, with or without
  8.  * modification, are permitted provided that the following conditions
  9.  * are met:
  10.  * 1. Redistributions of source code must retain the above copyright
  11.  *    notice, this list of conditions and the following disclaimer.
  12.  * 2. Redistributions in binary form must reproduce the above copyright
  13.  *    notice, this list of conditions and the following disclaimer in the
  14.  *    documentation and/or other materials provided with the distribution.
  15.  * 3. All advertising materials mentioning features or use of this software
  16.  *    must display the following acknowledgement:
  17.  *    This product includes software developed by the University of
  18.  *    California, Berkeley and its contributors.
  19.  * 4. Neither the name of the University nor the names of its contributors
  20.  *    may be used to endorse or promote products derived from this software
  21.  *    without specific prior written permission.
  22.  *
  23.  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  24.  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  25.  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  26.  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  27.  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  28.  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  29.  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  30.  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  31.  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  32.  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  33.  * SUCH DAMAGE.
  34.  */
  35.  
  36. #ifndef lint
  37. char copyright[] =
  38. "@(#) Copyright (c) 1989 The Regents of the University of California.\n\
  39.  All rights reserved.\n";
  40. #endif /* not lint */
  41.  
  42. #ifndef lint
  43. static char sccsid[] = "@(#)pom.c    5.3 (Berkeley) 2/28/91";
  44. #endif /* not lint */
  45.  
  46. /*
  47.  * Phase of the Moon.  Calculates the current phase of the moon.
  48.  * Based on routines from `Practical Astronomy with Your Calculator',
  49.  * by Duffett-Smith.  Comments give the section from the book that
  50.  * particular piece of code was adapted from.
  51.  *
  52.  * -- Keith E. Brandt  VIII 1984
  53.  *
  54.  */
  55.  
  56. #include <sys/time.h>
  57. #include <stdio.h>
  58. #include <tzfile.h>
  59. #include <math.h>
  60.  
  61. #define    PI      3.141592654
  62. #define    EPOCH      85
  63. #define    EPSILONg  279.611371    /* solar ecliptic long at EPOCH */
  64. #define    RHOg      282.680403    /* solar ecliptic long of perigee at EPOCH */
  65. #define    ECCEN      0.01671542    /* solar orbit eccentricity */
  66. #define    lzero      18.251907    /* lunar mean long at EPOCH */
  67. #define    Pzero      192.917585    /* lunar mean long of perigee at EPOCH */
  68. #define    Nzero      55.204723    /* lunar mean long of node at EPOCH */
  69.  
  70. double dtor(), potm(), adj360();
  71.  
  72. main()
  73. {
  74.     extern int errno;
  75.     struct timeval tp;
  76.     struct timezone tzp;
  77.     struct tm *GMT, *gmtime();
  78.     double days, today, tomorrow;
  79.     int cnt;
  80.     char *strerror();
  81.  
  82.     if (gettimeofday(&tp,&tzp)) {
  83.         (void)fprintf(stderr, "pom: %s\n", strerror(errno));
  84.         exit(1);
  85.     }
  86.     GMT = gmtime(&tp.tv_sec);
  87.     days = (GMT->tm_yday + 1) + ((GMT->tm_hour +
  88.         (GMT->tm_min / 60.0) + (GMT->tm_sec / 3600.0)) / 24.0);
  89.     for (cnt = EPOCH; cnt < GMT->tm_year; ++cnt)
  90.         days += isleap(cnt) ? 366 : 365;
  91.     today = potm(days) + .5;
  92.     (void)printf("The Moon is ");
  93.     if ((int)today == 100)
  94.         (void)printf("Full\n");
  95.     else if (!(int)today)
  96.         (void)printf("New\n");
  97.     else {
  98.         tomorrow = potm(days + 1);
  99.         if ((int)today == 50)
  100.             (void)printf("%s\n", tomorrow > today ?
  101.                 "at the First Quarter" : "at the Last Quarter");
  102.         else {
  103.             (void)printf("%s ", tomorrow > today ?
  104.                 "Waxing" : "Waning");
  105.             if (today > 50)
  106.                 (void)printf("Gibbous (%1.0f%% of Full)\n",
  107.                     today);
  108.             else if (today < 50)
  109.                 (void)printf("Crescent (%1.0f%% of Full)\n",
  110.                     today);
  111.         }
  112.     }
  113. }
  114.  
  115. /*
  116.  * potm --
  117.  *    return phase of the moon
  118.  */
  119. double
  120. potm(days)
  121.     double days;
  122. {
  123.     double N, Msol, Ec, LambdaSol, l, Mm, Ev, Ac, A3, Mmprime;
  124.     double A4, lprime, V, ldprime, D, Nm;
  125.  
  126.     N = 360 * days / 365.2422;                /* sec 42 #3 */
  127.     adj360(&N);
  128.     Msol = N + EPSILONg - RHOg;                /* sec 42 #4 */
  129.     adj360(&Msol);
  130.     Ec = 360 / PI * ECCEN * sin(dtor(Msol));        /* sec 42 #5 */
  131.     LambdaSol = N + Ec + EPSILONg;                /* sec 42 #6 */
  132.     adj360(&LambdaSol);
  133.     l = 13.1763966 * days + lzero;                /* sec 61 #4 */
  134.     adj360(&l);
  135.     Mm = l - (0.1114041 * days) - Pzero;            /* sec 61 #5 */
  136.     adj360(&Mm);
  137.     Nm = Nzero - (0.0529539 * days);            /* sec 61 #6 */
  138.     adj360(&Nm);
  139.     Ev = 1.2739 * sin(dtor(2*(l - LambdaSol) - Mm));    /* sec 61 #7 */
  140.     Ac = 0.1858 * sin(dtor(Msol));                /* sec 61 #8 */
  141.     A3 = 0.37 * sin(dtor(Msol));
  142.     Mmprime = Mm + Ev - Ac - A3;                /* sec 61 #9 */
  143.     Ec = 6.2886 * sin(dtor(Mmprime));            /* sec 61 #10 */
  144.     A4 = 0.214 * sin(dtor(2 * Mmprime));            /* sec 61 #11 */
  145.     lprime = l + Ev + Ec - Ac + A4;                /* sec 61 #12 */
  146.     V = 0.6583 * sin(dtor(2 * (lprime - LambdaSol)));    /* sec 61 #13 */
  147.     ldprime = lprime + V;                    /* sec 61 #14 */
  148.     D = ldprime - LambdaSol;                /* sec 63 #2 */
  149.     return(50 * (1 - cos(dtor(D))));            /* sec 63 #3 */
  150. }
  151.  
  152. /*
  153.  * dtor --
  154.  *    convert degrees to radians
  155.  */
  156. double
  157. dtor(deg)
  158.     double deg;
  159. {
  160.     return(deg * PI / 180);
  161. }
  162.  
  163. /*
  164.  * adj360 --
  165.  *    adjust value so 0 <= deg <= 360
  166.  */
  167. double
  168. adj360(deg)
  169.     double *deg;
  170. {
  171.     for (;;)
  172.         if (*deg < 0)
  173.             *deg += 360;
  174.         else if (*deg > 360)
  175.             *deg -= 360;
  176.         else
  177.             break;
  178. }
  179.