xref: /freebsd-10.0-release/usr.bin/calendar/sunpos.c

/*-
 * Copyright (c) 2009-2010 Edwin Groothuis. All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 * SUCH DAMAGE.
 *
 */

#include <sys/cdefs.h>
__FBSDID("$FreeBSD: head/usr.bin/calendar/sunpos.c 205821 2010-03-29 06:49:20Z edwin $");

/*
 * This code is created to match the formulas available at:
 * Formula and examples obtained from "How to Calculate alt/az: SAAO" at
 * http://www.saao.ac.za/public-info/sun-moon-stars/sun-index/how-to-calculate-altaz/
 */

#include <stdio.h>
#include <stdlib.h>
#include <limits.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include "calendar.h"

#define D2R(m)	((m) / 180 * M_PI)
#define R2D(m)	((m) * 180 / M_PI)

#define	SIN(x)	(sin(D2R(x)))
#define	COS(x)	(cos(D2R(x)))
#define	TAN(x)	(tan(D2R(x)))
#define	ASIN(x)	(R2D(asin(x)))
#define	ATAN(x)	(R2D(atan(x)))

#ifdef NOTDEF
static void
comp(char *s, double v, double c)
{

	printf("%-*s %*g %*g %*g\n", 15, s, 15, v, 15, c, 15, v - c);
}

int expY;
double expZJ = 30.5;
double expUTHM = 8.5;
double expD = 34743.854;
double expT = 0.9512349;
double expL = 324.885;
double expM = 42.029;
double expepsilon = 23.4396;
double explambda = 326.186;
double expalpha = 328.428;
double expDEC = -12.789;
double expeastlongitude = 17.10;
double explatitude = -22.57;
double expHA = -37.673;
double expALT = 49.822;
double expAZ = 67.49;
#endif

static double
fixup(double *d)
{

	if (*d < 0) {
		while (*d < 0)
			*d += 360;
	} else {
		while (*d > 360)
			*d -= 360;
	}

	return (*d);
}

static double ZJtable[] = {
	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 };

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)
{
	int Y;
	double ZJ, D, T, M, epsilon, lambda, alpha, HA, UTHM;

	ZJ = ZJtable[inMM];
	if (inMM <= 2 && isleap(inYY))
		ZJ -= 1.0;

	UTHM = inHOUR + inMIN / FMINSPERHOUR + inSEC / FSECSPERHOUR - UTCOFFSET;
	Y = inYY - 1900;						/*  1 */
	D = floor(365.25 * Y) + ZJ + inDD + UTHM / FHOURSPERDAY;	/*  3 */
	T = D / 36525.0;						/*  4 */
	*L = 279.697 + 36000.769 * T;					/*  5 */
	fixup(L);
	M = 358.476 + 35999.050 * T;					/*  6 */
	fixup(&M);
	epsilon = 23.452 - 0.013 * T;					/*  7 */
	fixup(&epsilon);

	lambda = *L + (1.919 - 0.005 * T) * SIN(M) + 0.020 * SIN(2 * M);/*  8 */
	fixup(&lambda);
	alpha = ATAN(TAN(lambda) * COS(epsilon));			/*  9 */

	/* Alpha should be in the same quadrant as lamba */
	{
		int lssign = sin(D2R(lambda)) < 0 ? -1 : 1;
		int lcsign = cos(D2R(lambda)) < 0 ? -1 : 1;
		while (((sin(D2R(alpha)) < 0) ? -1 : 1) != lssign
		    || ((cos(D2R(alpha)) < 0) ? -1 : 1) != lcsign)
			alpha += 90.0;
	}
	fixup(&alpha);

	*DEC = ASIN(SIN(lambda) * SIN(epsilon));			/* 10 */
	fixup(DEC);
	fixup(&eastlongitude);
	HA = *L - alpha + 180 + 15 * UTHM + eastlongitude;		/* 12 */
	fixup(&HA);
	fixup(&latitude);
#ifdef NOTDEF
	printf("%02d/%02d %02d:%02d:%02d l:%g d:%g h:%g\n",
	    inMM, inDD, inHOUR, inMIN, inSEC, latitude, *DEC, HA);
#endif
	return;

	/*
	 * The following calculations are not used, so to save time
	 * they are not calculated.
	 */
#ifdef NOTDEF
	*ALT = ASIN(SIN(latitude) * SIN(*DEC) +
	    COS(latitude) * COS(*DEC) * COS(HA));			/* 13 */
	fixup(ALT);
	*AZ = ATAN(SIN(HA) /
	    (COS(HA) * SIN(latitude) - TAN(*DEC) * COS(latitude)));	/* 14 */

	if (*ALT > 180)
		*ALT -= 360;
	if (*ALT < -180)
		*ALT += 360;
	printf("a:%g a:%g\n", *ALT, *AZ);
#endif

#ifdef NOTDEF
	printf("Y:\t\t\t     %d\t\t     %d\t\t      %d\n", Y, expY, Y - expY);
	comp("ZJ", ZJ, expZJ);
	comp("UTHM", UTHM, expUTHM);
	comp("D", D, expD);
	comp("T", T, expT);
	comp("L", L, fixup(&expL));
	comp("M", M, fixup(&expM));
	comp("epsilon", epsilon, fixup(&expepsilon));
	comp("lambda", lambda, fixup(&explambda));
	comp("alpha", alpha, fixup(&expalpha));
	comp("DEC", DEC, fixup(&expDEC));
	comp("eastlongitude", eastlongitude, fixup(&expeastlongitude));
	comp("latitude", latitude, fixup(&explatitude));
	comp("HA", HA, fixup(&expHA));
	comp("ALT", ALT, fixup(&expALT));
	comp("AZ", AZ, fixup(&expAZ));
#endif
}


