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