/**************************************************************************** * 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 #endif #include "./suntimes.h" #include #include /* ------------------------------------------------------------------------------------------------ * 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; }