1331722Seadler/* 2205821Sedwin * Copyright (c) 1989, 1993 3205821Sedwin * The Regents of the University of California. All rights reserved. 4205821Sedwin * 5205821Sedwin * This code is derived from software posted to USENET. 6205821Sedwin * 7205821Sedwin * Redistribution and use in source and binary forms, with or without 8205821Sedwin * modification, are permitted provided that the following conditions 9205821Sedwin * are met: 10205821Sedwin * 1. Redistributions of source code must retain the above copyright 11205821Sedwin * notice, this list of conditions and the following disclaimer. 12205821Sedwin * 2. Redistributions in binary form must reproduce the above copyright 13205821Sedwin * notice, this list of conditions and the following disclaimer in the 14205821Sedwin * documentation and/or other materials provided with the distribution. 15205821Sedwin * 4. Neither the name of the University nor the names of its contributors 16205821Sedwin * may be used to endorse or promote products derived from this software 17205821Sedwin * without specific prior written permission. 18205821Sedwin * 19205821Sedwin * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND 20205821Sedwin * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 21205821Sedwin * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 22205821Sedwin * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE 23205821Sedwin * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 24205821Sedwin * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 25205821Sedwin * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 26205821Sedwin * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 27205821Sedwin * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 28205821Sedwin * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 29205821Sedwin * SUCH DAMAGE. 30205821Sedwin */ 31205821Sedwin 32205821Sedwin#if 0 33205821Sedwin#ifndef lint 34205821Sedwinstatic const char copyright[] = 35205821Sedwin"@(#) Copyright (c) 1989, 1993\n\ 36205821Sedwin The Regents of the University of California. All rights reserved.\n"; 37205821Sedwin#endif /* not lint */ 38205821Sedwin 39205821Sedwin#ifndef lint 40205821Sedwinstatic const char sccsid[] = "@(#)pom.c 8.1 (Berkeley) 5/31/93"; 41205821Sedwin#endif /* not lint */ 42205821Sedwin#endif 43205821Sedwin#include <sys/cdefs.h> 44205821Sedwin__FBSDID("$FreeBSD$"); 45205821Sedwin 46205821Sedwin/* 47205821Sedwin * Phase of the Moon. Calculates the current phase of the moon. 48205821Sedwin * Based on routines from `Practical Astronomy with Your Calculator', 49205821Sedwin * by Duffett-Smith. Comments give the section from the book that 50205821Sedwin * particular piece of code was adapted from. 51205821Sedwin * 52205821Sedwin * -- Keith E. Brandt VIII 1984 53205821Sedwin * 54205821Sedwin */ 55205821Sedwin 56205821Sedwin#include <stdio.h> 57205821Sedwin#include <stdlib.h> 58205821Sedwin#include <math.h> 59205821Sedwin#include <string.h> 60205821Sedwin#include <sysexits.h> 61205821Sedwin#include <time.h> 62205821Sedwin#include <unistd.h> 63205821Sedwin 64205821Sedwin#include "calendar.h" 65205821Sedwin 66205821Sedwin#ifndef PI 67205821Sedwin#define PI 3.14159265358979323846 68205821Sedwin#endif 69205821Sedwin#define EPOCH 85 70205821Sedwin#define EPSILONg 279.611371 /* solar ecliptic long at EPOCH */ 71205821Sedwin#define RHOg 282.680403 /* solar ecliptic long of perigee at EPOCH */ 72205821Sedwin#define ECCEN 0.01671542 /* solar orbit eccentricity */ 73205821Sedwin#define lzero 18.251907 /* lunar mean long at EPOCH */ 74205821Sedwin#define Pzero 192.917585 /* lunar mean long of perigee at EPOCH */ 75205821Sedwin#define Nzero 55.204723 /* lunar mean long of node at EPOCH */ 76205821Sedwin#define isleap(y) ((((y) % 4) == 0 && ((y) % 100) != 0) || ((y) % 400) == 0) 77205821Sedwin 78205821Sedwinstatic void adj360(double *); 79205821Sedwinstatic double dtor(double); 80205821Sedwinstatic double potm(double onday); 81205821Sedwinstatic double potm_minute(double onday, int olddir); 82205821Sedwin 83205821Sedwinvoid 84205821Sedwinpom(int year, double utcoffset, int *fms, int *nms) 85205821Sedwin{ 86205821Sedwin double ffms[MAXMOONS]; 87205821Sedwin double fnms[MAXMOONS]; 88205821Sedwin int i, j; 89205821Sedwin 90205821Sedwin fpom(year, utcoffset, ffms, fnms); 91205821Sedwin 92205821Sedwin j = 0; 93205821Sedwin for (i = 0; ffms[i] != 0; i++) 94205821Sedwin fms[j++] = round(ffms[i]); 95205821Sedwin fms[i] = -1; 96205821Sedwin for (i = 0; fnms[i] != 0; i++) 97205821Sedwin nms[i] = round(fnms[i]); 98205821Sedwin nms[i] = -1; 99205821Sedwin} 100205821Sedwin 101205821Sedwinvoid 102205821Sedwinfpom(int year, double utcoffset, double *ffms, double *fnms) 103205821Sedwin{ 104205821Sedwin time_t tt; 105205821Sedwin struct tm GMT, tmd_today, tmd_tomorrow; 106205821Sedwin double days_today, days_tomorrow, today, tomorrow; 107205821Sedwin int cnt, d; 108205821Sedwin int yeardays; 109205821Sedwin int olddir, newdir; 110205821Sedwin double *pfnms, *pffms, t; 111205821Sedwin 112205821Sedwin pfnms = fnms; 113205821Sedwin pffms = ffms; 114205821Sedwin 115205821Sedwin /* 116205821Sedwin * We take the phase of the moon one second before and one second 117205821Sedwin * after midnight. 118205821Sedwin */ 119205821Sedwin memset(&tmd_today, 0, sizeof(tmd_today)); 120205821Sedwin tmd_today.tm_year = year - 1900; 121205821Sedwin tmd_today.tm_mon = 0; 122205821Sedwin tmd_today.tm_mday = -1; /* 31 December */ 123205821Sedwin tmd_today.tm_hour = 23; 124205821Sedwin tmd_today.tm_min = 59; 125205821Sedwin tmd_today.tm_sec = 59; 126205821Sedwin memset(&tmd_tomorrow, 0, sizeof(tmd_tomorrow)); 127205821Sedwin tmd_tomorrow.tm_year = year - 1900; 128205821Sedwin tmd_tomorrow.tm_mon = 0; 129205821Sedwin tmd_tomorrow.tm_mday = 0; /* 01 January */ 130205821Sedwin tmd_tomorrow.tm_hour = 0; 131205821Sedwin tmd_tomorrow.tm_min = 0; 132205821Sedwin tmd_tomorrow.tm_sec = 1; 133205821Sedwin 134205821Sedwin tt = mktime(&tmd_today); 135205821Sedwin gmtime_r(&tt, &GMT); 136205821Sedwin yeardays = 0; 137205821Sedwin for (cnt = EPOCH; cnt < GMT.tm_year; ++cnt) 138205821Sedwin yeardays += isleap(1900 + cnt) ? DAYSPERLEAPYEAR : DAYSPERYEAR; 139205821Sedwin days_today = (GMT.tm_yday + 1) + ((GMT.tm_hour + 140205821Sedwin (GMT.tm_min / FSECSPERMINUTE) + (GMT.tm_sec / FSECSPERHOUR)) / 141205821Sedwin FHOURSPERDAY); 142205821Sedwin days_today += yeardays; 143205821Sedwin 144205821Sedwin tt = mktime(&tmd_tomorrow); 145205821Sedwin gmtime_r(&tt, &GMT); 146205821Sedwin yeardays = 0; 147205821Sedwin for (cnt = EPOCH; cnt < GMT.tm_year; ++cnt) 148205821Sedwin yeardays += isleap(1900 + cnt) ? DAYSPERLEAPYEAR : DAYSPERYEAR; 149205821Sedwin days_tomorrow = (GMT.tm_yday + 1) + ((GMT.tm_hour + 150205821Sedwin (GMT.tm_min / FSECSPERMINUTE) + (GMT.tm_sec / FSECSPERHOUR)) / 151205821Sedwin FHOURSPERDAY); 152205821Sedwin days_tomorrow += yeardays; 153205821Sedwin 154205821Sedwin today = potm(days_today); /* 30 December 23:59:59 */ 155205821Sedwin tomorrow = potm(days_tomorrow); /* 31 December 00:00:01 */ 156205821Sedwin olddir = today > tomorrow ? -1 : +1; 157205821Sedwin 158223923Sdelphij yeardays = 1 + (isleap(year) ? DAYSPERLEAPYEAR : DAYSPERYEAR); /* reuse */ 159205821Sedwin for (d = 0; d <= yeardays; d++) { 160205821Sedwin today = potm(days_today); 161205821Sedwin tomorrow = potm(days_tomorrow); 162205821Sedwin newdir = today > tomorrow ? -1 : +1; 163205821Sedwin if (olddir != newdir) { 164205821Sedwin t = potm_minute(days_today - 1, olddir) + 165205821Sedwin utcoffset / FHOURSPERDAY; 166205821Sedwin if (olddir == -1 && newdir == +1) { 167205821Sedwin *pfnms = d - 1 + t; 168205821Sedwin pfnms++; 169205821Sedwin } else if (olddir == +1 && newdir == -1) { 170205821Sedwin *pffms = d - 1 + t; 171205821Sedwin pffms++; 172205821Sedwin } 173205821Sedwin } 174205821Sedwin olddir = newdir; 175205821Sedwin days_today++; 176205821Sedwin days_tomorrow++; 177205821Sedwin } 178205821Sedwin *pffms = -1; 179205821Sedwin *pfnms = -1; 180205821Sedwin} 181205821Sedwin 182205821Sedwinstatic double 183205821Sedwinpotm_minute(double onday, int olddir) { 184205821Sedwin double period = FSECSPERDAY / 2.0; 185205821Sedwin double p1, p2; 186205821Sedwin double before, after; 187205821Sedwin int newdir; 188205821Sedwin 189205821Sedwin// printf("---> days:%g olddir:%d\n", days, olddir); 190205821Sedwin 191205821Sedwin p1 = onday + (period / SECSPERDAY); 192205821Sedwin period /= 2; 193205821Sedwin 194205821Sedwin while (period > 30) { /* half a minute */ 195205821Sedwin// printf("period:%g - p1:%g - ", period, p1); 196205821Sedwin p2 = p1 + (2.0 / SECSPERDAY); 197205821Sedwin before = potm(p1); 198205821Sedwin after = potm(p2); 199205821Sedwin// printf("before:%10.10g - after:%10.10g\n", before, after); 200205821Sedwin newdir = before < after ? -1 : +1; 201205821Sedwin if (olddir != newdir) 202205821Sedwin p1 += (period / SECSPERDAY); 203205821Sedwin else 204205821Sedwin p1 -= (period / SECSPERDAY); 205205821Sedwin period /= 2; 206205821Sedwin// printf("newdir:%d - p1:%10.10f - period:%g\n", 207205821Sedwin// newdir, p1, period); 208205821Sedwin } 209205821Sedwin p1 -= floor(p1); 210205821Sedwin //exit(0); 211205821Sedwin return (p1); 212205821Sedwin} 213205821Sedwin 214205821Sedwin/* 215205821Sedwin * potm -- 216205821Sedwin * return phase of the moon, as a percentage [0 ... 100] 217205821Sedwin */ 218205821Sedwinstatic double 219205821Sedwinpotm(double onday) 220205821Sedwin{ 221205821Sedwin double N, Msol, Ec, LambdaSol, l, Mm, Ev, Ac, A3, Mmprime; 222205821Sedwin double A4, lprime, V, ldprime, D, Nm; 223205821Sedwin 224205821Sedwin N = 360 * onday / 365.2422; /* sec 42 #3 */ 225205821Sedwin adj360(&N); 226205821Sedwin Msol = N + EPSILONg - RHOg; /* sec 42 #4 */ 227205821Sedwin adj360(&Msol); 228205821Sedwin Ec = 360 / PI * ECCEN * sin(dtor(Msol)); /* sec 42 #5 */ 229205821Sedwin LambdaSol = N + Ec + EPSILONg; /* sec 42 #6 */ 230205821Sedwin adj360(&LambdaSol); 231205821Sedwin l = 13.1763966 * onday + lzero; /* sec 61 #4 */ 232205821Sedwin adj360(&l); 233205821Sedwin Mm = l - (0.1114041 * onday) - Pzero; /* sec 61 #5 */ 234205821Sedwin adj360(&Mm); 235205821Sedwin Nm = Nzero - (0.0529539 * onday); /* sec 61 #6 */ 236205821Sedwin adj360(&Nm); 237205821Sedwin Ev = 1.2739 * sin(dtor(2*(l - LambdaSol) - Mm)); /* sec 61 #7 */ 238205821Sedwin Ac = 0.1858 * sin(dtor(Msol)); /* sec 61 #8 */ 239205821Sedwin A3 = 0.37 * sin(dtor(Msol)); 240205821Sedwin Mmprime = Mm + Ev - Ac - A3; /* sec 61 #9 */ 241205821Sedwin Ec = 6.2886 * sin(dtor(Mmprime)); /* sec 61 #10 */ 242205821Sedwin A4 = 0.214 * sin(dtor(2 * Mmprime)); /* sec 61 #11 */ 243205821Sedwin lprime = l + Ev + Ec - Ac + A4; /* sec 61 #12 */ 244205821Sedwin V = 0.6583 * sin(dtor(2 * (lprime - LambdaSol))); /* sec 61 #13 */ 245205821Sedwin ldprime = lprime + V; /* sec 61 #14 */ 246205821Sedwin D = ldprime - LambdaSol; /* sec 63 #2 */ 247205821Sedwin return(50 * (1 - cos(dtor(D)))); /* sec 63 #3 */ 248205821Sedwin} 249205821Sedwin 250205821Sedwin/* 251205821Sedwin * dtor -- 252205821Sedwin * convert degrees to radians 253205821Sedwin */ 254205821Sedwinstatic double 255205821Sedwindtor(double deg) 256205821Sedwin{ 257205821Sedwin 258205821Sedwin return(deg * PI / 180); 259205821Sedwin} 260205821Sedwin 261205821Sedwin/* 262205821Sedwin * adj360 -- 263205821Sedwin * adjust value so 0 <= deg <= 360 264205821Sedwin */ 265205821Sedwinstatic void 266205821Sedwinadj360(double *deg) 267205821Sedwin{ 268205821Sedwin 269205821Sedwin for (;;) 270205821Sedwin if (*deg < 0) 271205821Sedwin *deg += 360; 272205821Sedwin else if (*deg > 360) 273205821Sedwin *deg -= 360; 274205821Sedwin else 275205821Sedwin break; 276205821Sedwin} 277