1 /*-
2 * SPDX-License-Identifier: BSD-2-Clause-FreeBSD
3 *
4 * Copyright (c) 2009-2010 Edwin Groothuis <[email protected]>.
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
comp(char * s,double v,double c)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
fixup(double * d)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
sunpos(int inYY,int inMM,int inDD,double UTCOFFSET,int inHOUR,int inMIN,int inSEC,double eastlongitude,double latitude,double * L,double * DEC)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
equinoxsolstice(int year,double UTCoffset,int * equinoxdays,int * solsticedays)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
fequinoxsolstice(int year,double UTCoffset,double * equinoxdays,double * solsticedays)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
calculatesunlongitude30(int year,int degreeGMToffset,int * ichinesemonths)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
main(int argc,char ** argv)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