home *** CD-ROM | disk | FTP | other *** search
/ The Fred Fish Collection 1.5 / ffcollection-1-5-1992-11.iso / ff_progs / educ / planets709.lha / moonphase.c < prev    next >
C/C++ Source or Header  |  1992-08-05  |  7KB  |  220 lines
  1. /*******************************************************************************
  2. ** moonphase.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 (Obviously for Aztec C, as)               **
  11. **                                                                            **
  12. **     -- Conversion to SAS/C 5.1b, Clean Up and German Version               **
  13. **        Correct usage of ANSI-C time routines, especially `gmtime'.         **
  14. **        Andreas Scherer VII 1992                                            **
  15. **                                                                            **
  16. *******************************************************************************/
  17.  
  18. #include <stdio.h>
  19. #include <stdlib.h>
  20. #include <time.h>
  21. #include <math.h>
  22.  
  23. #define EPOCH    1983
  24. #define EPSILONg  279.103035   /* solar ecliptic long at EPOCH            */
  25. #define RHOg      282.648015   /* solar ecliptic long of perigee at EPOCH */
  26. #define e           0.01671626 /* solar orbit eccentricity                */
  27. #define lzero     106.306091   /* lunar mean long at EPOCH                */
  28. #define Pzero     111.481526   /* lunar mean long of perigee at EPOCH     */
  29. #define Nzero      93.913033   /* lunar mean long of node at EPOCH        */
  30.  
  31. #define RAD 0.01745329251994329651 /* PI / 180 */
  32.  
  33. char *_TZ = "CET-01CDT";
  34.  
  35. void main()
  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.    long   lo;        /* used in time call */
  41.    int    i = EPOCH;
  42.    int    i_phase;
  43.    struct tm *pt;    /* ptr to time structure */
  44.    double potm();
  45.    int    ly();
  46.  
  47. #ifdef AMIGA
  48.    tzset();
  49. #endif /* AMIGA */
  50.  
  51.    time(&lo);        /* get system time */
  52.    pt = gmtime(&lo); /* get ptr to gmt time struct */
  53.  
  54. /*******************************************************************************
  55. ** Calculate days since EPOCH.                                                **
  56. *******************************************************************************/
  57.    days = (pt->tm_yday + 1.) +
  58.       (((double)pt->tm_hour + (pt->tm_min / 60.) +
  59.        (pt->tm_sec / 3600.)) / 24.);
  60.    while(i < pt->tm_year + 1900)
  61.       days += 365. + (double)ly(i++);
  62.  
  63.    phase = potm(days);
  64.  
  65. #ifdef EUTSCH
  66.    printf("\nDer Mond ist ");
  67. #else
  68.    printf("\nThe Moon is ");
  69. #endif
  70.  
  71.    i_phase = (int)(phase + 0.5);
  72.  
  73.    if(i_phase == 100)
  74.  
  75. #ifdef EUTSCH
  76.       printf("voll\n\n");
  77. #else
  78.       printf("Full\n\n");
  79. #endif
  80.  
  81.    else {
  82.       if(i_phase == 0)
  83.  
  84. #ifdef EUTSCH
  85.          printf("neu\n\n");
  86. #else
  87.          printf("New\n\n");
  88. #endif
  89.  
  90.       else {
  91.          if(i_phase == 50)  {
  92.             phase2 = potm(++days);
  93.             if(phase2 > phase)
  94.  
  95. #ifdef EUTSCH
  96.                printf("im ersten Viertel\n\n");
  97. #else
  98.                printf("at the First Quarter\n\n");
  99. #endif
  100.  
  101.             else 
  102.  
  103. #ifdef EUTSCH
  104.                printf("im letzten Viertel\n\n");
  105. #else
  106.                printf("at the Last Quarter\n\n");
  107. #endif
  108.  
  109.             }
  110.          else {
  111.             if(i_phase > 50) {
  112.                phase2 = potm(++days);
  113.                if(phase2 > phase)
  114.  
  115. #ifdef EUTSCH
  116.                  printf("zunehmend");
  117. #else
  118.                  printf("Waxing");
  119. #endif
  120.  
  121.                else 
  122.  
  123. #ifdef EUTSCH
  124.                   printf("abnehmend");
  125. #else
  126.                   printf("Waning");
  127. #endif
  128.  
  129. #ifdef EUTSCH
  130.                printf(", Halbmond (%1.0f%% voll)\n\n", phase);
  131. #else
  132.                printf(", Gibbous (%1.0f%% of Full)\n\n", phase);
  133. #endif
  134.  
  135.                }
  136.             else {
  137.                if(i_phase < 50) {
  138.                   phase2 = potm(++days);
  139.                   if(phase2 > phase)
  140.  
  141. #ifdef EUTSCH
  142.                      printf("zunehmend");
  143. #else
  144.                      printf("Waxing");
  145. #endif
  146.  
  147.                   else
  148.  
  149. #ifdef EUTSCH
  150.                      printf("abnehmend");
  151. #else
  152.                      printf("Waning");
  153. #endif
  154.  
  155. #ifdef EUTSCH
  156.                   printf(", Halbmond (%1.0f%% voll)\n\n", phase);
  157. #else
  158.                   printf(", Crescent (%1.0f%% of Full)\n\n", phase);
  159. #endif
  160.  
  161.                   }
  162.                }
  163.             }
  164.          }
  165.       }
  166.    }
  167.  
  168. /*******************************************************************************
  169. ** Calculate the Position Of The Moon to a given day fraction.                **
  170. *******************************************************************************/
  171. double potm(days)
  172. double days;
  173.    {
  174.    double N,Msol,LambdaSol,Ev,Mmprime,lprime;
  175.    double dtor(),range();
  176.  
  177.    N = range(360. * days / 365.2422) + EPSILONg;                /* sec 42 #03 */
  178.  
  179.    Msol = sin(dtor(range(N - RHOg)));                           /* sec 42 #04 */
  180.  
  181.    LambdaSol = range(N + 360./PI * e * Msol);                   /* sec 42 #05 */
  182.                                                                 /* sec 42 #06 */
  183.  
  184.    lprime = range(13.1763966 * days + lzero);                   /* sec 61 #04 */
  185.  
  186.    Mmprime = range(lprime - (0.1114041 * days) - Pzero);        /* sec 61 #05 */
  187.  
  188.    Ev = 1.2739 * sin(dtor(2. * (lprime - LambdaSol) - Mmprime)) /* sec 61 #07 */
  189.       - 0.1858 * Msol;                                          /* sec 61 #08 */
  190.  
  191.    Mmprime += Ev - 0.3700 * Msol;                               /* sec 61 #09 */
  192.  
  193.    lprime += Ev                                                 /* sec 61 #10 */
  194.       + 6.2886 * sin(dtor(Mmprime))                             /* sec 61 #11 */
  195.       + 0.2140 * sin(dtor(2. * Mmprime));                       /* sec 61 #12 */
  196.  
  197.    return(50. * (1. - cos(dtor(lprime                           /* sec 61 #13 */
  198.       + 0.6583 * sin(dtor(2. * (lprime - LambdaSol)))           /* sec 61 #14 */
  199.       - LambdaSol))));                                          /* sec 63 #02 */
  200.                                                                 /* sec 63 #03 */
  201.    }
  202.  
  203. /*******************************************************************************
  204. ** Returns 1 if LeapYear, 0 otherwise, according to Gregorian calendar.       **
  205. *******************************************************************************/
  206. int ly(yr)
  207. int yr;
  208.    {
  209.    return(yr % 4 == 0 && yr % 100 != 0 || yr % 400 == 0);
  210.    }
  211.  
  212. /*******************************************************************************
  213. ** Convert Degrees TO Radians                                                 **
  214. *******************************************************************************/
  215. double dtor(deg)
  216. double deg;
  217.    {
  218.    return(deg * RAD);
  219.    }
  220.