#define	SIGN(a)	(((a) > 180) ? -1 : 1)
#define ANGLE(a, b) (((a) < (b)) ? 1 : -1)
#define SHOUR(s) ((s) / 3600)
#define SMIN(s) (((s) % 3600) / 60)
#define SSEC(s) ((s) % 60)
#define HOUR(h) ((h) / 4)
#define MIN(h) (15 * ((h) % 4))
#define SEC(h)	0
#define	DEBUG1(y, m, d, hh, mm, pdec, dec) \
	printf("%4d-%02d-%02d %02d:%02d:00 - %7.7g -> %7.7g\n", \
	    y, m, d, hh, mm, pdec, dec)
#define	DEBUG2(y, m, d, hh, mm, pdec, dec, pang, ang) \
	printf("%4d-%02d-%02d %02d:%02d:00 - %7.7g -> %7.7g - %d -> %d\n", \
	    y, m, d, hh, mm, pdec, dec, pang, ang)
void
equinoxsolstice(int year, double UTCoffset, int *equinoxdays, int *solsticedays)
{
	double fe[2], fs[2];

	fequinoxsolstice(year, UTCoffset, fe, fs);
	equinoxdays[0] = round(fe[0]);
	equinoxdays[1] = round(fe[1]);
	solsticedays[0] = round(fs[0]);
	solsticedays[1] = round(fs[1]);
}

