204 lines
6.1 KiB
C
204 lines
6.1 KiB
C
/****************************************************************************
|
|
* This file is part of the project AqHome.
|
|
* AqHome (c) by 2024 Martin Preuss, all rights reserved.
|
|
*
|
|
* The license for this file can be found in the file COPYING which you
|
|
* should have received along with this file.
|
|
****************************************************************************/
|
|
|
|
#ifdef HAVE_CONFIG_H
|
|
# include <config.h>
|
|
#endif
|
|
|
|
#include "./suntimes.h"
|
|
|
|
#include <stdio.h>
|
|
#include <math.h>
|
|
|
|
|
|
|
|
/* ------------------------------------------------------------------------------------------------
|
|
* forward declarations
|
|
* ------------------------------------------------------------------------------------------------
|
|
*/
|
|
|
|
static GWEN_TIME *_getSunTime(const GWEN_DATE *date, double lat, double lng, int calcSunrise);
|
|
static double _calcSunTime(const GWEN_DATE *date, double lat, double lng, int calcSunrise);
|
|
|
|
|
|
|
|
/* ------------------------------------------------------------------------------------------------
|
|
* implementations
|
|
* ------------------------------------------------------------------------------------------------
|
|
*/
|
|
|
|
GWEN_TIME *AQHomeReact_GetSunriseTimeForDateAndLoc(const GWEN_DATE *date, double lat, double lng)
|
|
{
|
|
return _getSunTime(date, lat, lng, 1);
|
|
}
|
|
|
|
|
|
GWEN_TIME *AQHomeReact_GetSunsetTimeForDateAndLoc(const GWEN_DATE *date, double lat, double lng)
|
|
{
|
|
return _getSunTime(date, lat, lng, 0);
|
|
}
|
|
|
|
|
|
GWEN_TIME *_getSunTime(const GWEN_DATE *date, double lat, double lng, int calcSunrise)
|
|
{
|
|
double timeInFloat=fmod(24+_calcSunTime(date, lat, lng, calcSunrise), 24);
|
|
double hours;
|
|
double minutes = modf(timeInFloat, &hours)*60;
|
|
GWEN_TIME *ti=GWEN_Time_new(GWEN_Date_GetYear(date), GWEN_Date_GetMonth(date)-1, GWEN_Date_GetDay(date),
|
|
hours, minutes, 0, 1);
|
|
return ti;
|
|
}
|
|
|
|
|
|
|
|
#if 0
|
|
int _test(void) {
|
|
double lat = 52.619425; // Latitude of location.
|
|
double lng = 10.087891; // Longitude of location.
|
|
int y=2024;
|
|
int m=4;
|
|
int d;
|
|
|
|
for (d=19; d<=24; d++) {
|
|
GWEN_DATE *date;
|
|
GWEN_TIME *sunrise;
|
|
GWEN_TIME *sunset;
|
|
int sunsetHours, sunsetMinutes, sunsetSecs;
|
|
int sunriseHours, sunriseMinutes, sunriseSecs;
|
|
|
|
date=GWEN_Date_fromGregorian(y, m, d);
|
|
sunrise=AQHomeReact_GetSunriseTimeForDateAndLoc(date, lat, lng);
|
|
sunset=AQHomeReact_GetSunsetTimeForDateAndLoc(date, lat, lng);
|
|
|
|
GWEN_Time_GetBrokenDownTime(sunset, &sunsetHours, &sunsetMinutes, &sunsetSecs);
|
|
GWEN_Time_GetBrokenDownTime(sunrise, &sunriseHours, &sunriseMinutes, &sunriseSecs);
|
|
|
|
fprintf(stdout, "date=%04d-%02d-%02d: sunrise=%02d:%02d sunset=%02d:%02d\n",
|
|
GWEN_Date_GetYear(date), GWEN_Date_GetMonth(date), GWEN_Date_GetDay(date),
|
|
sunriseHours, sunriseMinutes,
|
|
sunsetHours, sunsetMinutes);
|
|
GWEN_Time_free(sunset);
|
|
GWEN_Time_free(sunrise);
|
|
GWEN_Date_free(date);
|
|
}
|
|
}
|
|
#endif
|
|
|
|
|
|
|
|
/**
|
|
* Calculate time of sunset or sunrise at a given position and date according to an algorithm presented
|
|
* on this site:
|
|
* http://edwilliams.org/sunrise_sunset_algorithm.htm
|
|
*
|
|
* Original notes follow.
|
|
*
|
|
* -----------------------------
|
|
* Source:
|
|
* Almanac for Computers, 1990
|
|
* published by Nautical Almanac Office
|
|
* United States Naval Observatory
|
|
* Washington, DC 20392
|
|
*
|
|
* Inputs:
|
|
* date : date for which to calculate sun time
|
|
* lat, lng: latitude and longitude of the location to calculate for (west and north are positive values)
|
|
* ZENITH : Sun's zenith for sunrise/sunset
|
|
* official = 90 degrees 50'
|
|
* civil = 96 degrees
|
|
* nautical = 102 degrees
|
|
* astronomical = 108 degrees
|
|
*
|
|
* NOTE: longitude is positive for East and negative for West
|
|
* NOTE: the algorithm assumes the use of a calculator with the
|
|
* trig functions in "degree" (rather than "radian") mode. Most
|
|
* programming languages assume radian arguments, requiring back
|
|
* and forth convertions. The factor is 180/pi. So, for instance,
|
|
* the equation RA = atan(0.91764 * tan(L)) would be coded as RA
|
|
* = (180/pi)*atan(0.91764 * tan((pi/180)*L)) to give a degree
|
|
* answer with a degree input for L.
|
|
*/
|
|
double _calcSunTime(const GWEN_DATE *date, double lat, double lng, int calcSunrise)
|
|
{
|
|
const double PI=3.1415926535897932384626433832795;
|
|
const double ZENITH=-.83; /* accounts for atmospheric refraction (assuming height of 0m above horizon) */
|
|
double PiBy180=PI/180.0;
|
|
double N;
|
|
double lngHour;
|
|
double t;
|
|
double M;
|
|
double L;
|
|
double RA;
|
|
double Lquadrant;
|
|
double RAquadrant;
|
|
double sinDec;
|
|
double cosDec;
|
|
double cosH;
|
|
double H;
|
|
double T;
|
|
double UT;
|
|
|
|
/* 1. first calculate the day of the year */
|
|
N=GWEN_Date_DaysInYear(date);
|
|
|
|
/* 2. convert the longitude to hour value and calculate an approximate time */
|
|
lngHour = lng/15.0;
|
|
if (calcSunrise)
|
|
t = N+((6-lngHour)/24);
|
|
else
|
|
t = N+((18-lngHour)/24);
|
|
|
|
/* 3. calculate the Sun's mean anomaly */
|
|
M = (0.9856*t)-3.289;
|
|
|
|
/* 4. calculate the Sun's true longitude */
|
|
L = fmod(M+(1.916*sin(PiBy180*M)) + (0.020 * sin(2*PiBy180*M)) + 282.634, 360.0);
|
|
|
|
/* 5a. calculate the Sun's right ascension */
|
|
RA = fmod(180/PI*atan(0.91764*tan(PiBy180*L)), 360.0);
|
|
|
|
/* 5b. right ascension value needs to be in the same quadrant as L */
|
|
Lquadrant = floor(L/90)*90;
|
|
RAquadrant = floor(RA/90)*90;
|
|
RA = RA + (Lquadrant-RAquadrant);
|
|
|
|
/* 5c. right ascension value needs to be converted into hours */
|
|
RA = RA / 15;
|
|
|
|
/* 6. calculate the Sun's declination */
|
|
sinDec = 0.39782*sin(PiBy180*L);
|
|
cosDec = cos(asin(sinDec));
|
|
|
|
/* 7a. calculate the Sun's local hour angle */
|
|
cosH = (sin(PiBy180*ZENITH)-(sinDec*sin(PiBy180*lat))) / (cosDec*cos(PiBy180*lat));
|
|
/* cosH >1 : sun never rises on this location and date
|
|
* cosH <-1: sun never sets on this location and date
|
|
*/
|
|
|
|
/* 7b. finish calculating H and convert into hours */
|
|
if (calcSunrise)
|
|
H = 360-(180/PI)*acos(cosH);
|
|
else
|
|
H = (180/PI)*acos(cosH);
|
|
H = H/15;
|
|
|
|
/* 8. calculate local mean time of rising/setting */
|
|
T = H+RA-(0.06571*t)-6.622;
|
|
|
|
/* 9. adjust back to UTC */
|
|
UT = fmod(T-lngHour, 24.0);
|
|
|
|
/* we don't do step 10. here which would be to convert UT into local time as we want UTC */
|
|
return UT;
|
|
}
|
|
|
|
|
|
|
|
|
|
|