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 449