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