|
|
1.1 root 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 are permitted
8: * provided that: (1) source distributions retain this entire copyright
9: * notice and comment, and (2) distributions including binaries display
10: * the following acknowledgement: ``This product includes software
11: * developed by the University of California, Berkeley and its contributors''
12: * in the documentation or other materials provided with the distribution
13: * and in all advertising materials mentioning features or use of this
14: * software. Neither the name of the University nor the names of its
15: * contributors may be used to endorse or promote products derived
16: * from this software without specific prior written permission.
17: * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
18: * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
19: * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
20: */
21:
22: #ifndef lint
23: char copyright[] =
24: "@(#) Copyright (c) 1989 The Regents of the University of California.\n\
25: All rights reserved.\n";
26: #endif /* not lint */
27:
28: #ifndef lint
29: static char sccsid[] = "@(#)pom.c 5.2 (Berkeley) 6/1/90";
30: #endif /* not lint */
31:
32: /*
33: * Phase of the Moon. Calculates the current phase of the moon.
34: * Based on routines from `Practical Astronomy with Your Calculator',
35: * by Duffett-Smith. Comments give the section from the book that
36: * particular piece of code was adapted from.
37: *
38: * -- Keith E. Brandt VIII 1984
39: *
40: */
41:
42: #include <sys/time.h>
43: #include <stdio.h>
44: #include <tzfile.h>
45: #include <math.h>
46:
47: #define PI 3.141592654
48: #define EPOCH 85
49: #define EPSILONg 279.611371 /* solar ecliptic long at EPOCH */
50: #define RHOg 282.680403 /* solar ecliptic long of perigee at EPOCH */
51: #define ECCEN 0.01671542 /* solar orbit eccentricity */
52: #define lzero 18.251907 /* lunar mean long at EPOCH */
53: #define Pzero 192.917585 /* lunar mean long of perigee at EPOCH */
54: #define Nzero 55.204723 /* lunar mean long of node at EPOCH */
55:
56: main()
57: {
58: extern int errno;
59: struct timeval tp;
60: struct timezone tzp;
61: struct tm *GMT, *gmtime();
62: double days, today, tomorrow, dtor(), adj360(), potm();
63: int cnt;
64: char *strerror();
65:
66: if (gettimeofday(&tp,&tzp)) {
67: (void)fprintf(stderr, "pom: %s\n", strerror(errno));
68: exit(1);
69: }
70: GMT = gmtime(&tp.tv_sec);
71: days = (GMT->tm_yday + 1) + ((GMT->tm_hour +
72: (GMT->tm_min / 60.0) + (GMT->tm_sec / 3600.0)) / 24.0);
73: for (cnt = EPOCH; cnt < GMT->tm_year; ++cnt)
74: days += isleap(cnt) ? 366 : 365;
75: today = potm(days) + .5;
76: (void)printf("The Moon is ");
77: if ((int)today == 100)
78: (void)printf("Full\n");
79: else if (!(int)today)
80: (void)printf("New\n");
81: else {
82: tomorrow = potm(days + 1);
83: if ((int)today == 50)
84: (void)printf("%s\n", tomorrow > today ?
85: "at the First Quarter" : "at the Last Quarter");
86: else {
87: (void)printf("%s ", tomorrow > today ?
88: "Waxing" : "Waning");
89: if (today > 50)
90: (void)printf("Gibbous (%1.0f%% of Full)\n",
91: today);
92: else if (today < 50)
93: (void)printf("Crescent (%1.0f%% of Full)\n",
94: today);
95: }
96: }
97: }
98:
99: /*
100: * potm --
101: * return phase of the moon
102: */
103: double
104: potm(days)
105: double days;
106: {
107: double N, Msol, Ec, LambdaSol, l, Mm, Ev, Ac, A3, Mmprime;
108: double A4, lprime, V, ldprime, D, Nm;
109:
110: N = 360 * days / 365.2422; /* sec 42 #3 */
111: adj360(&N);
112: Msol = N + EPSILONg - RHOg; /* sec 42 #4 */
113: adj360(&Msol);
114: Ec = 360 / PI * ECCEN * sin(dtor(Msol)); /* sec 42 #5 */
115: LambdaSol = N + Ec + EPSILONg; /* sec 42 #6 */
116: adj360(&LambdaSol);
117: l = 13.1763966 * days + lzero; /* sec 61 #4 */
118: adj360(&l);
119: Mm = l - (0.1114041 * days) - Pzero; /* sec 61 #5 */
120: adj360(&Mm);
121: Nm = Nzero - (0.0529539 * days); /* sec 61 #6 */
122: adj360(&Nm);
123: Ev = 1.2739 * sin(dtor(2*(l - LambdaSol) - Mm)); /* sec 61 #7 */
124: Ac = 0.1858 * sin(dtor(Msol)); /* sec 61 #8 */
125: A3 = 0.37 * sin(dtor(Msol));
126: Mmprime = Mm + Ev - Ac - A3; /* sec 61 #9 */
127: Ec = 6.2886 * sin(dtor(Mmprime)); /* sec 61 #10 */
128: A4 = 0.214 * sin(dtor(2 * Mmprime)); /* sec 61 #11 */
129: lprime = l + Ev + Ec - Ac + A4; /* sec 61 #12 */
130: V = 0.6583 * sin(dtor(2 * (lprime - LambdaSol))); /* sec 61 #13 */
131: ldprime = lprime + V; /* sec 61 #14 */
132: D = ldprime - LambdaSol; /* sec 63 #2 */
133: return(50 * (1 - cos(dtor(D)))); /* sec 63 #3 */
134: }
135:
136: /*
137: * dtor --
138: * convert degrees to radians
139: */
140: double
141: dtor(deg)
142: double deg;
143: {
144: return(deg * PI / 180);
145: }
146:
147: /*
148: * adj360 --
149: * adjust value so 0 <= deg <= 360
150: */
151: double
152: adj360(deg)
153: double *deg;
154: {
155: for (;;)
156: if (*deg < 0)
157: *deg += 360;
158: else if (*deg > 360)
159: *deg -= 360;
160: else
161: break;
162: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.