pom/pom.c

   1 /*      $NetBSD: pom.c,v 1.20 2010/12/05 04:34:22 pgoyette Exp $        */
   2
   3 /*
   4  * Copyright (c) 1989, 1993
   5  *      The Regents of the University of California.  All rights reserved.
   6  *
   7  * This code is derived from software posted to USENET.
   8  *
   9  * Redistribution and use in source and binary forms, with or without
  10  * modification, are permitted provided that the following conditions
  11  * are met:
  12  * 1. Redistributions of source code must retain the above copyright
  13  *    notice, this list of conditions and the following disclaimer.
  14  * 2. Redistributions in binary form must reproduce the above copyright
  15  *    notice, this list of conditions and the following disclaimer in the
  16  *    documentation and/or other materials provided with the distribution.
  17  * 3. Neither the name of the University nor the names of its contributors
  18  *    may be used to endorse or promote products derived from this software
  19  *    without specific prior written permission.
  20  *
  21  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  22  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  23  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  24  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  25  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  26  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  27  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  28  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  29  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  30  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  31  * SUCH DAMAGE.
  32  */
  33
  34 #include <sys/cdefs.h>
  35 #ifndef lint
  36 __COPYRIGHT("@(#) Copyright (c) 1989, 1993\
  37  The Regents of the University of California.  All rights reserved.");
  38 #endif /* not lint */
  39
  40 #ifndef lint
  41 #if 0
  42 static char sccsid[] = "@(#)pom.c       8.1 (Berkeley) 5/31/93";
  43 #else
  44 __RCSID("$NetBSD: pom.c,v 1.20 2010/12/05 04:34:22 pgoyette Exp $");
  45 #endif
  46 #endif /* not lint */
  47
  48 /*
  49  * Phase of the Moon.  Calculates the current phase of the moon.
  50  * Based on routines from `Practical Astronomy with Your Calculator',
  51  * by Duffett-Smith.  Comments give the section from the book that
  52  * particular piece of code was adapted from.
  53  *
  54  * -- Keith E. Brandt  VIII 1984
  55  *
  56  * Updated to the Third Edition of Duffett-Smith's book, Paul Janzen, IX 1998
  57  *
  58  */
  59
  60 #include <ctype.h>
  61 #include <err.h>
  62 #include <math.h>
  63 #include <stdio.h>
  64 #include <string.h>
  65 #include <stdlib.h>
  66 #include <time.h>
  67 #include <unistd.h>
  68
  69 #ifndef PI
  70 #define PI        3.14159265358979323846
  71 #endif
  72
  73 /*
  74  * The EPOCH in the third edition of the book is 1990 Jan 0.0 TDT.
  75  * In this program, we do not bother to correct for the differences
  76  * between UTC (as shown by the UNIX clock) and TDT.  (TDT = TAI + 32.184s;
  77  * TAI-UTC = 32s in Jan 1999.)
  78  */
  79 #define EPOCH_MINUS_1970        (20 * 365 + 5 - 1) /* 20 years, 5 leaps, back 1 day to Jan 0 */
  80 #define EPSILONg  279.403303    /* solar ecliptic long at EPOCH */
  81 #define RHOg      282.768422    /* solar ecliptic long of perigee at EPOCH */
  82 #define ECCEN     0.016713      /* solar orbit eccentricity */
  83 #define lzero     318.351648    /* lunar mean long at EPOCH */
  84 #define Pzero     36.340410     /* lunar mean long of perigee at EPOCH */
  85 #define Nzero     318.510107    /* lunar mean long of node at EPOCH */
  86
  87 int     main(int, char *[]);
  88 static void adj360(double *);
  89 static double dtor(double);
  90 static double potm(double);
  91 static time_t parsetime(char *);
  92 static void badformat(void) __dead;
  93
  94 int
  95 main(int argc, char *argv[])
  96 {
  97         time_t tmpt, now;
  98         double days, today, tomorrow;
  99         char buf[1024];
 100
 101         if (time(&now) == (time_t)-1)
 102                 err(1, "time");
 103         if (argc > 1) {
 104                 tmpt = parsetime(argv[1]);
 105                 strftime(buf, sizeof(buf), "%a %Y %b %e %H:%M:%S (%Z)",
 106                         localtime(&tmpt));
 107                 printf("%s:  ", buf);
 108         } else {
 109                 tmpt = now;
 110         }
 111         days = (tmpt - EPOCH_MINUS_1970 * 86400) / 86400.0;
 112         today = potm(days) + .5;
 113         if (tmpt < now)
 114                 (void)printf("The Moon was ");
 115         else if (tmpt == now)
 116                 (void)printf("The Moon is ");
 117         else
 118                 (void)printf("The Moon will be ");
 119         if ((int)today == 100)
 120                 (void)printf("Full\n");
 121         else if (!(int)today)
 122                 (void)printf("New\n");
 123         else {
 124                 tomorrow = potm(days + 1);
 125                 if ((int)today == 50)
 126                         (void)printf("%s\n", tomorrow > today ?
 127                             "at the First Quarter" : "at the Last Quarter");
 128                         /* today is 0.5 too big, but it doesn't matter here
 129                          * since the phase is changing fast enough
 130                          */
 131                 else {
 132                         today -= 0.5;           /* Now it might matter */
 133                         (void)printf("%s ", tomorrow > today ?
 134                             "Waxing" : "Waning");
 135                         if (today > 50)
 136                                 (void)printf("Gibbous (%1.0f%% of Full)\n",
 137                                     today);
 138                         else if (today < 50)
 139                                 (void)printf("Crescent (%1.0f%% of Full)\n",
 140                                     today);
 141                 }
 142         }
 143
 144         return EXIT_SUCCESS;
 145 }
 146
 147 /*
 148  * potm --
 149  *      return phase of the moon
 150  */
 151 static double
 152 potm(double days)
 153 {
 154         double N, Msol, Ec, LambdaSol, l, Mm, Ev, Ac, A3, Mmprime;
 155         double A4, lprime, V, ldprime, D, Nm;
 156
 157         N = 360 * days / 365.242191;                            /* sec 46 #3 */
 158         adj360(&N);
 159         Msol = N + EPSILONg - RHOg;                             /* sec 46 #4 */
 160         adj360(&Msol);
 161         Ec = 360 / PI * ECCEN * sin(dtor(Msol));                /* sec 46 #5 */
 162         LambdaSol = N + Ec + EPSILONg;                          /* sec 46 #6 */
 163         adj360(&LambdaSol);
 164         l = 13.1763966 * days + lzero;                          /* sec 65 #4 */
 165         adj360(&l);
 166         Mm = l - (0.1114041 * days) - Pzero;                    /* sec 65 #5 */
 167         adj360(&Mm);
 168         Nm = Nzero - (0.0529539 * days);                        /* sec 65 #6 */
 169         adj360(&Nm);
 170         Ev = 1.2739 * sin(dtor(2*(l - LambdaSol) - Mm));        /* sec 65 #7 */
 171         Ac = 0.1858 * sin(dtor(Msol));                          /* sec 65 #8 */
 172         A3 = 0.37 * sin(dtor(Msol));
 173         Mmprime = Mm + Ev - Ac - A3;                            /* sec 65 #9 */
 174         Ec = 6.2886 * sin(dtor(Mmprime));                       /* sec 65 #10 */
 175         A4 = 0.214 * sin(dtor(2 * Mmprime));                    /* sec 65 #11 */
 176         lprime = l + Ev + Ec - Ac + A4;                         /* sec 65 #12 */
 177         V = 0.6583 * sin(dtor(2 * (lprime - LambdaSol)));       /* sec 65 #13 */
 178         ldprime = lprime + V;                                   /* sec 65 #14 */
 179         D = ldprime - LambdaSol;                                /* sec 67 #2 */
 180         return(50.0 * (1 - cos(dtor(D))));                      /* sec 67 #3 */
 181 }
 182
 183 /*
 184  * dtor --
 185  *      convert degrees to radians
 186  */
 187 static double
 188 dtor(double deg)
 189 {
 190         return(deg * PI / 180);
 191 }
 192
 193 /*
 194  * adj360 --
 195  *      adjust value so 0 <= deg <= 360
 196  */
 197 static void
 198 adj360(double *deg)
 199 {
 200         for (;;)
 201                 if (*deg < 0)
 202                         *deg += 360;
 203                 else if (*deg > 360)
 204                         *deg -= 360;
 205                 else
 206                         break;
 207 }
 208
 209 #define ATOI2(ar)       ((ar)[0] - '0') * 10 + ((ar)[1] - '0'); (ar) += 2;
 210 static time_t
 211 parsetime(char *p)
 212 {
 213         struct tm *lt;
 214         int bigyear;
 215         int yearset = 0;
 216         time_t tval;
 217         char *t;
 218
 219         for (t = p; *t; ++t) {
 220                 if (isdigit((unsigned char) *t))
 221                         continue;
 222                 badformat();
 223         }
 224
 225         tval = time(NULL);
 226         lt = localtime(&tval);
 227         lt->tm_sec = 0;
 228         lt->tm_min = 0;
 229
 230         switch (strlen(p)) {
 231         case 10:                                /* yyyy */
 232                 bigyear = ATOI2(p);
 233                 lt->tm_year = bigyear * 100 - 1900;
 234                 yearset = 1;
 235                 /* FALLTHROUGH */
 236         case 8:                                 /* yy */
 237                 if (yearset) {
 238                         lt->tm_year += ATOI2(p);
 239                 } else {
 240                         lt->tm_year = ATOI2(p);
 241                         if (lt->tm_year < 69)           /* hack for 2000 */
 242                                 lt->tm_year += 100;
 243                 }
 244                 /* FALLTHROUGH */
 245         case 6:                                 /* mm */
 246                 lt->tm_mon = ATOI2(p);
 247                 if ((lt->tm_mon > 12) || !lt->tm_mon)
 248                         badformat();
 249                 --lt->tm_mon;                   /* time struct is 0 - 11 */
 250                 /* FALLTHROUGH */
 251         case 4:                                 /* dd */
 252                 lt->tm_mday = ATOI2(p);
 253                 if ((lt->tm_mday > 31) || !lt->tm_mday)
 254                         badformat();
 255                 /* FALLTHROUGH */
 256         case 2:                                 /* HH */
 257                 lt->tm_hour = ATOI2(p);
 258                 if (lt->tm_hour > 23)
 259                         badformat();
 260                 break;
 261         default:
 262                 badformat();
 263         }
 264         /* The calling code needs a valid tm_ydays and this is the easiest
 265          * way to get one */
 266         if ((tval = mktime(lt)) == -1)
 267                 errx(1, "specified date is outside allowed range");
 268         return (tval);
 269 }
 270
 271 static void
 272 badformat(void)
 273 {
 274         warnx("illegal time format");
 275         (void)fprintf(stderr, "usage: %s [[[[[cc]yy]mm]dd]HH]\n",
 276             getprogname());
 277         exit(EXIT_FAILURE);
 278 }