usr.bin/calendar/sunpos.c

   1 /*-
   2  * Copyright (c) 2009-2010 Edwin Groothuis <edwin@FreeBSD.org>.
   3  * All rights reserved.
   4  *
   5  * Redistribution and use in source and binary forms, with or without
   6  * modification, are permitted provided that the following conditions
   7  * are met:
   8  * 1. Redistributions of source code must retain the above copyright
   9  *    notice, this list of conditions and the following disclaimer.
  10  * 2. Redistributions in binary form must reproduce the above copyright
  11  *    notice, this list of conditions and the following disclaimer in the
  12  *    documentation and/or other materials provided with the distribution.
  13  *
  14  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
  15  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  16  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  17  * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
  18  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  19  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  20  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  21  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  22  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  23  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  24  * SUCH DAMAGE.
  25  *
  26  */
  27
  28 #include <sys/cdefs.h>
  29 __FBSDID("$FreeBSD$");
  30
  31 /*
  32  * This code is created to match the formulas available at:
  33  * Formula and examples obtained from "How to Calculate alt/az: SAAO" at
  34  * http://www.saao.ac.za/public-info/sun-moon-stars/sun-index/how-to-calculate-altaz/
  35  */
  36
  37 #include <stdio.h>
  38 #include <stdlib.h>
  39 #include <limits.h>
  40 #include <math.h>
  41 #include <string.h>
  42 #include <time.h>
  43 #include "calendar.h"
  44
  45 #define D2R(m)  ((m) / 180 * M_PI)
  46 #define R2D(m)  ((m) * 180 / M_PI)
  47
  48 #define SIN(x)  (sin(D2R(x)))
  49 #define COS(x)  (cos(D2R(x)))
  50 #define TAN(x)  (tan(D2R(x)))
  51 #define ASIN(x) (R2D(asin(x)))
  52 #define ATAN(x) (R2D(atan(x)))
  53
  54 #ifdef NOTDEF
  55 static void
  56 comp(char *s, double v, double c)
  57 {
  58
  59         printf("%-*s %*g %*g %*g\n", 15, s, 15, v, 15, c, 15, v - c);
  60 }
  61
  62 int expY;
  63 double expZJ = 30.5;
  64 double expUTHM = 8.5;
  65 double expD = 34743.854;
  66 double expT = 0.9512349;
  67 double expL = 324.885;
  68 double expM = 42.029;
  69 double expepsilon = 23.4396;
  70 double explambda = 326.186;
  71 double expalpha = 328.428;
  72 double expDEC = -12.789;
  73 double expeastlongitude = 17.10;
  74 double explatitude = -22.57;
  75 double expHA = -37.673;
  76 double expALT = 49.822;
  77 double expAZ = 67.49;
  78 #endif
  79
  80 static double
  81 fixup(double *d)
  82 {
  83
  84         if (*d < 0) {
  85                 while (*d < 0)
  86                         *d += 360;
  87         } else {
  88                 while (*d > 360)
  89                         *d -= 360;
  90         }
  91
  92         return (*d);
  93 }
  94
  95 static double ZJtable[] = {
  96         0, -0.5, 30.5, 58.5, 89.5, 119.5, 150.5, 180.5, 211.5, 242.5, 272.5, 303.5, 333.5 };
  97
  98 static void
  99 sunpos(int inYY, int inMM, int inDD, double UTCOFFSET, int inHOUR, int inMIN,
 100     int inSEC, double eastlongitude, double latitude, double *L, double *DEC)
 101 {
 102         int Y;
 103         double ZJ, D, T, M, epsilon, lambda, alpha, HA, UTHM;
 104
 105         ZJ = ZJtable[inMM];
 106         if (inMM <= 2 && isleap(inYY))
 107                 ZJ -= 1.0;
 108
 109         UTHM = inHOUR + inMIN / FMINSPERHOUR + inSEC / FSECSPERHOUR - UTCOFFSET;
 110         Y = inYY - 1900;                                                /*  1 */
 111         D = floor(365.25 * Y) + ZJ + inDD + UTHM / FHOURSPERDAY;        /*  3 */
 112         T = D / 36525.0;                                                /*  4 */
 113         *L = 279.697 + 36000.769 * T;                                   /*  5 */
 114         fixup(L);
 115         M = 358.476 + 35999.050 * T;                                    /*  6 */
 116         fixup(&M);
 117         epsilon = 23.452 - 0.013 * T;                                   /*  7 */
 118         fixup(&epsilon);
 119
 120         lambda = *L + (1.919 - 0.005 * T) * SIN(M) + 0.020 * SIN(2 * M);/*  8 */
 121         fixup(&lambda);
 122         alpha = ATAN(TAN(lambda) * COS(epsilon));                       /*  9 */
 123
 124         /* Alpha should be in the same quadrant as lamba */
 125         {
 126                 int lssign = sin(D2R(lambda)) < 0 ? -1 : 1;
 127                 int lcsign = cos(D2R(lambda)) < 0 ? -1 : 1;
 128                 while (((sin(D2R(alpha)) < 0) ? -1 : 1) != lssign
 129                     || ((cos(D2R(alpha)) < 0) ? -1 : 1) != lcsign)
 130                         alpha += 90.0;
 131         }
 132         fixup(&alpha);
 133
 134         *DEC = ASIN(SIN(lambda) * SIN(epsilon));                        /* 10 */
 135         fixup(DEC);
 136         fixup(&eastlongitude);
 137         HA = *L - alpha + 180 + 15 * UTHM + eastlongitude;              /* 12 */
 138         fixup(&HA);
 139         fixup(&latitude);
 140 #ifdef NOTDEF
 141         printf("%02d/%02d %02d:%02d:%02d l:%g d:%g h:%g\n",
 142             inMM, inDD, inHOUR, inMIN, inSEC, latitude, *DEC, HA);
 143 #endif
 144         return;
 145
 146         /*
 147          * The following calculations are not used, so to save time
 148          * they are not calculated.
 149          */
 150 #ifdef NOTDEF
 151         *ALT = ASIN(SIN(latitude) * SIN(*DEC) +
 152             COS(latitude) * COS(*DEC) * COS(HA));                       /* 13 */
 153         fixup(ALT);
 154         *AZ = ATAN(SIN(HA) /
 155             (COS(HA) * SIN(latitude) - TAN(*DEC) * COS(latitude)));     /* 14 */
 156
 157         if (*ALT > 180)
 158                 *ALT -= 360;
 159         if (*ALT < -180)
 160                 *ALT += 360;
 161         printf("a:%g a:%g\n", *ALT, *AZ);
 162 #endif
 163
 164 #ifdef NOTDEF
 165         printf("Y:\t\t\t     %d\t\t     %d\t\t      %d\n", Y, expY, Y - expY);
 166         comp("ZJ", ZJ, expZJ);
 167         comp("UTHM", UTHM, expUTHM);
 168         comp("D", D, expD);
 169         comp("T", T, expT);
 170         comp("L", L, fixup(&expL));
 171         comp("M", M, fixup(&expM));
 172         comp("epsilon", epsilon, fixup(&expepsilon));
 173         comp("lambda", lambda, fixup(&explambda));
 174         comp("alpha", alpha, fixup(&expalpha));
 175         comp("DEC", DEC, fixup(&expDEC));
 176         comp("eastlongitude", eastlongitude, fixup(&expeastlongitude));
 177         comp("latitude", latitude, fixup(&explatitude));
 178         comp("HA", HA, fixup(&expHA));
 179         comp("ALT", ALT, fixup(&expALT));
 180         comp("AZ", AZ, fixup(&expAZ));
 181 #endif
 182 }
 183
 184
 185 #define SIGN(a) (((a) > 180) ? -1 : 1)
 186 #define ANGLE(a, b) (((a) < (b)) ? 1 : -1)
 187 #define SHOUR(s) ((s) / 3600)
 188 #define SMIN(s) (((s) % 3600) / 60)
 189 #define SSEC(s) ((s) % 60)
 190 #define HOUR(h) ((h) / 4)
 191 #define MIN(h) (15 * ((h) % 4))
 192 #define SEC(h)  0
 193 #define DEBUG1(y, m, d, hh, mm, pdec, dec) \
 194         printf("%4d-%02d-%02d %02d:%02d:00 - %7.7g -> %7.7g\n", \
 195             y, m, d, hh, mm, pdec, dec)
 196 #define DEBUG2(y, m, d, hh, mm, pdec, dec, pang, ang) \
 197         printf("%4d-%02d-%02d %02d:%02d:00 - %7.7g -> %7.7g - %d -> %d\n", \
 198             y, m, d, hh, mm, pdec, dec, pang, ang)
 199 void
 200 equinoxsolstice(int year, double UTCoffset, int *equinoxdays, int *solsticedays)
 201 {
 202         double fe[2], fs[2];
 203
 204         fequinoxsolstice(year, UTCoffset, fe, fs);
 205         equinoxdays[0] = round(fe[0]);
 206         equinoxdays[1] = round(fe[1]);
 207         solsticedays[0] = round(fs[0]);
 208         solsticedays[1] = round(fs[1]);
 209 }
 210
 211 void
 212 fequinoxsolstice(int year, double UTCoffset, double *equinoxdays, double *solsticedays)
 213 {
 214         double dec, prevdec, L;
 215         int h, d, prevangle, angle;
 216         int found = 0;
 217
 218         double decleft, decright, decmiddle;
 219         int dial, s;
 220
 221         int *cumdays;
 222         cumdays = cumdaytab[isleap(year)];
 223
 224         /*
 225          * Find the first equinox, somewhere in March:
 226          * It happens when the returned value "dec" goes from
 227          * [350 ... 360> -> [0 ... 10]
 228          */
 229         for (d = 18; d < 31; d++) {
 230                 /* printf("Comparing day %d to %d.\n", d, d+1); */
 231                 sunpos(year, 3, d, UTCoffset, 0, 0, 0, 0.0, 0.0, &L, &decleft);
 232                 sunpos(year, 3, d + 1, UTCoffset, 0, 0, 0, 0.0, 0.0,
 233                     &L, &decright);
 234                 /* printf("Found %g and %g.\n", decleft, decright); */
 235                 if (SIGN(decleft) == SIGN(decright))
 236                         continue;
 237
 238                 dial = SECSPERDAY;
 239                 s = SECSPERDAY / 2;
 240                 while (s > 0) {
 241                         /* printf("Obtaining %d (%02d:%02d)\n",
 242                             dial, SHOUR(dial), SMIN(dial)); */
 243                         sunpos(year, 3, d, UTCoffset,
 244                             SHOUR(dial), SMIN(dial), SSEC(dial),
 245                             0.0, 0.0, &L, &decmiddle);
 246                         /* printf("Found %g\n", decmiddle); */
 247                         if (SIGN(decleft) == SIGN(decmiddle)) {
 248                                 decleft = decmiddle;
 249                                 dial += s;
 250                         } else {
 251                                 decright = decmiddle;
 252                                 dial -= s;
 253                         }
 254                         /*
 255                          printf("New boundaries: %g - %g\n", decleft, decright);
 256                         */
 257
 258                         s /= 2;
 259                 }
 260                 equinoxdays[0] = 1 + cumdays[3] + d + (dial / FSECSPERDAY);
 261                 break;
 262         }
 263
 264         /* Find the second equinox, somewhere in September:
 265          * It happens when the returned value "dec" goes from
 266          * [10 ... 0] -> <360 ... 350]
 267          */
 268         for (d = 18; d < 31; d++) {
 269                 /* printf("Comparing day %d to %d.\n", d, d+1); */
 270                 sunpos(year, 9, d, UTCoffset, 0, 0, 0, 0.0, 0.0, &L, &decleft);
 271                 sunpos(year, 9, d + 1, UTCoffset, 0, 0, 0, 0.0, 0.0,
 272                     &L, &decright);
 273                 /* printf("Found %g and %g.\n", decleft, decright); */
 274                 if (SIGN(decleft) == SIGN(decright))
 275                         continue;
 276
 277                 dial = SECSPERDAY;
 278                 s = SECSPERDAY / 2;
 279                 while (s > 0) {
 280                         /* printf("Obtaining %d (%02d:%02d)\n",
 281                             dial, SHOUR(dial), SMIN(dial)); */
 282                         sunpos(year, 9, d, UTCoffset,
 283                             SHOUR(dial), SMIN(dial), SSEC(dial),
 284                             0.0, 0.0, &L, &decmiddle);
 285                         /* printf("Found %g\n", decmiddle); */
 286                         if (SIGN(decleft) == SIGN(decmiddle)) {
 287                                 decleft = decmiddle;
 288                                 dial += s;
 289                         } else {
 290                                 decright = decmiddle;
 291                                 dial -= s;
 292                         }
 293                         /*
 294                         printf("New boundaries: %g - %g\n", decleft, decright);
 295                         */
 296
 297                         s /= 2;
 298                 }
 299                 equinoxdays[1] = 1 + cumdays[9] + d + (dial / FSECSPERDAY);
 300                 break;
 301         }
 302
 303         /*
 304          * Find the first solstice, somewhere in June:
 305          * It happens when the returned value "dec" peaks
 306          * [40 ... 45] -> [45 ... 40]
 307          */
 308         found = 0;
 309         prevdec = 0;
 310         prevangle = 1;
 311         for (d = 18; d < 31; d++) {
 312                 for (h = 0; h < 4 * HOURSPERDAY; h++) {
 313                         sunpos(year, 6, d, UTCoffset, HOUR(h), MIN(h), SEC(h),
 314                             0.0, 0.0, &L, &dec);
 315                         angle = ANGLE(prevdec, dec);
 316                         if (prevangle != angle) {
 317 #ifdef NOTDEF
 318                                 DEBUG2(year, 6, d, HOUR(h), MIN(h),
 319                                     prevdec, dec, prevangle, angle);
 320 #endif
 321                                 solsticedays[0] = 1 + cumdays[6] + d +
 322                                     ((h / 4.0) / 24.0);
 323                                 found = 1;
 324                                 break;
 325                         }
 326                         prevdec = dec;
 327                         prevangle = angle;
 328                 }
 329                 if (found)
 330                         break;
 331         }
 332
 333         /*
 334          * Find the second solstice, somewhere in December:
 335          * It happens when the returned value "dec" peaks
 336          * [315 ... 310] -> [310 ... 315]
 337          */
 338         found = 0;
 339         prevdec = 360;
 340         prevangle = -1;
 341         for (d = 18; d < 31; d++) {
 342                 for (h = 0; h < 4 * HOURSPERDAY; h++) {
 343                         sunpos(year, 12, d, UTCoffset, HOUR(h), MIN(h), SEC(h),
 344                             0.0, 0.0, &L, &dec);
 345                         angle = ANGLE(prevdec, dec);
 346                         if (prevangle != angle) {
 347 #ifdef NOTDEF
 348                                 DEBUG2(year, 12, d, HOUR(h), MIN(h),
 349                                     prevdec, dec, prevangle, angle);
 350 #endif
 351                                 solsticedays[1] = 1 + cumdays[12] + d +
 352                                     ((h / 4.0) / 24.0);
 353                                 found = 1;
 354                                 break;
 355                         }
 356                         prevdec = dec;
 357                         prevangle = angle;
 358                 }
 359                 if (found)
 360                         break;
 361         }
 362
 363         return;
 364 }
 365
 366 int
 367 calculatesunlongitude30(int year, int degreeGMToffset, int *ichinesemonths)
 368 {
 369         int m, d, h;
 370         double dec;
 371         double curL, prevL;
 372         int *pichinesemonths, *monthdays, *cumdays, i;
 373         int firstmonth330 = -1;
 374
 375         cumdays = cumdaytab[isleap(year)];
 376         monthdays = monthdaytab[isleap(year)];
 377         pichinesemonths = ichinesemonths;
 378
 379         h = 0;
 380         sunpos(year - 1, 12, 31,
 381             -24 * (degreeGMToffset / 360.0),
 382             HOUR(h), MIN(h), SEC(h), 0.0, 0.0, &prevL, &dec);
 383
 384         for (m = 1; m <= 12; m++) {
 385                 for (d = 1; d <= monthdays[m]; d++) {
 386                         for (h = 0; h < 4 * HOURSPERDAY; h++) {
 387                                 sunpos(year, m, d,
 388                                     -24 * (degreeGMToffset / 360.0),
 389                                     HOUR(h), MIN(h), SEC(h),
 390                                     0.0, 0.0, &curL, &dec);
 391                                 if (curL < 180 && prevL > 180) {
 392                                         *pichinesemonths = cumdays[m] + d;
 393 #ifdef DEBUG
 394 printf("%04d-%02d-%02d %02d:%02d - %d %g\n",
 395     year, m, d, HOUR(h), MIN(h), *pichinesemonths, curL);
 396 #endif
 397                                             pichinesemonths++;
 398                                 } else {
 399                                         for (i = 0; i <= 360; i += 30)
 400                                                 if (curL > i && prevL < i) {
 401                                                         *pichinesemonths =
 402                                                             cumdays[m] + d;
 403 #ifdef DEBUG
 404 printf("%04d-%02d-%02d %02d:%02d - %d %g\n",
 405     year, m, d, HOUR(h), MIN(h), *pichinesemonths, curL);
 406 #endif
 407                                                         if (i == 330)
 408                                                                 firstmonth330 = *pichinesemonths;
 409                                                         pichinesemonths++;
 410                                                 }
 411                                 }
 412                                 prevL = curL;
 413                         }
 414                 }
 415         }
 416         *pichinesemonths = -1;
 417         return (firstmonth330);
 418 }
 419
 420 #ifdef NOTDEF
 421 int
 422 main(int argc, char **argv)
 423 {
 424 /*
 425         year      Mar        June       Sept       Dec
 426              day   time  day   time  day time  day time
 427         2004  20   06:49  21   00:57  22  16:30 21  12:42
 428         2005  20   12:33  21   06:46  22  22:23 21  18:35
 429         2006  20   18:26  21   12:26  23  04:03 22  00:22
 430         2007  21   00:07  21   18:06  23  09:51 22  06:08
 431         2008  20   05:48  20   23:59  22  15:44 21  12:04
 432         2009  20   11:44  21   05:45  22  21:18 21  17:47
 433         2010  20   17:32  21   11:28  23  03:09 21  23:38
 434         2011  20   23:21  21   17:16  23  09:04 22  05:30
 435         2012  20   05:14  20   23:09  22  14:49 21  11:11
 436         2013  20   11:02  21   05:04  22  20:44 21  17:11
 437         2014  20   16:57  21   10:51  23  02:29 21  23:03
 438         2015  20   22:45  21   16:38  23  08:20 22  04:48
 439         2016  20   04:30  20   22:34  22  14:21 21  10:44
 440         2017  20   10:28  21   04:24  22  20:02 21  16:28
 441 */
 442
 443         int eq[2], sol[2];
 444         equinoxsolstice(strtol(argv[1], NULL, 10), 0.0, eq, sol);
 445         printf("%d - %d - %d - %d\n", eq[0], sol[0], eq[1], sol[1]);
 446         return(0);
 447 }
 448 #endif