void
fequinoxsolstice(int year, double UTCoffset, double *equinoxdays, double *solsticedays)
{
	double dec, prevdec, L;
	int h, d, prevangle, angle;
	int found = 0;

	double decleft, decright, decmiddle;
	int dial, s;

	int *cumdays;
	cumdays = cumdaytab[isleap(year)];

	/*
	 * Find the first equinox, somewhere in March:
	 * It happens when the returned value "dec" goes from
	 * [350 ... 360> -> [0 ... 10]
	 */
	found = 0;
	prevdec = 350;
	for (d = 18; d < 31; d++) {
//		printf("Comparing day %d to %d.\n", d, d+1);
		sunpos(year, 3, d, UTCoffset, 0, 0, 0, 0.0, 0.0, &L, &decleft);
		sunpos(year, 3, d + 1, UTCoffset, 0, 0, 0, 0.0, 0.0,
		    &L, &decright);
//		printf("Found %g and %g.\n", decleft, decright);
		if (SIGN(decleft) == SIGN(decright))
			continue;

		dial = SECSPERDAY;
		s = SECSPERDAY / 2;
		while (s > 0) {
//			printf("Obtaining %d (%02d:%02d)\n",
//			    dial, SHOUR(dial), SMIN(dial));
			sunpos(year, 3, d, UTCoffset,
			    SHOUR(dial), SMIN(dial), SSEC(dial),
			    0.0, 0.0, &L, &decmiddle);
//			printf("Found %g\n", decmiddle);
			if (SIGN(decleft) == SIGN(decmiddle)) {
				decleft = decmiddle;
				dial += s;
			} else {
				decright = decmiddle;
				dial -= s;
			}
//			printf("New boundaries: %g - %g\n", decleft, decright);

			s /= 2;
		}
		equinoxdays[0] = 1 + cumdays[3] + d + (dial / FSECSPERDAY);
		break;
	}

	/* Find the second equinox, somewhere in September:
	 * It happens when the returned value "dec" goes from
	 * [10 ... 0] -> <360 ... 350]
	 */
	found = 0;
	prevdec = 10;
	for (d = 18; d < 31; d++) {
//		printf("Comparing day %d to %d.\n", d, d+1);
		sunpos(year, 9, d, UTCoffset, 0, 0, 0, 0.0, 0.0, &L, &decleft);
		sunpos(year, 9, d + 1, UTCoffset, 0, 0, 0, 0.0, 0.0,
		    &L, &decright);
//		printf("Found %g and %g.\n", decleft, decright);
		if (SIGN(decleft) == SIGN(decright))
			continue;

		dial = SECSPERDAY;
		s = SECSPERDAY / 2;
		while (s > 0) {
//			printf("Obtaining %d (%02d:%02d)\n",
//			    dial, SHOUR(dial), SMIN(dial));
			sunpos(year, 9, d, UTCoffset,
			    SHOUR(dial), SMIN(dial), SSEC(dial),
			    0.0, 0.0, &L, &decmiddle);
//			printf("Found %g\n", decmiddle);
			if (SIGN(decleft) == SIGN(decmiddle)) {
				decleft = decmiddle;
				dial += s;
			} else {
				decright = decmiddle;
				dial -= s;
			}
//			printf("New boundaries: %g - %g\n", decleft, decright);

			s /= 2;
		}
		equinoxdays[1] = 1 + cumdays[9] + d + (dial / FSECSPERDAY);
		break;
	}

	/*
	 * Find the first solstice, somewhere in June:
	 * It happens when the returned value "dec" peaks
	 * [40 ... 45] -> [45 ... 40]
	 */
	found = 0;
	prevdec = 0;
	prevangle = 1;
	for (d = 18; d < 31; d++) {
		for (h = 0; h < 4 * HOURSPERDAY; h++) {
			sunpos(year, 6, d, UTCoffset, HOUR(h), MIN(h), SEC(h),
			    0.0, 0.0, &L, &dec);
			angle = ANGLE(prevdec, dec);
			if (prevangle != angle) {
#ifdef NOTDEF
				DEBUG2(year, 6, d, HOUR(h), MIN(h),
				    prevdec, dec, prevangle, angle);
#endif
				solsticedays[0] = 1 + cumdays[6] + d +
				    ((h / 4.0) / 24.0);
				found = 1;
				break;
			}
			prevdec = dec;
			prevangle = angle;
		}
		if (found)
			break;
	}

	/*
	 * Find the second solstice, somewhere in December:
	 * It happens when the returned value "dec" peaks
	 * [315 ... 310] -> [310 ... 315]
	 */
	found = 0;
	prevdec = 360;
	prevangle = -1;
	for (d = 18; d < 31; d++) {
		for (h = 0; h < 4 * HOURSPERDAY; h++) {
			sunpos(year, 12, d, UTCoffset, HOUR(h), MIN(h), SEC(h),
			    0.0, 0.0, &L, &dec);
			angle = ANGLE(prevdec, dec);
			if (prevangle != angle) {
#ifdef NOTDEF
				DEBUG2(year, 12, d, HOUR(h), MIN(h),
				    prevdec, dec, prevangle, angle);
#endif
				solsticedays[1] = 1 + cumdays[12] + d +
				    ((h / 4.0) / 24.0);
				found = 1;
				break;
			}
			prevdec = dec;
			prevangle = angle;
		}
		if (found)
			break;
	}

	return;
}

