usr.bin/calendar/sunpos.c

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