12490Sjkh/* 22490Sjkh * Copyright (c) 1989, 1993 32490Sjkh * The Regents of the University of California. All rights reserved. 42490Sjkh * 52490Sjkh * This code is derived from software posted to USENET. 62490Sjkh * 72490Sjkh * Redistribution and use in source and binary forms, with or without 82490Sjkh * modification, are permitted provided that the following conditions 92490Sjkh * are met: 102490Sjkh * 1. Redistributions of source code must retain the above copyright 112490Sjkh * notice, this list of conditions and the following disclaimer. 122490Sjkh * 2. Redistributions in binary form must reproduce the above copyright 132490Sjkh * notice, this list of conditions and the following disclaimer in the 142490Sjkh * documentation and/or other materials provided with the distribution. 15203932Simp * 3. Neither the name of the University nor the names of its contributors 162490Sjkh * may be used to endorse or promote products derived from this software 172490Sjkh * without specific prior written permission. 182490Sjkh * 192490Sjkh * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND 202490Sjkh * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 212490Sjkh * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 222490Sjkh * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE 232490Sjkh * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 242490Sjkh * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 252490Sjkh * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 262490Sjkh * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 272490Sjkh * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 282490Sjkh * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 292490Sjkh * SUCH DAMAGE. 302490Sjkh */ 312490Sjkh 32114725Sobrien#if 0 332490Sjkh#ifndef lint 3415949Sachestatic const char copyright[] = 352490Sjkh"@(#) Copyright (c) 1989, 1993\n\ 362490Sjkh The Regents of the University of California. All rights reserved.\n"; 372490Sjkh#endif /* not lint */ 382490Sjkh 392490Sjkh#ifndef lint 4015949Sachestatic const char sccsid[] = "@(#)pom.c 8.1 (Berkeley) 5/31/93"; 41114725Sobrien#endif /* not lint */ 4253920Sbillf#endif 43114725Sobrien#include <sys/cdefs.h> 44114725Sobrien__FBSDID("$FreeBSD$"); 452490Sjkh 462490Sjkh/* 472490Sjkh * Phase of the Moon. Calculates the current phase of the moon. 482490Sjkh * Based on routines from `Practical Astronomy with Your Calculator', 492490Sjkh * by Duffett-Smith. Comments give the section from the book that 502490Sjkh * particular piece of code was adapted from. 512490Sjkh * 522490Sjkh * -- Keith E. Brandt VIII 1984 532490Sjkh * 542490Sjkh */ 552490Sjkh 562490Sjkh#include <stdio.h> 57201613Sedwin#include <stdlib.h> 582490Sjkh#include <math.h> 59201613Sedwin#include <string.h> 60201613Sedwin#include <sysexits.h> 61201613Sedwin#include <time.h> 62201613Sedwin#include <unistd.h> 632490Sjkh 6445801Ssteve#ifndef PI 6545801Ssteve#define PI 3.14159265358979323846 6645801Ssteve#endif 672490Sjkh#define EPOCH 85 682490Sjkh#define EPSILONg 279.611371 /* solar ecliptic long at EPOCH */ 692490Sjkh#define RHOg 282.680403 /* solar ecliptic long of perigee at EPOCH */ 702490Sjkh#define ECCEN 0.01671542 /* solar orbit eccentricity */ 712490Sjkh#define lzero 18.251907 /* lunar mean long at EPOCH */ 722490Sjkh#define Pzero 192.917585 /* lunar mean long of perigee at EPOCH */ 732490Sjkh#define Nzero 55.204723 /* lunar mean long of node at EPOCH */ 7415949Sache#define isleap(y) ((((y) % 4) == 0 && ((y) % 100) != 0) || ((y) % 400) == 0) 752490Sjkh 7690828Simpstatic void adj360(double *); 7790828Simpstatic double dtor(double); 7890828Simpstatic double potm(double); 79201613Sedwinstatic void usage(char *progname); 802490Sjkh 8117201Sjoergint 82201613Sedwinmain(int argc, char **argv) 832490Sjkh{ 8415949Sache time_t tt; 85201613Sedwin struct tm GMT, tmd; 862490Sjkh double days, today, tomorrow; 87210089Semaste int ch, cnt, pflag = 0; 88201627Sdelphij char *odate = NULL, *otime = NULL; 89230644Smaxim char *progname = argv[0]; 902490Sjkh 91210089Semaste while ((ch = getopt(argc, argv, "d:pt:")) != -1) 92201613Sedwin switch (ch) { 93201613Sedwin case 'd': 94201613Sedwin odate = optarg; 95201613Sedwin break; 96210089Semaste case 'p': 97210089Semaste pflag = 1; 98210089Semaste break; 99201613Sedwin case 't': 100201613Sedwin otime = optarg; 101201613Sedwin break; 102201613Sedwin default: 103230644Smaxim usage(progname); 104201613Sedwin } 105201613Sedwin 106201613Sedwin argc -= optind; 107201613Sedwin argv += optind; 108201613Sedwin 109201613Sedwin if (argc) 110230644Smaxim usage(progname); 111201613Sedwin 112201613Sedwin /* Adjust based on users preferences */ 113201613Sedwin time(&tt); 114201613Sedwin if (otime != NULL || odate != NULL) { 115201613Sedwin /* Save today in case -d isn't specified */ 116201613Sedwin localtime_r(&tt, &tmd); 117201613Sedwin 118201613Sedwin if (odate != NULL) { 119201613Sedwin tmd.tm_year = strtol(odate, NULL, 10) - 1900; 120201613Sedwin tmd.tm_mon = strtol(odate + 5, NULL, 10) - 1; 121201613Sedwin tmd.tm_mday = strtol(odate + 8, NULL, 10); 122201613Sedwin /* Use midnight as the middle of the night */ 123201613Sedwin tmd.tm_hour = 0; 124201613Sedwin tmd.tm_min = 0; 125201613Sedwin tmd.tm_sec = 0; 126201613Sedwin } 127201613Sedwin if (otime != NULL) { 128201613Sedwin tmd.tm_hour = strtol(otime, NULL, 10); 129201613Sedwin tmd.tm_min = strtol(otime + 3, NULL, 10); 130201613Sedwin tmd.tm_sec = strtol(otime + 6, NULL, 10); 131201613Sedwin } 132201613Sedwin tt = mktime(&tmd); 133201613Sedwin } 134201613Sedwin 135201613Sedwin gmtime_r(&tt, &GMT); 136201613Sedwin days = (GMT.tm_yday + 1) + ((GMT.tm_hour + 137201613Sedwin (GMT.tm_min / 60.0) + (GMT.tm_sec / 3600.0)) / 24.0); 138201613Sedwin for (cnt = EPOCH; cnt < GMT.tm_year; ++cnt) 13937868Simp days += isleap(1900 + cnt) ? 366 : 365; 1402490Sjkh today = potm(days) + .5; 141210089Semaste if (pflag) { 142210089Semaste (void)printf("%1.0f\n", today); 143210089Semaste return (0); 144210089Semaste } 1452490Sjkh (void)printf("The Moon is "); 1462490Sjkh if ((int)today == 100) 1472490Sjkh (void)printf("Full\n"); 1482490Sjkh else if (!(int)today) 1492490Sjkh (void)printf("New\n"); 1502490Sjkh else { 1512490Sjkh tomorrow = potm(days + 1); 1522490Sjkh if ((int)today == 50) 1532490Sjkh (void)printf("%s\n", tomorrow > today ? 1542490Sjkh "at the First Quarter" : "at the Last Quarter"); 1552490Sjkh else { 1562490Sjkh (void)printf("%s ", tomorrow > today ? 1572490Sjkh "Waxing" : "Waning"); 1582490Sjkh if (today > 50) 1592490Sjkh (void)printf("Gibbous (%1.0f%% of Full)\n", 1602490Sjkh today); 1612490Sjkh else if (today < 50) 1622490Sjkh (void)printf("Crescent (%1.0f%% of Full)\n", 1632490Sjkh today); 1642490Sjkh } 1652490Sjkh } 16617201Sjoerg 16717201Sjoerg return 0; 1682490Sjkh} 1692490Sjkh 1702490Sjkh/* 1712490Sjkh * potm -- 1722490Sjkh * return phase of the moon 1732490Sjkh */ 17417201Sjoergstatic double 175145782Sstefanfpotm(double days) 1762490Sjkh{ 1772490Sjkh double N, Msol, Ec, LambdaSol, l, Mm, Ev, Ac, A3, Mmprime; 1782490Sjkh double A4, lprime, V, ldprime, D, Nm; 1792490Sjkh 1802490Sjkh N = 360 * days / 365.2422; /* sec 42 #3 */ 1812490Sjkh adj360(&N); 1822490Sjkh Msol = N + EPSILONg - RHOg; /* sec 42 #4 */ 1832490Sjkh adj360(&Msol); 1842490Sjkh Ec = 360 / PI * ECCEN * sin(dtor(Msol)); /* sec 42 #5 */ 1852490Sjkh LambdaSol = N + Ec + EPSILONg; /* sec 42 #6 */ 1862490Sjkh adj360(&LambdaSol); 1872490Sjkh l = 13.1763966 * days + lzero; /* sec 61 #4 */ 1882490Sjkh adj360(&l); 1892490Sjkh Mm = l - (0.1114041 * days) - Pzero; /* sec 61 #5 */ 1902490Sjkh adj360(&Mm); 1912490Sjkh Nm = Nzero - (0.0529539 * days); /* sec 61 #6 */ 1922490Sjkh adj360(&Nm); 1932490Sjkh Ev = 1.2739 * sin(dtor(2*(l - LambdaSol) - Mm)); /* sec 61 #7 */ 1942490Sjkh Ac = 0.1858 * sin(dtor(Msol)); /* sec 61 #8 */ 1952490Sjkh A3 = 0.37 * sin(dtor(Msol)); 1962490Sjkh Mmprime = Mm + Ev - Ac - A3; /* sec 61 #9 */ 1972490Sjkh Ec = 6.2886 * sin(dtor(Mmprime)); /* sec 61 #10 */ 1982490Sjkh A4 = 0.214 * sin(dtor(2 * Mmprime)); /* sec 61 #11 */ 1992490Sjkh lprime = l + Ev + Ec - Ac + A4; /* sec 61 #12 */ 2002490Sjkh V = 0.6583 * sin(dtor(2 * (lprime - LambdaSol))); /* sec 61 #13 */ 2012490Sjkh ldprime = lprime + V; /* sec 61 #14 */ 2022490Sjkh D = ldprime - LambdaSol; /* sec 63 #2 */ 2032490Sjkh return(50 * (1 - cos(dtor(D)))); /* sec 63 #3 */ 2042490Sjkh} 2052490Sjkh 2062490Sjkh/* 2072490Sjkh * dtor -- 2082490Sjkh * convert degrees to radians 2092490Sjkh */ 21017201Sjoergstatic double 211145782Sstefanfdtor(double deg) 2122490Sjkh{ 213201613Sedwin 2142490Sjkh return(deg * PI / 180); 2152490Sjkh} 2162490Sjkh 2172490Sjkh/* 2182490Sjkh * adj360 -- 2192490Sjkh * adjust value so 0 <= deg <= 360 2202490Sjkh */ 22117201Sjoergstatic void 222145782Sstefanfadj360(double *deg) 2232490Sjkh{ 224201613Sedwin 2252490Sjkh for (;;) 2262490Sjkh if (*deg < 0) 2272490Sjkh *deg += 360; 2282490Sjkh else if (*deg > 360) 2292490Sjkh *deg -= 360; 2302490Sjkh else 2312490Sjkh break; 2322490Sjkh} 233201613Sedwin 234201613Sedwinstatic void 235201613Sedwinusage(char *progname) 236201613Sedwin{ 237201613Sedwin 238210089Semaste fprintf(stderr, "Usage: %s [-p] [-d yyyy.mm.dd] [-t hh:mm:ss]\n", 239210089Semaste progname); 240201613Sedwin exit(EX_USAGE); 241201613Sedwin} 242