2 * SPDX-License-Identifier: BSD-2-Clause-FreeBSD
4 * Copyright (c) 2009-2010 Edwin Groothuis <edwin@FreeBSD.org>.
7 * Redistribution and use in source and binary forms, with or without
8 * modification, are permitted provided that the following conditions
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.
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
30 #include <sys/cdefs.h>
31 __FBSDID("$FreeBSD$");
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/
47 #define D2R(m) ((m) / 180 * M_PI)
48 #define R2D(m) ((m) * 180 / M_PI)
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)))
58 comp(char *s, double v, double c)
61 printf("%-*s %*g %*g %*g\n", 15, s, 15, v, 15, c, 15, v - c);
67 double expD = 34743.854;
68 double expT = 0.9512349;
69 double expL = 324.885;
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;
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 };
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)
105 double ZJ, D, T, M, epsilon, lambda, alpha, HA, UTHM;
108 if (inMM <= 2 && isleap(inYY))
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 */
117 M = 358.476 + 35999.050 * T; /* 6 */
119 epsilon = 23.452 - 0.013 * T; /* 7 */
122 lambda = *L + (1.919 - 0.005 * T) * SIN(M) + 0.020 * SIN(2 * M);/* 8 */
124 alpha = ATAN(TAN(lambda) * COS(epsilon)); /* 9 */
126 /* Alpha should be in the same quadrant as lamba */
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)
136 *DEC = ASIN(SIN(lambda) * SIN(epsilon)); /* 10 */
138 fixup(&eastlongitude);
139 HA = *L - alpha + 180 + 15 * UTHM + eastlongitude; /* 12 */
143 printf("%02d/%02d %02d:%02d:%02d l:%g d:%g h:%g\n",
144 inMM, inDD, inHOUR, inMIN, inSEC, latitude, *DEC, HA);
149 * The following calculations are not used, so to save time
150 * they are not calculated.
153 *ALT = ASIN(SIN(latitude) * SIN(*DEC) +
154 COS(latitude) * COS(*DEC) * COS(HA)); /* 13 */
157 (COS(HA) * SIN(latitude) - TAN(*DEC) * COS(latitude))); /* 14 */
163 printf("a:%g a:%g\n", *ALT, *AZ);
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);
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));
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))
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)
202 equinoxsolstice(int year, double UTCoffset, int *equinoxdays, int *solsticedays)
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]);
214 fequinoxsolstice(int year, double UTCoffset, double *equinoxdays, double *solsticedays)
216 double dec, prevdec, L;
217 int h, d, prevangle, angle;
220 double decleft, decright, decmiddle;
224 cumdays = cumdaytab[isleap(year)];
227 * Find the first equinox, somewhere in March:
228 * It happens when the returned value "dec" goes from
229 * [350 ... 360> -> [0 ... 10]
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,
236 /* printf("Found %g and %g.\n", decleft, decright); */
237 if (SIGN(decleft) == SIGN(decright))
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)) {
253 decright = decmiddle;
257 printf("New boundaries: %g - %g\n", decleft, decright);
262 equinoxdays[0] = 1 + cumdays[3] + d + (dial / FSECSPERDAY);
266 /* Find the second equinox, somewhere in September:
267 * It happens when the returned value "dec" goes from
268 * [10 ... 0] -> <360 ... 350]
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,
275 /* printf("Found %g and %g.\n", decleft, decright); */
276 if (SIGN(decleft) == SIGN(decright))
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)) {
292 decright = decmiddle;
296 printf("New boundaries: %g - %g\n", decleft, decright);
301 equinoxdays[1] = 1 + cumdays[9] + d + (dial / FSECSPERDAY);
306 * Find the first solstice, somewhere in June:
307 * It happens when the returned value "dec" peaks
308 * [40 ... 45] -> [45 ... 40]
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),
317 angle = ANGLE(prevdec, dec);
318 if (prevangle != angle) {
320 DEBUG2(year, 6, d, HOUR(h), MIN(h),
321 prevdec, dec, prevangle, angle);
323 solsticedays[0] = 1 + cumdays[6] + d +
336 * Find the second solstice, somewhere in December:
337 * It happens when the returned value "dec" peaks
338 * [315 ... 310] -> [310 ... 315]
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),
347 angle = ANGLE(prevdec, dec);
348 if (prevangle != angle) {
350 DEBUG2(year, 12, d, HOUR(h), MIN(h),
351 prevdec, dec, prevangle, angle);
353 solsticedays[1] = 1 + cumdays[12] + d +
369 calculatesunlongitude30(int year, int degreeGMToffset, int *ichinesemonths)
374 int *pichinesemonths, *monthdays, *cumdays, i;
375 int firstmonth330 = -1;
377 cumdays = cumdaytab[isleap(year)];
378 monthdays = monthdaytab[isleap(year)];
379 pichinesemonths = ichinesemonths;
382 sunpos(year - 1, 12, 31,
383 -24 * (degreeGMToffset / 360.0),
384 HOUR(h), MIN(h), SEC(h), 0.0, 0.0, &prevL, &dec);
386 for (m = 1; m <= 12; m++) {
387 for (d = 1; d <= monthdays[m]; d++) {
388 for (h = 0; h < 4 * HOURSPERDAY; h++) {
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;
396 printf("%04d-%02d-%02d %02d:%02d - %d %g\n",
397 year, m, d, HOUR(h), MIN(h), *pichinesemonths, curL);
401 for (i = 0; i <= 360; i += 30)
402 if (curL > i && prevL < i) {
406 printf("%04d-%02d-%02d %02d:%02d - %d %g\n",
407 year, m, d, HOUR(h), MIN(h), *pichinesemonths, curL);
410 firstmonth330 = *pichinesemonths;
418 *pichinesemonths = -1;
419 return (firstmonth330);
424 main(int argc, char **argv)
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
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]);