int
calculatesunlongitude30(int year, int degreeGMToffset, int *ichinesemonths)
{
	int m, d, h;
	double dec;
	double curL, prevL;
	int *pichinesemonths, *monthdays, *cumdays, i;
	int firstmonth330 = -1;

	cumdays = cumdaytab[isleap(year)];
	monthdays = mondaytab[isleap(year)];
	pichinesemonths = ichinesemonths;

	h = 0;
	sunpos(year - 1, 12, 31,
	    -24 * (degreeGMToffset / 360.0),
	    HOUR(h), MIN(h), SEC(h), 0.0, 0.0, &prevL, &dec);

	for (m = 1; m <= 12; m++) {
		for (d = 1; d <= monthdays[m]; d++) {
			for (h = 0; h < 4 * HOURSPERDAY; h++) {
				sunpos(year, m, d,
				    -24 * (degreeGMToffset / 360.0),
				    HOUR(h), MIN(h), SEC(h),
				    0.0, 0.0, &curL, &dec);
				if (curL < 180 && prevL > 180) {
					*pichinesemonths = cumdays[m] + d;
#ifdef DEBUG
printf("%04d-%02d-%02d %02d:%02d - %d %g\n",
    year, m, d, HOUR(h), MIN(h), *pichinesemonths, curL);
#endif
					    pichinesemonths++;
				} else {
					for (i = 0; i <= 360; i += 30)
						if (curL > i && prevL < i) {
							*pichinesemonths =
							    cumdays[m] + d;
#ifdef DEBUG
printf("%04d-%02d-%02d %02d:%02d - %d %g\n",
    year, m, d, HOUR(h), MIN(h), *pichinesemonths, curL);
#endif
							if (i == 330)
								firstmonth330 = *pichinesemonths;
							pichinesemonths++;
						}
				}
				prevL = curL;
			}
		}
	}
	*pichinesemonths = -1;
	return (firstmonth330);
}

#ifdef NOTDEF
int
main(int argc, char **argv)
{
/*
	year      Mar        June       Sept       Dec
	     day   time  day   time  day time  day time
	2004  20   06:49  21   00:57  22  16:30 21  12:42
	2005  20   12:33  21   06:46  22  22:23 21  18:35
	2006  20   18:26  21   12:26  23  04:03 22  00:22
	2007  21   00:07  21   18:06  23  09:51 22  06:08
	2008  20   05:48  20   23:59  22  15:44 21  12:04
	2009  20   11:44  21   05:45  22  21:18 21  17:47
	2010  20   17:32  21   11:28  23  03:09 21  23:38
	2011  20   23:21  21   17:16  23  09:04 22  05:30
	2012  20   05:14  20   23:09  22  14:49 21  11:11
	2013  20   11:02  21   05:04  22  20:44 21  17:11
	2014  20   16:57  21   10:51  23  02:29 21  23:03
	2015  20   22:45  21   16:38  23  08:20 22  04:48
	2016  20   04:30  20   22:34  22  14:21 21  10:44
	2017  20   10:28  21   04:24  22  20:02 21  16:28
*/

	int eq[2], sol[2];
	equinoxsolstice(strtol(argv[1], NULL, 10), 0.0, eq, sol);
	printf("%d - %d - %d - %d\n", eq[0], sol[0], eq[1], sol[1]);
	return(0);
}
#endif