pom.c revision 2490
119370Spst/* 2130803Smarcel * Copyright (c) 1989, 1993 3130803Smarcel * The Regents of the University of California. All rights reserved. 4130803Smarcel * 5130803Smarcel * This code is derived from software posted to USENET. 619370Spst * 719370Spst * Redistribution and use in source and binary forms, with or without 819370Spst * modification, are permitted provided that the following conditions 998944Sobrien * are met: 1019370Spst * 1. Redistributions of source code must retain the above copyright 1198944Sobrien * notice, this list of conditions and the following disclaimer. 1298944Sobrien * 2. Redistributions in binary form must reproduce the above copyright 1398944Sobrien * notice, this list of conditions and the following disclaimer in the 1498944Sobrien * documentation and/or other materials provided with the distribution. 1519370Spst * 3. All advertising materials mentioning features or use of this software 1698944Sobrien * must display the following acknowledgement: 1798944Sobrien * This product includes software developed by the University of 1898944Sobrien * California, Berkeley and its contributors. 1998944Sobrien * 4. Neither the name of the University nor the names of its contributors 2019370Spst * may be used to endorse or promote products derived from this software 2198944Sobrien * without specific prior written permission. 2298944Sobrien * 2398944Sobrien * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND 2419370Spst * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 2519370Spst * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 26130803Smarcel * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE 27130803Smarcel * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 28130803Smarcel * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 2919370Spst * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 30130803Smarcel * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 31130803Smarcel * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 32130803Smarcel * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 33130803Smarcel * SUCH DAMAGE. 34130803Smarcel */ 35130803Smarcel 36130803Smarcel#ifndef lint 37130803Smarcelstatic char copyright[] = 38130803Smarcel"@(#) Copyright (c) 1989, 1993\n\ 39130803Smarcel The Regents of the University of California. All rights reserved.\n"; 40130803Smarcel#endif /* not lint */ 41130803Smarcel 42130803Smarcel#ifndef lint 43130803Smarcelstatic char sccsid[] = "@(#)pom.c 8.1 (Berkeley) 5/31/93"; 44130803Smarcel#endif /* not lint */ 45130803Smarcel 46130803Smarcel/* 47130803Smarcel * Phase of the Moon. Calculates the current phase of the moon. 48130803Smarcel * Based on routines from `Practical Astronomy with Your Calculator', 49130803Smarcel * by Duffett-Smith. Comments give the section from the book that 50130803Smarcel * particular piece of code was adapted from. 51130803Smarcel * 52130803Smarcel * -- Keith E. Brandt VIII 1984 53130803Smarcel * 54130803Smarcel */ 55130803Smarcel 56130803Smarcel#include <sys/time.h> 57130803Smarcel#include <stdio.h> 58130803Smarcel#include <tzfile.h> 59130803Smarcel#include <math.h> 60130803Smarcel 61130803Smarcel#define PI 3.141592654 62130803Smarcel#define EPOCH 85 63130803Smarcel#define EPSILONg 279.611371 /* solar ecliptic long at EPOCH */ 64130803Smarcel#define RHOg 282.680403 /* solar ecliptic long of perigee at EPOCH */ 65130803Smarcel#define ECCEN 0.01671542 /* solar orbit eccentricity */ 66130803Smarcel#define lzero 18.251907 /* lunar mean long at EPOCH */ 67130803Smarcel#define Pzero 192.917585 /* lunar mean long of perigee at EPOCH */ 68130803Smarcel#define Nzero 55.204723 /* lunar mean long of node at EPOCH */ 69130803Smarcel 70130803Smarceldouble dtor(), potm(), adj360(); 71130803Smarcel 72130803Smarcelmain() 73130803Smarcel{ 74130803Smarcel extern int errno; 75130803Smarcel struct timeval tp; 76130803Smarcel struct timezone tzp; 77130803Smarcel struct tm *GMT, *gmtime(); 78130803Smarcel double days, today, tomorrow; 79130803Smarcel int cnt; 80130803Smarcel char *strerror(); 81130803Smarcel 82130803Smarcel if (gettimeofday(&tp,&tzp)) { 83130803Smarcel (void)fprintf(stderr, "pom: %s\n", strerror(errno)); 84130803Smarcel exit(1); 8598944Sobrien } 8698944Sobrien GMT = gmtime(&tp.tv_sec); 8719370Spst days = (GMT->tm_yday + 1) + ((GMT->tm_hour + 8898944Sobrien (GMT->tm_min / 60.0) + (GMT->tm_sec / 3600.0)) / 24.0); 8998944Sobrien for (cnt = EPOCH; cnt < GMT->tm_year; ++cnt) 9098944Sobrien days += isleap(cnt) ? 366 : 365; 9198944Sobrien today = potm(days) + .5; 9298944Sobrien (void)printf("The Moon is "); 9319370Spst if ((int)today == 100) 9498944Sobrien (void)printf("Full\n"); 9598944Sobrien else if (!(int)today) 9619370Spst (void)printf("New\n"); 9798944Sobrien else { 9898944Sobrien tomorrow = potm(days + 1); 9919370Spst if ((int)today == 50) 10098944Sobrien (void)printf("%s\n", tomorrow > today ? 10119370Spst "at the First Quarter" : "at the Last Quarter"); 10219370Spst else { 10319370Spst (void)printf("%s ", tomorrow > today ? 10419370Spst "Waxing" : "Waning"); 10519370Spst if (today > 50) 10619370Spst (void)printf("Gibbous (%1.0f%% of Full)\n", 10719370Spst today); 10819370Spst else if (today < 50) 10998944Sobrien (void)printf("Crescent (%1.0f%% of Full)\n", 11019370Spst today); 11119370Spst } 11219370Spst } 11319370Spst} 11419370Spst 11519370Spst/* 11619370Spst * potm -- 11719370Spst * return phase of the moon 11819370Spst */ 11919370Spstdouble 12019370Spstpotm(days) 12119370Spst double days; 12219370Spst{ 123130803Smarcel double N, Msol, Ec, LambdaSol, l, Mm, Ev, Ac, A3, Mmprime; 124130803Smarcel double A4, lprime, V, ldprime, D, Nm; 12519370Spst 126130803Smarcel N = 360 * days / 365.2422; /* sec 42 #3 */ 127130803Smarcel adj360(&N); 128130803Smarcel Msol = N + EPSILONg - RHOg; /* sec 42 #4 */ 129130803Smarcel adj360(&Msol); 13019370Spst Ec = 360 / PI * ECCEN * sin(dtor(Msol)); /* sec 42 #5 */ 131130803Smarcel LambdaSol = N + Ec + EPSILONg; /* sec 42 #6 */ 132130803Smarcel adj360(&LambdaSol); 13319370Spst l = 13.1763966 * days + lzero; /* sec 61 #4 */ 134130803Smarcel adj360(&l); 135130803Smarcel Mm = l - (0.1114041 * days) - Pzero; /* sec 61 #5 */ 136130803Smarcel adj360(&Mm); 137130803Smarcel Nm = Nzero - (0.0529539 * days); /* sec 61 #6 */ 13819370Spst adj360(&Nm); 139130803Smarcel Ev = 1.2739 * sin(dtor(2*(l - LambdaSol) - Mm)); /* sec 61 #7 */ 140130803Smarcel Ac = 0.1858 * sin(dtor(Msol)); /* sec 61 #8 */ 14119370Spst A3 = 0.37 * sin(dtor(Msol)); 142130803Smarcel Mmprime = Mm + Ev - Ac - A3; /* sec 61 #9 */ 143130803Smarcel Ec = 6.2886 * sin(dtor(Mmprime)); /* sec 61 #10 */ 144130803Smarcel A4 = 0.214 * sin(dtor(2 * Mmprime)); /* sec 61 #11 */ 145130803Smarcel lprime = l + Ev + Ec - Ac + A4; /* sec 61 #12 */ 14619370Spst V = 0.6583 * sin(dtor(2 * (lprime - LambdaSol))); /* sec 61 #13 */ 147130803Smarcel ldprime = lprime + V; /* sec 61 #14 */ 148130803Smarcel D = ldprime - LambdaSol; /* sec 63 #2 */ 14919370Spst return(50 * (1 - cos(dtor(D)))); /* sec 63 #3 */ 150130803Smarcel} 151130803Smarcel 152130803Smarcel/* 153130803Smarcel * dtor -- 15419370Spst * convert degrees to radians 15519370Spst */ 15619370Spstdouble 15719370Spstdtor(deg) 15819370Spst double deg; 15919370Spst{ 16019370Spst return(deg * PI / 180); 16119370Spst} 16219370Spst 16319370Spst/* 16419370Spst * adj360 -- 16519370Spst * adjust value so 0 <= deg <= 360 16619370Spst */ 16719370Spstdouble 16819370Spstadj360(deg) 16919370Spst double *deg; 17019370Spst{ 17119370Spst for (;;) 17219370Spst if (*deg < 0) 17319370Spst *deg += 360; 17419370Spst else if (*deg > 360) 17519370Spst *deg -= 360; 17619370Spst else 17719370Spst break; 17819370Spst} 17919370Spst