449 lines
12 KiB
C
449 lines
12 KiB
C
/*-
|
|
* Copyright (c) 2009-2010 Edwin Groothuis <edwin@FreeBSD.org>.
|
|
* 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$");
|
|
|
|
/*
|
|
* 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]
|
|
*/
|
|
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]
|
|
*/
|
|
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
|