home *** CD-ROM | disk | FTP | other *** search
/ Club Amiga de Montreal - CAM / CAM_CD_1.iso / files / 383.lha / Planets / moonphase.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-05-25  |  4.0 KB  |  178 lines
  1. /****************************************************************************
  2.  pom.c
  3.  
  4.      Phase of the Moon. Calculates the current phase of the moon.
  5.      Based on routines from `Practical Astronomy with Your Calculator',
  6.         by Duffett-Smith.
  7.      Comments give the section from the book that particular piece
  8.         of code was adapted from.
  9.  
  10.      -- Keith E. Brandt  VIII 1984
  11.  
  12.  ****************************************************************************/
  13.  
  14. #include <stdio.h>
  15. #include <sys/time.h>
  16. #include <math.h>
  17. #define PI         3.141592654
  18. #define EPOCH   1983
  19. #define EPSILONg 279.103035     /* solar ecliptic long at EPOCH */
  20. #define RHOg     282.648015     /* solar ecliptic long of perigee at EPOCH */
  21. #define e          0.01671626   /* solar orbit eccentricity */
  22. #define lzero    106.306091     /* lunar mean long at EPOCH */
  23. #define Pzero    111.481526     /* lunar mean long of perigee at EPOCH */
  24. #define Nzero     93.913033     /* lunar mean long of node at EPOCH */
  25.  
  26. main()  {
  27.  
  28. double dtor();
  29. double adj360();
  30. double potm();
  31. int i_phase;
  32.  
  33.     long * calloc();
  34.     long *lo = calloc (1, sizeof(long));  /* used by time calls */
  35.     struct tm *pt; /* ptr to time structure */
  36.  
  37.     double days;   /* days since EPOCH */
  38.     double phase;  /* percent of lunar surface illuminated */
  39.     double phase2; /* percent of lunar surface illuminated one day later */
  40.     int i = EPOCH;
  41.  
  42.     time (lo);  /* get system time */
  43.     pt = gmtime(lo);  /* get ptr to gmt time struct */
  44.  
  45.     /* calculate days since EPOCH */
  46.     days = (pt->tm_yday +1) + ((pt->tm_hour + (pt->tm_min / 60.0)
  47.            + (pt->tm_sec / 3600.0)) / 24.0);
  48.     while (i < pt->tm_year + 1900)
  49.        days = days + 365 + ly(i++);
  50.  
  51.     phase = potm(days);
  52.     printf("The Moon is ");
  53.  
  54.     i_phase = phase + 0.5;
  55.  
  56.  
  57.     if (i_phase == 100) {
  58.         printf("Full\n");
  59.     } else {
  60.         if (i_phase == 0) 
  61.             printf("New\n");
  62.         else {
  63.             if (i_phase == 50)  {
  64.                 phase2 = potm(++days);
  65.                 if (phase2 > phase)
  66.                     printf("at the First Quarter\n");
  67.                 else 
  68.                     printf("at the Last Quarter\n");
  69.             } else {
  70.                 if (i_phase > 50) {
  71.                     phase2 = potm(++days);
  72.                     if (phase2 > phase)
  73.                         printf("Waxing ");
  74.                     else 
  75.                         printf("Waning ");
  76.                     printf("Gibbous (%1.0f%% of Full)\n", phase);
  77.                 } else {
  78.                     if (i_phase < 50) {
  79.                         phase2 = potm(++days);
  80.                         if (phase2 > phase)
  81.                             printf("Waxing ");
  82.                         else
  83.                             printf("Waning ");
  84.                         printf("Crescent (%1.0f%% of Full)\n", phase);
  85.                     }
  86.                 }
  87.             }
  88.         }
  89.     }
  90. }
  91.  
  92. double potm(days)
  93. double days;
  94. {
  95. double N;
  96. double Msol;
  97. double Ec;
  98. double LambdaSol;
  99. double l;
  100. double Mm;
  101. double Ev;
  102. double Ac;
  103. double A3;
  104. double Mmprime;
  105. double A4;
  106. double lprime;
  107. double V;
  108. double ldprime;
  109. double D;
  110. double Nm;
  111.  
  112. N = 360 * days / 365.2422;  /* sec 42 #3 */
  113. adj360(&N);
  114.  
  115. Msol = N + EPSILONg - RHOg; /* sec 42 #4 */
  116. adj360(&Msol);
  117.  
  118. Ec = 360 / PI * e * sin(dtor(Msol)); /* sec 42 #5 */
  119.  
  120. LambdaSol = N + Ec + EPSILONg;       /* sec 42 #6 */
  121. adj360(&LambdaSol);
  122.  
  123. l = 13.1763966 * days + lzero;       /* sec 61 #4 */
  124. adj360(&l);
  125.  
  126. Mm = l - (0.1114041 * days) - Pzero; /* sec 61 #5 */
  127. adj360(&Mm);
  128.  
  129. Nm = Nzero - (0.0529539 * days);     /* sec 61 #6 */
  130. adj360(&Nm);
  131.  
  132. Ev = 1.2739 * sin(dtor(2*(l - LambdaSol) - Mm)); /* sec 61 #7 */
  133.  
  134. Ac = 0.1858 * sin(dtor(Msol));       /* sec 61 #8 */
  135. A3 = 0.37 * sin(dtor(Msol));
  136.  
  137. Mmprime = Mm + Ev - Ac - A3;         /* sec 61 #9 */
  138.  
  139. Ec = 6.2886 * sin(dtor(Mmprime));    /* sec 61 #10 */
  140.  
  141. A4 = 0.214 * sin(dtor(2 * Mmprime)); /* sec 61 #11 */
  142.  
  143. lprime = l + Ev + Ec - Ac + A4;      /* sec 61 #12 */
  144.  
  145. V = 0.6583 * sin(dtor(2 * (lprime - LambdaSol))); /* sec 61 #13 */
  146.  
  147. ldprime = lprime + V;                /* sec 61 #14 */
  148.  
  149. D = ldprime - LambdaSol;             /* sec 63 #2 */
  150.  
  151. return (50 * (1 - cos(dtor(D))));    /* sec 63 #3 */
  152. }
  153.  
  154. ly(yr)
  155. int yr;
  156. {
  157. /* returns 1 if leapyear, 0 otherwise */
  158. return (yr % 4 == 0 && yr % 100 != 0 || yr % 400 == 0);
  159. }
  160.  
  161. double dtor(deg)
  162. double deg;
  163. {
  164. /* convert degrees to radians */
  165. return (deg * PI / 180);
  166. }
  167.  
  168. double adj360(deg)
  169. double *deg;
  170. {
  171. /* adjust value so 0 <= deg <= 360 */
  172. do if (*deg < 0)
  173.    *deg += 360;
  174. else if (*deg > 360)
  175.    *deg -= 360;
  176. while (*deg < 0 || *deg > 360);
  177. }
  178.