pom/pom.c

   1 /*      $NetBSD: pom.c,v 1.13 2003/08/07 09:37:32 agc 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\n\
  37         The Regents of the University of California.  All rights reserved.\n");
  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.13 2003/08/07 09:37:32 agc 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 void    adj360 __P((double *));
  88 double  dtor __P((double));
  89 int     main __P((int, char *[]));
  90 double  potm __P((double));
  91 time_t  parsetime __P((char *));
  92 void    badformat __P((void)) __attribute__((__noreturn__));
  93
  94 int
  95 main(argc, argv)
  96         int argc;
  97         char *argv[];
  98 {
  99         time_t tmpt, now;
 100         double days, today, tomorrow;
 101         char buf[1024];
 102
 103         if (time(&now) == (time_t)-1)
 104                 err(1, "time");
 105         if (argc > 1) {
 106                 tmpt = parsetime(argv[1]);
 107                 strftime(buf, sizeof(buf), "%a %Y %b %e %H:%M:%S (%Z)",
 108                         localtime(&tmpt));
 109                 printf("%s:  ", buf);
 110         } else {
 111                 tmpt = now;
 112         }
 113         days = (tmpt - EPOCH_MINUS_1970 * 86400) / 86400.0;
 114         today = potm(days) + .5;
 115         if (tmpt < now)
 116                 (void)printf("The Moon was ");
 117         else if (tmpt == now)
 118                 (void)printf("The Moon is ");
 119         else
 120                 (void)printf("The Moon will be ");
 121         if ((int)today == 100)
 122                 (void)printf("Full\n");
 123         else if (!(int)today)
 124                 (void)printf("New\n");
 125         else {
 126                 tomorrow = potm(days + 1);
 127                 if ((int)today == 50)
 128                         (void)printf("%s\n", tomorrow > today ?
 129                             "at the First Quarter" : "at the Last Quarter");
 130                         /* today is 0.5 too big, but it doesn't matter here
 131                          * since the phase is changing fast enough
 132                          */
 133                 else {
 134                         today -= 0.5;           /* Now it might matter */
 135                         (void)printf("%s ", tomorrow > today ?
 136                             "Waxing" : "Waning");
 137                         if (today > 50)
 138                                 (void)printf("Gibbous (%1.0f%% of Full)\n",
 139                                     today);
 140                         else if (today < 50)
 141                                 (void)printf("Crescent (%1.0f%% of Full)\n",
 142                                     today);
 143                 }
 144         }
 145         exit(0);
 146 }
 147
 148 /*
 149  * potm --
 150  *      return phase of the moon
 151  */
 152 double
 153 potm(days)
 154         double days;
 155 {
 156         double N, Msol, Ec, LambdaSol, l, Mm, Ev, Ac, A3, Mmprime;
 157         double A4, lprime, V, ldprime, D, Nm;
 158
 159         N = 360 * days / 365.242191;                            /* sec 46 #3 */
 160         adj360(&N);
 161         Msol = N + EPSILONg - RHOg;                             /* sec 46 #4 */
 162         adj360(&Msol);
 163         Ec = 360 / PI * ECCEN * sin(dtor(Msol));                /* sec 46 #5 */
 164         LambdaSol = N + Ec + EPSILONg;                          /* sec 46 #6 */
 165         adj360(&LambdaSol);
 166         l = 13.1763966 * days + lzero;                          /* sec 65 #4 */
 167         adj360(&l);
 168         Mm = l - (0.1114041 * days) - Pzero;                    /* sec 65 #5 */
 169         adj360(&Mm);
 170         Nm = Nzero - (0.0529539 * days);                        /* sec 65 #6 */
 171         adj360(&Nm);
 172         Ev = 1.2739 * sin(dtor(2*(l - LambdaSol) - Mm));        /* sec 65 #7 */
 173         Ac = 0.1858 * sin(dtor(Msol));                          /* sec 65 #8 */
 174         A3 = 0.37 * sin(dtor(Msol));
 175         Mmprime = Mm + Ev - Ac - A3;                            /* sec 65 #9 */
 176         Ec = 6.2886 * sin(dtor(Mmprime));                       /* sec 65 #10 */
 177         A4 = 0.214 * sin(dtor(2 * Mmprime));                    /* sec 65 #11 */
 178         lprime = l + Ev + Ec - Ac + A4;                         /* sec 65 #12 */
 179         V = 0.6583 * sin(dtor(2 * (lprime - LambdaSol)));       /* sec 65 #13 */
 180         ldprime = lprime + V;                                   /* sec 65 #14 */
 181         D = ldprime - LambdaSol;                                /* sec 67 #2 */
 182         return(50.0 * (1 - cos(dtor(D))));                      /* sec 67 #3 */
 183 }
 184
 185 /*
 186  * dtor --
 187  *      convert degrees to radians
 188  */
 189 double
 190 dtor(deg)
 191         double deg;
 192 {
 193         return(deg * PI / 180);
 194 }
 195
 196 /*
 197  * adj360 --
 198  *      adjust value so 0 <= deg <= 360
 199  */
 200 void
 201 adj360(deg)
 202         double *deg;
 203 {
 204         for (;;)
 205                 if (*deg < 0)
 206                         *deg += 360;
 207                 else if (*deg > 360)
 208                         *deg -= 360;
 209                 else
 210                         break;
 211 }
 212
 213 #define ATOI2(ar)       ((ar)[0] - '0') * 10 + ((ar)[1] - '0'); (ar) += 2;
 214 time_t
 215 parsetime(p)
 216         char *p;
 217 {
 218         struct tm *lt;
 219         int bigyear;
 220         int yearset = 0;
 221         time_t tval;
 222         unsigned char *t;
 223
 224         for (t = p; *t; ++t) {
 225                 if (isdigit(*t))
 226                         continue;
 227                 badformat();
 228         }
 229
 230         tval = time(NULL);
 231         lt = localtime(&tval);
 232         lt->tm_sec = 0;
 233         lt->tm_min = 0;
 234
 235         switch (strlen(p)) {
 236         case 10:                                /* yyyy */
 237                 bigyear = ATOI2(p);
 238                 lt->tm_year = bigyear * 100 - 1900;
 239                 yearset = 1;
 240                 /* FALLTHROUGH */
 241         case 8:                                 /* yy */
 242                 if (yearset) {
 243                         lt->tm_year += ATOI2(p);
 244                 } else {
 245                         lt->tm_year = ATOI2(p);
 246                         if (lt->tm_year < 69)           /* hack for 2000 */
 247                                 lt->tm_year += 100;
 248                 }
 249                 /* FALLTHROUGH */
 250         case 6:                                 /* mm */
 251                 lt->tm_mon = ATOI2(p);
 252                 if ((lt->tm_mon > 12) || !lt->tm_mon)
 253                         badformat();
 254                 --lt->tm_mon;                   /* time struct is 0 - 11 */
 255                 /* FALLTHROUGH */
 256         case 4:                                 /* dd */
 257                 lt->tm_mday = ATOI2(p);
 258                 if ((lt->tm_mday > 31) || !lt->tm_mday)
 259                         badformat();
 260                 /* FALLTHROUGH */
 261         case 2:                                 /* HH */
 262                 lt->tm_hour = ATOI2(p);
 263                 if (lt->tm_hour > 23)
 264                         badformat();
 265                 break;
 266         default:
 267                 badformat();
 268         }
 269         /* The calling code needs a valid tm_ydays and this is the easiest
 270          * way to get one */
 271         if ((tval = mktime(lt)) == -1)
 272                 errx(1, "specified date is outside allowed range");
 273         return (tval);
 274 }
 275
 276 void
 277 badformat()
 278 {
 279         warnx("illegal time format");
 280         (void)fprintf(stderr, "usage: pom [[[[[cc]yy]mm]dd]HH]\n");
 281         exit(1);
 282 }