sunpos.c revision 330449 
1/*-
2 * SPDX-License-Identifier: BSD-2-Clause-FreeBSD
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__FBSDID("$FreeBSD: stable/11/usr.bin/calendar/sunpos.c 330449 2018-03-05 07:26:05Z eadler $");
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
57static void
58comp(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
64int expY;
65double expZJ = 30.5;
66double expUTHM = 8.5;
67double expD = 34743.854;
68double expT = 0.9512349;
69double expL = 324.885;
70double expM = 42.029;
71double expepsilon = 23.4396;
72double explambda = 326.186;
73double expalpha = 328.428;
74double expDEC = -12.789;
75double expeastlongitude = 17.10;
76double explatitude = -22.57;
77double expHA = -37.673;
78double expALT = 49.822;
79double expAZ = 67.49;
80#endif
81
82static double
83fixup(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
97static 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
100static void
101sunpos(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)
201void
202equinoxsolstice(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
213void
214fequinoxsolstice(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
368int
369calculatesunlongitude30(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
396printf("%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
406printf("%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
423int
424main(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