* Copyright (c) 1989 The Regents of the University of California.
* This code is derived from software posted to USENET.
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
* 3. All advertising materials mentioning features or use of this software
* must display the following acknowledgement:
* This product includes software developed by the University of
* California, Berkeley and its contributors.
* 4. Neither the name of the University nor the names of its contributors
* may be used to endorse or promote products derived from this software
* without specific prior written permission.
* THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
"@(#) Copyright (c) 1989 The Regents of the University of California.\n\
static char sccsid
[] = "@(#)pom.c 5.3 (Berkeley) 2/28/91";
* Phase of the Moon. Calculates the current phase of the moon.
* Based on routines from `Practical Astronomy with Your Calculator',
* by Duffett-Smith. Comments give the section from the book that
* particular piece of code was adapted from.
* -- Keith E. Brandt VIII 1984
#define EPSILONg 279.611371 /* solar ecliptic long at EPOCH */
#define RHOg 282.680403 /* solar ecliptic long of perigee at EPOCH */
#define ECCEN 0.01671542 /* solar orbit eccentricity */
#define lzero 18.251907 /* lunar mean long at EPOCH */
#define Pzero 192.917585 /* lunar mean long of perigee at EPOCH */
#define Nzero 55.204723 /* lunar mean long of node at EPOCH */
double dtor(), potm(), adj360();
struct tm
*GMT
, *gmtime();
double days
, today
, tomorrow
;
if (gettimeofday(&tp
,&tzp
)) {
(void)fprintf(stderr
, "pom: %s\n", strerror(errno
));
GMT
= gmtime(&tp
.tv_sec
);
days
= (GMT
->tm_yday
+ 1) + ((GMT
->tm_hour
+
(GMT
->tm_min
/ 60.0) + (GMT
->tm_sec
/ 3600.0)) / 24.0);
for (cnt
= EPOCH
; cnt
< GMT
->tm_year
; ++cnt
)
days
+= isleap(cnt
) ? 366 : 365;
(void)printf("The Moon is ");
tomorrow
= potm(days
+ 1);
(void)printf("%s\n", tomorrow
> today
?
"at the First Quarter" : "at the Last Quarter");
(void)printf("%s ", tomorrow
> today
?
(void)printf("Gibbous (%1.0f%% of Full)\n",
(void)printf("Crescent (%1.0f%% of Full)\n",
* return phase of the moon
double N
, Msol
, Ec
, LambdaSol
, l
, Mm
, Ev
, Ac
, A3
, Mmprime
;
double A4
, lprime
, V
, ldprime
, D
, Nm
;
N
= 360 * days
/ 365.2422; /* sec 42 #3 */
Msol
= N
+ EPSILONg
- RHOg
; /* sec 42 #4 */
Ec
= 360 / PI
* ECCEN
* sin(dtor(Msol
)); /* sec 42 #5 */
LambdaSol
= N
+ Ec
+ EPSILONg
; /* sec 42 #6 */
l
= 13.1763966 * days
+ lzero
; /* sec 61 #4 */
Mm
= l
- (0.1114041 * days
) - Pzero
; /* sec 61 #5 */
Nm
= Nzero
- (0.0529539 * days
); /* sec 61 #6 */
Ev
= 1.2739 * sin(dtor(2*(l
- LambdaSol
) - Mm
)); /* sec 61 #7 */
Ac
= 0.1858 * sin(dtor(Msol
)); /* sec 61 #8 */
A3
= 0.37 * sin(dtor(Msol
));
Mmprime
= Mm
+ Ev
- Ac
- A3
; /* sec 61 #9 */
Ec
= 6.2886 * sin(dtor(Mmprime
)); /* sec 61 #10 */
A4
= 0.214 * sin(dtor(2 * Mmprime
)); /* sec 61 #11 */
lprime
= l
+ Ev
+ Ec
- Ac
+ A4
; /* sec 61 #12 */
V
= 0.6583 * sin(dtor(2 * (lprime
- LambdaSol
))); /* sec 61 #13 */
ldprime
= lprime
+ V
; /* sec 61 #14 */
D
= ldprime
- LambdaSol
; /* sec 63 #2 */
return(50 * (1 - cos(dtor(D
)))); /* sec 63 #3 */
* convert degrees to radians
* adjust value so 0 <= deg <= 360