1100280Sgordon/* 2100280Sgordon * Copyright (c) 1989, 1993 3118638Sfjoe * The Regents of the University of California. All rights reserved. 466830Sobrien * 566830Sobrien * This code is derived from software posted to USENET. 666830Sobrien * 766830Sobrien * Redistribution and use in source and binary forms, with or without 866830Sobrien * modification, are permitted provided that the following conditions 966830Sobrien * are met: 1066830Sobrien * 1. Redistributions of source code must retain the above copyright 1166830Sobrien * notice, this list of conditions and the following disclaimer. 1266830Sobrien * 2. Redistributions in binary form must reproduce the above copyright 1366830Sobrien * notice, this list of conditions and the following disclaimer in the 1466830Sobrien * documentation and/or other materials provided with the distribution. 1566830Sobrien * 3. Neither the name of the University nor the names of its contributors 1666830Sobrien * may be used to endorse or promote products derived from this software 1766830Sobrien * without specific prior written permission. 1866830Sobrien * 1966830Sobrien * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND 2066830Sobrien * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 2166830Sobrien * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 2266830Sobrien * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE 2366830Sobrien * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 2466830Sobrien * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 2566830Sobrien * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 2666830Sobrien * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 2751231Ssheldonh * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 2851231Ssheldonh * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 29100280Sgordon * SUCH DAMAGE. 30126744Spjd */ 31100280Sgordon 32121067Sdougb#if 0 33108191Sdillon#ifndef lint 34108191Sdillonstatic const char copyright[] = 35108191Sdillon"@(#) Copyright (c) 1989, 1993\n\ 36108191Sdillon The Regents of the University of California. All rights reserved.\n"; 37108191Sdillon#endif /* not lint */ 38108191Sdillon 39108191Sdillon#ifndef lint 40108191Sdillonstatic const char sccsid[] = "@(#)pom.c 8.1 (Berkeley) 5/31/93"; 41108191Sdillon#endif /* not lint */ 42108191Sdillon#endif 43108191Sdillon#include <sys/cdefs.h> 44108191Sdillon__FBSDID("$FreeBSD$"); 45108191Sdillon 46108191Sdillon/* 47108191Sdillon * Phase of the Moon. Calculates the current phase of the moon. 48108191Sdillon * Based on routines from `Practical Astronomy with Your Calculator', 49108191Sdillon * by Duffett-Smith. Comments give the section from the book that 50108191Sdillon * particular piece of code was adapted from. 51108191Sdillon * 52126787Sphk * -- Keith E. Brandt VIII 1984 53126787Sphk * 54126787Sphk */ 55126787Sphk 56126787Sphk#include <stdio.h> 57108191Sdillon#include <stdlib.h> 58108191Sdillon#include <math.h> 59108191Sdillon#include <string.h> 60108191Sdillon#include <sysexits.h> 61108191Sdillon#include <time.h> 62126787Sphk#include <unistd.h> 63126787Sphk 64108191Sdillon#ifndef PI 65108191Sdillon#define PI 3.14159265358979323846 66108191Sdillon#endif 67121014Skris#define EPOCH 85 68108191Sdillon#define EPSILONg 279.611371 /* solar ecliptic long at EPOCH */ 69121014Skris#define RHOg 282.680403 /* solar ecliptic long of perigee at EPOCH */ 70108191Sdillon#define ECCEN 0.01671542 /* solar orbit eccentricity */ 71108191Sdillon#define lzero 18.251907 /* lunar mean long at EPOCH */ 72108191Sdillon#define Pzero 192.917585 /* lunar mean long of perigee at EPOCH */ 73108191Sdillon#define Nzero 55.204723 /* lunar mean long of node at EPOCH */ 74108191Sdillon#define isleap(y) ((((y) % 4) == 0 && ((y) % 100) != 0) || ((y) % 400) == 0) 75108191Sdillon 76108191Sdillonstatic void adj360(double *); 77108191Sdillonstatic double dtor(double); 78108191Sdillonstatic double potm(double); 79108191Sdillonstatic void usage(char *progname); 80108191Sdillon 81108191Sdillonint 82108191Sdillonmain(int argc, char **argv) 83108191Sdillon{ 84126868Sbrooks time_t tt; 85126868Sbrooks struct tm GMT, tmd; 86108191Sdillon double days, today, tomorrow; 87100280Sgordon int ch, cnt, pflag = 0; 88126787Sphk char *odate = NULL, *otime = NULL; 89100280Sgordon char *progname = argv[0]; 9043803Sdillon 9143803Sdillon while ((ch = getopt(argc, argv, "d:pt:")) != -1) 9243803Sdillon switch (ch) { 9343803Sdillon case 'd': 9443803Sdillon odate = optarg; 9543803Sdillon break; 9651231Ssheldonh case 'p': 97108191Sdillon pflag = 1; 98108191Sdillon break; 99108191Sdillon case 't': 100108191Sdillon otime = optarg; 101108191Sdillon break; 102108191Sdillon default: 103108191Sdillon usage(progname); 104108191Sdillon } 105108191Sdillon 106108191Sdillon argc -= optind; 10743803Sdillon argv += optind; 10843803Sdillon 109108191Sdillon if (argc) 110108191Sdillon usage(progname); 111108191Sdillon 112108191Sdillon /* Adjust based on users preferences */ 11375931Simp time(&tt); 11475931Simp if (otime != NULL || odate != NULL) { 115108191Sdillon /* Save today in case -d isn't specified */ 116108191Sdillon localtime_r(&tt, &tmd); 117108191Sdillon 118108191Sdillon if (odate != NULL) { 119110942Sjhay tmd.tm_year = strtol(odate, NULL, 10) - 1900; 120121014Skris tmd.tm_mon = strtol(odate + 5, NULL, 10) - 1; 121108191Sdillon tmd.tm_mday = strtol(odate + 8, NULL, 10); 122108191Sdillon /* Use midnight as the middle of the night */ 123108191Sdillon tmd.tm_hour = 0; 124108191Sdillon tmd.tm_min = 0; 125108191Sdillon tmd.tm_sec = 0; 126108191Sdillon } 127108191Sdillon if (otime != NULL) { 128108191Sdillon tmd.tm_hour = strtol(otime, NULL, 10); 129108191Sdillon tmd.tm_min = strtol(otime + 3, NULL, 10); 13043803Sdillon tmd.tm_sec = strtol(otime + 6, NULL, 10); 13143803Sdillon } 13255520Sluigi tt = mktime(&tmd); 13343803Sdillon } 13451231Ssheldonh 13543803Sdillon gmtime_r(&tt, &GMT); 13655520Sluigi days = (GMT.tm_yday + 1) + ((GMT.tm_hour + 13755520Sluigi (GMT.tm_min / 60.0) + (GMT.tm_sec / 3600.0)) / 24.0); 13855520Sluigi for (cnt = EPOCH; cnt < GMT.tm_year; ++cnt) 139126787Sphk days += isleap(1900 + cnt) ? 366 : 365; 140126787Sphk today = potm(days) + .5; 141126787Sphk if (pflag) { 142126787Sphk (void)printf("%1.0f\n", today); 143126787Sphk return (0); 144126787Sphk } 145126787Sphk (void)printf("The Moon is "); 146126787Sphk if ((int)today == 100) 147126787Sphk (void)printf("Full\n"); 148126787Sphk else if (!(int)today) 149126787Sphk (void)printf("New\n"); 150126787Sphk else { 151126787Sphk tomorrow = potm(days + 1); 152126787Sphk if ((int)today == 50) 153126787Sphk (void)printf("%s\n", tomorrow > today ? 154126787Sphk "at the First Quarter" : "at the Last Quarter"); 155126787Sphk else { 156127657Sluigi (void)printf("%s ", tomorrow > today ? 157127657Sluigi "Waxing" : "Waning"); 158127657Sluigi if (today > 50) 159127657Sluigi (void)printf("Gibbous (%1.0f%% of Full)\n", 160127657Sluigi today); 161126787Sphk else if (today < 50) 162126787Sphk (void)printf("Crescent (%1.0f%% of Full)\n", 16343803Sdillon today); 164117087Sbrooks } 165121067Sdougb } 166117087Sbrooks 167117087Sbrooks return 0; 168117087Sbrooks} 169117087Sbrooks 170117087Sbrooks/* 171117087Sbrooks * potm -- 172117087Sbrooks * return phase of the moon 173117087Sbrooks */ 174117087Sbrooksstatic double 175127657Sluigipotm(double days) 176127657Sluigi{ 177127657Sluigi double N, Msol, Ec, LambdaSol, l, Mm, Ev, Ac, A3, Mmprime; 178127657Sluigi double A4, lprime, V, ldprime, D, Nm; 179127657Sluigi 180127657Sluigi N = 360 * days / 365.2422; /* sec 42 #3 */ 181126868Sbrooks adj360(&N); 182126868Sbrooks Msol = N + EPSILONg - RHOg; /* sec 42 #4 */ 183126868Sbrooks adj360(&Msol); 184126868Sbrooks Ec = 360 / PI * ECCEN * sin(dtor(Msol)); /* sec 42 #5 */ 185126868Sbrooks LambdaSol = N + Ec + EPSILONg; /* sec 42 #6 */ 186126868Sbrooks adj360(&LambdaSol); 187126868Sbrooks l = 13.1763966 * days + lzero; /* sec 61 #4 */ 188126868Sbrooks adj360(&l); 189126868Sbrooks Mm = l - (0.1114041 * days) - Pzero; /* sec 61 #5 */ 190127657Sluigi adj360(&Mm); 191126868Sbrooks Nm = Nzero - (0.0529539 * days); /* sec 61 #6 */ 192126868Sbrooks adj360(&Nm); 193108191Sdillon Ev = 1.2739 * sin(dtor(2*(l - LambdaSol) - Mm)); /* sec 61 #7 */ 194121067Sdougb Ac = 0.1858 * sin(dtor(Msol)); /* sec 61 #8 */ 195108191Sdillon A3 = 0.37 * sin(dtor(Msol)); 19643803Sdillon Mmprime = Mm + Ev - Ac - A3; /* sec 61 #9 */ 197108191Sdillon Ec = 6.2886 * sin(dtor(Mmprime)); /* sec 61 #10 */ 198108191Sdillon A4 = 0.214 * sin(dtor(2 * Mmprime)); /* sec 61 #11 */ 199108191Sdillon lprime = l + Ev + Ec - Ac + A4; /* sec 61 #12 */ 200108191Sdillon V = 0.6583 * sin(dtor(2 * (lprime - LambdaSol))); /* sec 61 #13 */ 201108191Sdillon ldprime = lprime + V; /* sec 61 #14 */ 202108191Sdillon D = ldprime - LambdaSol; /* sec 63 #2 */ 203108191Sdillon return(50 * (1 - cos(dtor(D)))); /* sec 63 #3 */ 204108191Sdillon} 205108191Sdillon 206108191Sdillon/* 207108191Sdillon * dtor -- 208117087Sbrooks * convert degrees to radians 209117087Sbrooks */ 210117087Sbrooksstatic double 211117087Sbrooksdtor(double deg) 212108191Sdillon{ 213127657Sluigi 214108191Sdillon return(deg * PI / 180); 215108191Sdillon} 216108191Sdillon 217108191Sdillon/* 218108191Sdillon * adj360 -- 219108191Sdillon * adjust value so 0 <= deg <= 360 220108191Sdillon */ 22143803Sdillonstatic void 222126787Sphkadj360(double *deg) 223126787Sphk{ 224126787Sphk 225126787Sphk for (;;) 226126787Sphk if (*deg < 0) 227126787Sphk *deg += 360; 228127657Sluigi else if (*deg > 360) 229126787Sphk *deg -= 360; 230126787Sphk else 231108191Sdillon break; 232108191Sdillon} 233108191Sdillon 234108191Sdillonstatic void 235117087Sbrooksusage(char *progname) 236117087Sbrooks{ 237117087Sbrooks 238108191Sdillon fprintf(stderr, "Usage: %s [-p] [-d yyyy.mm.dd] [-t hh:mm:ss]\n", 239108191Sdillon progname); 240127657Sluigi exit(EX_USAGE); 24175746Sbsd} 242108191Sdillon