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 451