home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / CLIPPER / MISC / XBBS7200.ZIP / XBBS7200.TAR / today / moontx.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-05-29  |  3.6 KB  |  165 lines
  1. /****************************************************************************
  2.  moon.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    "moontx.h"
  15.  
  16. struct tm *gmtime();
  17.  
  18. moontxt(buf)
  19. char    buf[];
  20. {
  21.     char *cp=buf;
  22.  
  23. double dtor();
  24. double potm();
  25.  
  26. long *lo = (long *) calloc (1, sizeof(long)); /* used by time calls */
  27. struct tm *pt; /* ptr to time structure */
  28.  
  29. double days;   /* days since EPOCH */
  30. double phase;  /* percent of lunar surface illuminated */
  31. double phase2; /* percent of lunar surface illuminated one day later */
  32. int i = EPOCH;
  33.  
  34. time (lo);  /* get system time */
  35. pt = gmtime(lo);  /* get ptr to gmt time struct */
  36. cfree(lo);
  37.  
  38. /* calculate days since EPOCH */
  39. days = (pt->tm_yday +1.0) + ((pt->tm_hour + (pt->tm_min / 60.0)
  40.        + (pt->tm_sec / 3600.0)) / 24.0);
  41. while (i < pt->tm_year + 1900)
  42.    days = days + 365 + ly(i++);
  43.  
  44. phase = potm(days);
  45. sprintf(cp,"The Moon is ");
  46. cp += strlen(buf);
  47. if ((int)(phase + .5) == 100) {
  48.    sprintf(cp,"Full");
  49.    }
  50. else if ((int)(phase + 0.5) == 0) 
  51.    sprintf(cp,"New");
  52. else if ((int)(phase + 0.5) == 50)  {
  53.    phase2 = potm(++days);
  54.    if (phase2 > phase)
  55.       sprintf(cp,"at the First Quarter");
  56.    else 
  57.       sprintf(cp,"at the Last Quarter");
  58.    }
  59. else if ((int)(phase + 0.5) > 50) {
  60.    phase2 = potm(++days);
  61.    if (phase2 > phase)
  62.       sprintf(cp,"Waxing ");
  63.    else 
  64.       sprintf(cp,"Waning ");
  65.    cp = buf + strlen(buf);
  66.    sprintf(cp,"Gibbous (%1.0f%% of Full)", phase);
  67.    }
  68. else if ((int)(phase + 0.5) < 50) {
  69.    phase2 = potm(++days);
  70.    if (phase2 > phase)
  71.       sprintf(cp,"Waxing ");
  72.    else
  73.       sprintf(cp,"Waning ");
  74.    cp = buf + strlen(buf);
  75.    sprintf(cp,"Crescent (%1.0f%% of Full)", phase);
  76.    }
  77. }
  78.  
  79. double potm(days)
  80. double days;
  81. {
  82. double N;
  83. double Msol;
  84. double Ec;
  85. double LambdaSol;
  86. double l;
  87. double Mm;
  88. double Ev;
  89. double Ac;
  90. double A3;
  91. double Mmprime;
  92. double A4;
  93. double lprime;
  94. double V;
  95. double ldprime;
  96. double D;
  97. double Nm;
  98.  
  99. N = 360 * days / 365.2422;  /* sec 42 #3 */
  100. adj360(&N);
  101.  
  102. Msol = N + EPSILONg - RHOg; /* sec 42 #4 */
  103. adj360(&Msol);
  104.  
  105. Ec = 360 / PI * e * sin(dtor(Msol)); /* sec 42 #5 */
  106.  
  107. LambdaSol = N + Ec + EPSILONg;       /* sec 42 #6 */
  108. adj360(&LambdaSol);
  109.  
  110. l = 13.1763966 * days + lzero;       /* sec 61 #4 */
  111. adj360(&l);
  112.  
  113. Mm = l - (0.1114041 * days) - Pzero; /* sec 61 #5 */
  114. adj360(&Mm);
  115.  
  116. Nm = Nzero - (0.0529539 * days);     /* sec 61 #6 */
  117. adj360(&Nm);
  118.  
  119. Ev = 1.2739 * sin(dtor(2*(l - LambdaSol) - Mm)); /* sec 61 #7 */
  120.  
  121. Ac = 0.1858 * sin(dtor(Msol));       /* sec 61 #8 */
  122. A3 = 0.37 * sin(dtor(Msol));
  123.  
  124. Mmprime = Mm + Ev - Ac - A3;         /* sec 61 #9 */
  125.  
  126. Ec = 6.2886 * sin(dtor(Mmprime));    /* sec 61 #10 */
  127.  
  128. A4 = 0.214 * sin(dtor(2 * Mmprime)); /* sec 61 #11 */
  129.  
  130. lprime = l + Ev + Ec - Ac + A4;      /* sec 61 #12 */
  131.  
  132. V = 0.6583 * sin(dtor(2 * (lprime - LambdaSol))); /* sec 61 #13 */
  133.  
  134. ldprime = lprime + V;                /* sec 61 #14 */
  135.  
  136. D = ldprime - LambdaSol;             /* sec 63 #2 */
  137.  
  138. return (50 * (1 - cos(dtor(D))));    /* sec 63 #3 */
  139. }
  140.  
  141. ly(yr)
  142. int yr;
  143. {
  144. /* returns 1 if leapyear, 0 otherwise */
  145. return (yr % 4 == 0 && yr % 100 != 0 || yr % 400 == 0);
  146. }
  147.  
  148. double dtor(deg)
  149. double deg;
  150. {
  151. /* convert degrees to radians */
  152. return (deg * PI / 180.0);
  153. }
  154.  
  155. adj360(deg)
  156. double *deg;
  157. {
  158. /* adjust value so 0 <= deg <= 360 */
  159. do if (*deg < 0.0)
  160.    *deg += 360.0;
  161. else if (*deg > 360.0)
  162.    *deg -= 360.0;
  163. while (*deg < 0.0 || *deg > 360.0);
  164. }
  165.