source: buchla-68k/orig/DATE/MOONTX.C@ 0170798

Last change on this file since 0170798 was 3ae31e9, checked in by Thomas Lopatic <thomas@…>, 7 years ago

Imported original source code.

  • Property mode set to 100755
File size: 4.6 KB
Line 
1/*
2 =============================================================================
3 moontx.c -- Calculate the current phase of the moon.
4 Version 2 -- 1989-01-10 -- D.N. Lynx Crowe
5
6 Original code by Keith E. Brandt VIII 1984
7
8 Based on routines from "Practical Astronomy with Your Calculator",
9 by Peter Duffett-Smith.
10
11 Comments give the section of the book that particular piece
12 of code was adapted from.
13 =============================================================================
14*/
15
16#define MAINPROG 1
17
18#include "stdio.h"
19
20#ifdef VAX
21#include "sys/time.h"
22#else VAX
23#include "time.h"
24#endif VAX
25
26#include "math.h"
27
28#define EPOCH 1985
29#define EPSILONg 279.611371 /* solar ecliptic long at EPOCH */
30#define RHOg 282.680403 /* solar ecliptic long of perigee at EPOCH */
31#define e 0.01671542 /* solar orbit eccentricity */
32#define lzero 18.251907 /* lunar mean long at EPOCH */
33#define Pzero 192.917585 /* lunar mean long of perigee at EPOCH */
34#define Nzero 55.204723 /* lunar mean long of node at EPOCH */
35
36#ifndef PI
37#define PI 3.141592654
38#endif
39
40int
41ly(yr)
42int yr;
43{
44 /* returns 1 if leapyear, 0 otherwise */
45
46 return (yr % 4 == 0 && yr % 100 != 0 || yr % 400 == 0);
47}
48
49double
50dtor(deg)
51double deg;
52{
53 /* convert degrees to radians */
54
55 return (deg * PI / 180);
56}
57
58double
59adj360(deg)
60double *deg;
61{
62 /* adjust value so 0 <= deg <= 360 */
63
64 do {
65
66 if (*deg < 0)
67 *deg += 360;
68 else if (*deg > 360)
69 *deg -= 360;
70
71 } while (*deg < 0 || *deg > 360);
72}
73
74double
75potm(days)
76double days;
77{
78 double N;
79 double Msol;
80 double Ec;
81 double LambdaSol;
82 double l;
83 double Mm;
84 double Ev;
85 double Ac;
86 double A3;
87 double Mmprime;
88 double A4;
89 double lprime;
90 double V;
91 double ldprime;
92 double D;
93 double Nm;
94
95 N = 360 * days / 365.2422; /* sec 42 #3 */
96 adj360(&N);
97
98 Msol = N + EPSILONg - RHOg; /* sec 42 #4 */
99 adj360(&Msol);
100
101 Ec = 360 / PI * e * sin(dtor(Msol)); /* sec 42 #5 */
102
103 LambdaSol = N + Ec + EPSILONg; /* sec 42 #6 */
104 adj360(&LambdaSol);
105
106 l = 13.1763966 * days + lzero; /* sec 61 #4 */
107 adj360(&l);
108
109 Mm = l - (0.1114041 * days) - Pzero; /* sec 61 #5 */
110 adj360(&Mm);
111
112 Nm = Nzero - (0.0529539 * days); /* sec 61 #6 */
113 adj360(&Nm);
114
115 Ev = 1.2739 * sin(dtor(2*(l - LambdaSol) - Mm)); /* sec 61 #7 */
116
117 Ac = 0.1858 * sin(dtor(Msol)); /* sec 61 #8 */
118 A3 = 0.37 * sin(dtor(Msol));
119
120 Mmprime = Mm + Ev - Ac - A3; /* sec 61 #9 */
121
122 Ec = 6.2886 * sin(dtor(Mmprime)); /* sec 61 #10 */
123
124 A4 = 0.214 * sin(dtor(2 * Mmprime)); /* sec 61 #11 */
125
126 lprime = l + Ev + Ec - Ac + A4; /* sec 61 #12 */
127
128 V = 0.6583 * sin(dtor(2 * (lprime - LambdaSol))); /* sec 61 #13 */
129
130 ldprime = lprime + V; /* sec 61 #14 */
131
132 D = ldprime - LambdaSol; /* sec 63 #2 */
133
134 return(50 * (1 - cos(dtor(D)))); /* sec 63 #3 */
135}
136
137moontxt(buf)
138char buf[];
139{
140 char *cp = buf;
141 long *lo; /* used by time calls */
142 long curtime;
143 double days; /* days since EPOCH */
144 double phase; /* percent of lunar surface illuminated */
145 double phase2; /* percent of lunar surface illuminated one day later */
146 int i = EPOCH;
147 struct tm *pt; /* ptr to time structure */
148
149 lo = &curtime;
150 time(lo); /* get system time */
151 pt = gmtime(lo); /* get ptr to gmt time struct */
152
153 /* calculate days since EPOCH */
154
155 days = (pt->tm_yday +1) + ((pt->tm_hour + (pt->tm_min / 60.0)
156 + (pt->tm_sec / 3600.0)) / 24.0);
157
158 while (i < pt->tm_year + 1900)
159 days = days + 365 + ly(i++);
160
161 phase = potm(days);
162 sprintf(cp,"The Moon is ");
163 cp += strlen(buf);
164
165 if ((int)(phase + .5) == 100)
166 sprintf(cp,"Full");
167 else if ((int)(phase + 0.5) == 0)
168 sprintf(cp,"New");
169 else if ((int)(phase + 0.5) == 50) {
170
171 phase2 = potm(++days);
172
173 if (phase2 > phase)
174 sprintf(cp,"at the First Quarter");
175 else
176 sprintf(cp,"at the Last Quarter");
177
178 } else if ((int)(phase + 0.5) > 50) {
179
180 phase2 = potm(++days);
181
182 if (phase2 > phase)
183 sprintf(cp,"Waxing ");
184 else
185 sprintf(cp,"Waning ");
186
187 cp = buf + strlen(buf);
188 sprintf(cp,"Gibbous (%1.0f%% of Full)", phase);
189
190 } else if ((int)(phase + 0.5) < 50) {
191
192 phase2 = potm(++days);
193
194 if (phase2 > phase)
195 sprintf(cp,"Waxing ");
196 else
197 sprintf(cp,"Waning ");
198
199 cp = buf + strlen(buf);
200 sprintf(cp,"Crescent (%1.0f%% of Full)", phase);
201 }
202}
203
204#if MAINPROG
205
206char mstring[128];
207
208main()
209{
210 long curtime;
211
212 moontxt(mstring);
213 time(&curtime);
214
215 printf("%s", ctime(&curtime));
216 printf("%s\n", mstring);
217}
218
219#endif
Note: See TracBrowser for help on using the repository browser.