1/* 2 * Copyright (c) 1989, 1993 3 * The Regents of the University of California. All rights reserved. 4 * 5 * This code is derived from software posted to USENET. 6 * 7 * Redistribution and use in source and binary forms, with or without 8 * modification, are permitted provided that the following conditions 9 * are met: 10 * 1. Redistributions of source code must retain the above copyright 11 * notice, this list of conditions and the following disclaimer. 12 * 2. Redistributions in binary form must reproduce the above copyright 13 * notice, this list of conditions and the following disclaimer in the 14 * documentation and/or other materials provided with the distribution. 15 * 4. Neither the name of the University nor the names of its contributors 16 * may be used to endorse or promote products derived from this software 17 * without specific prior written permission. 18 * 19 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND 20 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 21 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 22 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE 23 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 24 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 25 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 26 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 27 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 28 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 29 * SUCH DAMAGE. 30 */ 31 32#if 0 33#ifndef lint 34static const char copyright[] = 35"@(#) Copyright (c) 1989, 1993\n\ 36 The Regents of the University of California. All rights reserved.\n"; 37#endif /* not lint */ 38 39#ifndef lint 40static const char sccsid[] = "@(#)pom.c 8.1 (Berkeley) 5/31/93"; 41#endif /* not lint */ 42#endif 43#include <sys/cdefs.h> 44__FBSDID("$FreeBSD$"); 45 46/* 47 * Phase of the Moon. Calculates the current phase of the moon. 48 * Based on routines from `Practical Astronomy with Your Calculator', 49 * by Duffett-Smith. Comments give the section from the book that 50 * particular piece of code was adapted from. 51 * 52 * -- Keith E. Brandt VIII 1984 53 * 54 */ 55 56#include <stdio.h> 57#include <stdlib.h> 58#include <math.h> 59#include <string.h> 60#include <sysexits.h> 61#include <time.h> 62#include <unistd.h> 63 64#include "calendar.h" 65 66#ifndef PI 67#define PI 3.14159265358979323846 68#endif 69#define EPOCH 85 70#define EPSILONg 279.611371 /* solar ecliptic long at EPOCH */ 71#define RHOg 282.680403 /* solar ecliptic long of perigee at EPOCH */ 72#define ECCEN 0.01671542 /* solar orbit eccentricity */ 73#define lzero 18.251907 /* lunar mean long at EPOCH */ 74#define Pzero 192.917585 /* lunar mean long of perigee at EPOCH */ 75#define Nzero 55.204723 /* lunar mean long of node at EPOCH */ 76#define isleap(y) ((((y) % 4) == 0 && ((y) % 100) != 0) || ((y) % 400) == 0) 77 78static void adj360(double *); 79static double dtor(double); 80static double potm(double onday); 81static double potm_minute(double onday, int olddir); 82 83void 84pom(int year, double utcoffset, int *fms, int *nms) 85{ 86 double ffms[MAXMOONS]; 87 double fnms[MAXMOONS]; 88 int i, j; 89 90 fpom(year, utcoffset, ffms, fnms); 91 92 j = 0; 93 for (i = 0; ffms[i] != 0; i++) 94 fms[j++] = round(ffms[i]); 95 fms[i] = -1; 96 for (i = 0; fnms[i] != 0; i++) 97 nms[i] = round(fnms[i]); 98 nms[i] = -1; 99} 100 101void 102fpom(int year, double utcoffset, double *ffms, double *fnms) 103{ 104 time_t tt; 105 struct tm GMT, tmd_today, tmd_tomorrow; 106 double days_today, days_tomorrow, today, tomorrow; 107 int cnt, d; 108 int yeardays; 109 int olddir, newdir; 110 double *pfnms, *pffms, t; 111 112 pfnms = fnms; 113 pffms = ffms; 114 115 /* 116 * We take the phase of the moon one second before and one second 117 * after midnight. 118 */ 119 memset(&tmd_today, 0, sizeof(tmd_today)); 120 tmd_today.tm_year = year - 1900; 121 tmd_today.tm_mon = 0; 122 tmd_today.tm_mday = -1; /* 31 December */ 123 tmd_today.tm_hour = 23; 124 tmd_today.tm_min = 59; 125 tmd_today.tm_sec = 59; 126 memset(&tmd_tomorrow, 0, sizeof(tmd_tomorrow)); 127 tmd_tomorrow.tm_year = year - 1900; 128 tmd_tomorrow.tm_mon = 0; 129 tmd_tomorrow.tm_mday = 0; /* 01 January */ 130 tmd_tomorrow.tm_hour = 0; 131 tmd_tomorrow.tm_min = 0; 132 tmd_tomorrow.tm_sec = 1; 133 134 tt = mktime(&tmd_today); 135 gmtime_r(&tt, &GMT); 136 yeardays = 0; 137 for (cnt = EPOCH; cnt < GMT.tm_year; ++cnt) 138 yeardays += isleap(1900 + cnt) ? DAYSPERLEAPYEAR : DAYSPERYEAR; 139 days_today = (GMT.tm_yday + 1) + ((GMT.tm_hour + 140 (GMT.tm_min / FSECSPERMINUTE) + (GMT.tm_sec / FSECSPERHOUR)) / 141 FHOURSPERDAY); 142 days_today += yeardays; 143 144 tt = mktime(&tmd_tomorrow); 145 gmtime_r(&tt, &GMT); 146 yeardays = 0; 147 for (cnt = EPOCH; cnt < GMT.tm_year; ++cnt) 148 yeardays += isleap(1900 + cnt) ? DAYSPERLEAPYEAR : DAYSPERYEAR; 149 days_tomorrow = (GMT.tm_yday + 1) + ((GMT.tm_hour + 150 (GMT.tm_min / FSECSPERMINUTE) + (GMT.tm_sec / FSECSPERHOUR)) / 151 FHOURSPERDAY); 152 days_tomorrow += yeardays; 153 154 today = potm(days_today); /* 30 December 23:59:59 */ 155 tomorrow = potm(days_tomorrow); /* 31 December 00:00:01 */ 156 olddir = today > tomorrow ? -1 : +1; 157 158 yeardays = 1 + (isleap(year) ? DAYSPERLEAPYEAR : DAYSPERYEAR); /* reuse */ 159 for (d = 0; d <= yeardays; d++) { 160 today = potm(days_today); 161 tomorrow = potm(days_tomorrow); 162 newdir = today > tomorrow ? -1 : +1; 163 if (olddir != newdir) { 164 t = potm_minute(days_today - 1, olddir) + 165 utcoffset / FHOURSPERDAY; 166 if (olddir == -1 && newdir == +1) { 167 *pfnms = d - 1 + t; 168 pfnms++; 169 } else if (olddir == +1 && newdir == -1) { 170 *pffms = d - 1 + t; 171 pffms++; 172 } 173 } 174 olddir = newdir; 175 days_today++; 176 days_tomorrow++; 177 } 178 *pffms = -1; 179 *pfnms = -1; 180} 181 182static double 183potm_minute(double onday, int olddir) { 184 double period = FSECSPERDAY / 2.0; 185 double p1, p2; 186 double before, after; 187 int newdir; 188 189// printf("---> days:%g olddir:%d\n", days, olddir); 190 191 p1 = onday + (period / SECSPERDAY); 192 period /= 2; 193 194 while (period > 30) { /* half a minute */ 195// printf("period:%g - p1:%g - ", period, p1); 196 p2 = p1 + (2.0 / SECSPERDAY); 197 before = potm(p1); 198 after = potm(p2); 199// printf("before:%10.10g - after:%10.10g\n", before, after); 200 newdir = before < after ? -1 : +1; 201 if (olddir != newdir) 202 p1 += (period / SECSPERDAY); 203 else 204 p1 -= (period / SECSPERDAY); 205 period /= 2; 206// printf("newdir:%d - p1:%10.10f - period:%g\n", 207// newdir, p1, period); 208 } 209 p1 -= floor(p1); 210 //exit(0); 211 return (p1); 212} 213 214/* 215 * potm -- 216 * return phase of the moon, as a percentage [0 ... 100] 217 */ 218static double 219potm(double onday) 220{ 221 double N, Msol, Ec, LambdaSol, l, Mm, Ev, Ac, A3, Mmprime; 222 double A4, lprime, V, ldprime, D, Nm; 223 224 N = 360 * onday / 365.2422; /* sec 42 #3 */ 225 adj360(&N); 226 Msol = N + EPSILONg - RHOg; /* sec 42 #4 */ 227 adj360(&Msol); 228 Ec = 360 / PI * ECCEN * sin(dtor(Msol)); /* sec 42 #5 */ 229 LambdaSol = N + Ec + EPSILONg; /* sec 42 #6 */ 230 adj360(&LambdaSol); 231 l = 13.1763966 * onday + lzero; /* sec 61 #4 */ 232 adj360(&l); 233 Mm = l - (0.1114041 * onday) - Pzero; /* sec 61 #5 */ 234 adj360(&Mm); 235 Nm = Nzero - (0.0529539 * onday); /* sec 61 #6 */ 236 adj360(&Nm); 237 Ev = 1.2739 * sin(dtor(2*(l - LambdaSol) - Mm)); /* sec 61 #7 */ 238 Ac = 0.1858 * sin(dtor(Msol)); /* sec 61 #8 */ 239 A3 = 0.37 * sin(dtor(Msol)); 240 Mmprime = Mm + Ev - Ac - A3; /* sec 61 #9 */ 241 Ec = 6.2886 * sin(dtor(Mmprime)); /* sec 61 #10 */ 242 A4 = 0.214 * sin(dtor(2 * Mmprime)); /* sec 61 #11 */ 243 lprime = l + Ev + Ec - Ac + A4; /* sec 61 #12 */ 244 V = 0.6583 * sin(dtor(2 * (lprime - LambdaSol))); /* sec 61 #13 */ 245 ldprime = lprime + V; /* sec 61 #14 */ 246 D = ldprime - LambdaSol; /* sec 63 #2 */ 247 return(50 * (1 - cos(dtor(D)))); /* sec 63 #3 */ 248} 249 250/* 251 * dtor -- 252 * convert degrees to radians 253 */ 254static double 255dtor(double deg) 256{ 257 258 return(deg * PI / 180); 259} 260 261/* 262 * adj360 -- 263 * adjust value so 0 <= deg <= 360 264 */ 265static void 266adj360(double *deg) 267{ 268 269 for (;;) 270 if (*deg < 0) 271 *deg += 360; 272 else if (*deg > 360) 273 *deg -= 360; 274 else 275 break; 276} 277