I have a problem with calculating the sun cycles, based on the current time and the latitude and longitude data.
I'm using an ESP32 with Arduino framework. My suntimes always offset for one or one and a half hour. I got my geolocation data from an external API, which will give me my latitude and longitude. This is somewhat reliable, it puts me in the same city where iam.
The current time is coming from an NTP server, if there is no NTP, i have an RTC module alongside my ESP32 and i got the time from there. It is 98% accurate that way.
So the calculation. My current dst offset in my country is 0. If it would be +1 hour now, the calculation would be almost correct. ( a couple of minutes off )
So i start with requesting my geo data from the API like this:
void gSystem::checkConfigGeoData(){
if( hsh_timeSystem.timeIsOk() && gotGeoData && hsh_fileSystem.config.lastGeoDay != hsh_timeSystem.getDayOfWeek() ){
gotGeoData = false;
sunTimesCalculated = false;
geoAttempts = 0;
}
}
void gSystem::getInfo(){
checkConfigGeoData();
if( !gotGeoData && hsh_networkSystem.isConnected() && hsh_networkSystem.connectionMode != AP_CONNECTION &&
( geoAttempts < MAX_GEO_ATTEMPT_COUNT ) && ( millis() - lastGeoAttemptMS >= GEO_ATTEMPT_RETRY_MS ) ){
lastGeoAttemptMS = millis();
if( hsh_fileSystem.config.latitude != 0 && hsh_timeSystem.getUnixSec() - hsh_fileSystem.config.lastGeoEpoch <= GEO_EPOCH_SEC_DAY ){
#if GEO_DEBUG_MODE
Serial.printf("Geo data has been requested less then a day. Abort request...\n");
#endif
gotGeoData = true;
sunTimesCalculated = false;
geoAttempts = MAX_GEO_ATTEMPT_COUNT + 1;
return;
}
geoAttempts++;
HTTPClient http;
http.begin( ipStackURL ); // url is coming from there: http://api.ipstack.com/
int httpResponseCode = http.GET();
if (httpResponseCode > 0) {
DynamicJsonDocument doc(GEO_RESPONSE_JSON_SIZE);
DeserializationError error = deserializeJson(doc,http.getString());
if(error){
char errorMsg[128];
sprintf(errorMsg,"geolocateAPI response deserialization failed with error: %s",error.c_str());
hsh_fileSystem.logToFile(errorMsg,"error","gSystem");
#if GEO_DEBUG_MODE
Serial.println( errorMsg );
#endif
}else{
if( doc["success"] == false ){
char errorMsg[300];
sprintf( errorMsg,"Error: %s",doc["error"]["info"].as<const char*>() );
hsh_fileSystem.logToFile(errorMsg,"error","gSystem" );
#if GEO_DEBUG_MODE
Serial.println( errorMsg );
#endif
}else{
sunTimesCalculated = false;
gotGeoData = true;
hsh_fileSystem.config.publicIP.fromString(doc["ip"].as<const char*>());
strncpy(hsh_fileSystem.config.city, doc["city"], sizeof(hsh_fileSystem.config.city));
hsh_fileSystem.config.latitude = doc["latitude"].as<float>();
hsh_fileSystem.config.longitude = doc["longitude"].as<float>();
hsh_fileSystem.config.lastGeoDay = hsh_timeSystem.getDayOfWeek();
hsh_fileSystem.config.lastGeoEpoch = hsh_timeSystem.getUnixSec();
#if GEO_DEBUG_MODE
Serial.printf("\n**** GOT GEO DATA ****\n");
Serial.printf("GS - City: %s\n",hsh_fileSystem.config.city);
Serial.printf("GS - Public IP: %s\n",hsh_fileSystem.config.publicIP.toString().c_str());
Serial.printf("GS - Latitude: %f\n",hsh_fileSystem.config.latitude);
Serial.printf("GS - Longitude: %f\n\n",hsh_fileSystem.config.longitude);
#endif
hsh_fileSystem.makeConfig();
geoAttempts = 0;
}
}
}
}
}
If i got the data and my time is okay, i start to calculate the sunrise and sunset with this function:
void gSystem::startCalculateSunTimes(){
if( hsh_timeSystem.timeIsOk() && gotGeoData && !sunTimesCalculated ){
boolean timesIsOk = true;
time_t seconds;
time_t tseconds;
struct tm *ptm = NULL;
struct tm tm;
int year = hsh_timeSystem.getYear(),
month = hsh_timeSystem.getMonth(),
day = hsh_timeSystem.getDayOfMonth(),
hour = hsh_timeSystem.getHour(),
min = hsh_timeSystem.getMinute(),
sec = hsh_timeSystem.getSecond();
float JD = calcJD(year, month, day);
tm.tm_year = year - 1900;
tm.tm_mon = month - 1;
tm.tm_mday = day;
tm.tm_hour = hour; //0;
tm.tm_min = min; //0;
tm.tm_sec = sec; //0;
#if GEO_DEBUG_MODE
Serial.printf("\n**** START CALCULATE SUN TIMES ****\n");
Serial.printf("GS - JD: %f\n",JD);
Serial.printf("GS - Time: %d-%d-%d %d:%d:%d\n",year,month,day,hour,min,sec);
Serial.printf("GS - Latitude: %f\n",hsh_fileSystem.config.latitude);
Serial.printf("GS - Longitude: %f\n",hsh_fileSystem.config.longitude);
#endif
seconds = mktime(&tm);
int delta;
ptm = gmtime(&seconds);
delta = ptm->tm_hour;
tseconds = seconds;
seconds = seconds + calcSunriseUTC(JD, hsh_fileSystem.config.latitude, -hsh_fileSystem.config.longitude) * 60;
seconds = seconds - delta * 3600;
ptm = gmtime(&seconds);
int calculatedYear = ptm->tm_year + 1900;
if( calculatedYear == hsh_timeSystem.getYear() ){
hsh_fileSystem.config.sunRiseHour = ptm->tm_hour + hsh_fileSystem.config.dst;
hsh_fileSystem.config.sunRiseMinute = ptm->tm_min;
}else{
timesIsOk = false;
}
seconds = tseconds;
seconds += calcSunsetUTC(JD, hsh_fileSystem.config.latitude, -hsh_fileSystem.config.longitude) * 60;
seconds = seconds - delta * 3600;
ptm = gmtime(&seconds);
calculatedYear = ptm->tm_year + 1900;
if( calculatedYear == hsh_timeSystem.getYear() ){
hsh_fileSystem.config.sunSetHour = ptm->tm_hour + hsh_fileSystem.config.dst;
hsh_fileSystem.config.sunSetMinute = ptm->tm_min;
}else{
timesIsOk = false;
}
if(timesIsOk){
sunTimesCalculated = true;
hsh_fileSystem.makeConfig();
#if GEO_DEBUG_MODE
Serial.println("\n**** Sun times calculated ****");
Serial.printf("GS - Sun Rise Info: %02d:%02d\n",hsh_fileSystem.config.sunRiseHour,hsh_fileSystem.config.sunRiseMinute);
Serial.printf("GS - Sun Set Info: %02d:%02d\n\n",hsh_fileSystem.config.sunSetHour,hsh_fileSystem.config.sunSetMinute);
#endif
}
}
}
This two function is inside a task, and an infinite loop. like this:
void gSystemLoopTask(void* parameter) {
for (;;) {
hsh_GeoSystem.getInfo();
hsh_GeoSystem.startCalculateSunTimes();
vTaskDelay(1000);
}
}
The rest of the calculation is here ( i got that from an open source C project ):
/* Convert degree angle to radians */
double gSystem::degToRad(double angleDeg) {
return (PI * angleDeg / 180.0);
}
double gSystem::radToDeg(double angleRad) {
return (180.0 * angleRad / PI);
}
double gSystem::calcMeanObliquityOfEcliptic(double t) {
double seconds = 21.448 - t * (46.8150 + t * (0.00059 - t * (0.001813)));
double e0 = 23.0 + (26.0 + (seconds / 60.0)) / 60.0;
return e0; // in degrees
}
double gSystem::calcGeomMeanLongSun(double t) {
double L = 280.46646 + t * (36000.76983 + 0.0003032 * t);
while ((int)L > 360) {
L -= 360.0;
}
while (L < 0) {
L += 360.0;
}
return L; // in degrees
}
double gSystem::calcObliquityCorrection(double t) {
double e0 = calcMeanObliquityOfEcliptic(t);
double omega = 125.04 - 1934.136 * t;
double e = e0 + 0.00256 * cos(degToRad(omega));
return e; // in degrees
}
double gSystem::calcEccentricityEarthOrbit(double t) {
double e = 0.016708634 - t * (0.000042037 + 0.0000001267 * t);
return e; // unitless
}
double gSystem::calcGeomMeanAnomalySun(double t) {
double M = 357.52911 + t * (35999.05029 - 0.0001537 * t);
return M; // in degrees
}
double gSystem::calcEquationOfTime(double t) {
double epsilon = calcObliquityCorrection(t);
double l0 = calcGeomMeanLongSun(t);
double e = calcEccentricityEarthOrbit(t);
double m = calcGeomMeanAnomalySun(t);
double y = tan(degToRad(epsilon) / 2.0);
y *= y;
double sin2l0 = sin(2.0 * degToRad(l0));
double sinm = sin(degToRad(m));
double cos2l0 = cos(2.0 * degToRad(l0));
double sin4l0 = sin(4.0 * degToRad(l0));
double sin2m = sin(2.0 * degToRad(m));
double Etime = y * sin2l0 - 2.0 * e * sinm + 4.0 * e * y * sinm * cos2l0 - 0.5 * y * y * sin4l0 - 1.25 * e * e * sin2m;
return radToDeg(Etime) * 4.0; // in minutes of time
}
double gSystem::calcTimeJulianCent(double jd) {
double T = (jd - 2451545.0) / 36525.0;
return T;
}
double gSystem::calcSunTrueLong(double t) {
double l0 = calcGeomMeanLongSun(t);
double c = calcSunEqOfCenter(t);
double O = l0 + c;
return O; // in degrees
}
double gSystem::calcSunApparentLong(double t) {
double o = calcSunTrueLong(t);
double omega = 125.04 - 1934.136 * t;
double lambda = o - 0.00569 - 0.00478 * sin(degToRad(omega));
return lambda; // in degrees
}
double gSystem::calcSunDeclination(double t) {
double e = calcObliquityCorrection(t);
double lambda = calcSunApparentLong(t);
double sint = sin(degToRad(e)) * sin(degToRad(lambda));
double theta = radToDeg(asin(sint));
return theta; // in degrees
}
double gSystem::calcHourAngleSunrise(double lat, double solarDec) {
double latRad = degToRad(lat);
double sdRad = degToRad(solarDec);
double HA = (acos(cos(degToRad(90.833)) / (cos(latRad) * cos(sdRad)) - tan(latRad) * tan(sdRad)));
return HA; // in radians
}
double gSystem::calcHourAngleSunset(double lat, double solarDec) {
double latRad = degToRad(lat);
double sdRad = degToRad(solarDec);
double HA = (acos(cos(degToRad(90.833)) / (cos(latRad) * cos(sdRad)) - tan(latRad) * tan(sdRad)));
return -HA; // in radians
}
double gSystem::calcJD(int year, int month, int day) {
if (month <= 2) {
year -= 1;
month += 12;
}
int A = floor(year / 100);
int B = 2 - A + floor(A / 4);
double JD = floor(365.25 * (year + 4716)) + floor(30.6001 * (month + 1)) + day + B - 1524.5;
return JD;
}
double gSystem::calcJDFromJulianCent(double t) {
double JD = t * 36525.0 + 2451545.0;
return JD;
}
double gSystem::calcSunEqOfCenter(double t) {
double m = calcGeomMeanAnomalySun(t);
double mrad = degToRad(m);
double sinm = sin(mrad);
double sin2m = sin(mrad + mrad);
double sin3m = sin(mrad + mrad + mrad);
double C = sinm * (1.914602 - t * (0.004817 + 0.000014 * t)) + sin2m * (0.019993 - 0.000101 * t) + sin3m * 0.000289;
return C; // in degrees
}
double gSystem::calcSunriseUTC(double JD, double latitude, double longitude) {
double t = calcTimeJulianCent(JD);
double eqTime = calcEquationOfTime(t);
double solarDec = calcSunDeclination(t);
double hourAngle = calcHourAngleSunrise(latitude, solarDec);
double delta = longitude - radToDeg(hourAngle);
double timeDiff = 4 * delta; // in minutes of time
double timeUTC = 720 + timeDiff - eqTime; // in minutes
double newt = calcTimeJulianCent(calcJDFromJulianCent(t) + timeUTC / 1440.0);
eqTime = calcEquationOfTime(newt);
solarDec = calcSunDeclination(newt);
hourAngle = calcHourAngleSunrise(latitude, solarDec);
delta = longitude - radToDeg(hourAngle);
timeDiff = 4 * delta;
timeUTC = 720 + timeDiff - eqTime;
return timeUTC;
}
double gSystem::calcSunsetUTC(double JD, double latitude, double longitude) {
double t = calcTimeJulianCent(JD);
double eqTime = calcEquationOfTime(t);
double solarDec = calcSunDeclination(t);
double hourAngle = calcHourAngleSunset(latitude, solarDec);
double delta = longitude - radToDeg(hourAngle);
double timeDiff = 4 * delta; // in minutes of time
double timeUTC = 720 + timeDiff - eqTime; // in minutes
double newt = calcTimeJulianCent(calcJDFromJulianCent(t) + timeUTC / 1440.0);
eqTime = calcEquationOfTime(newt);
solarDec = calcSunDeclination(newt);
hourAngle = calcHourAngleSunset(latitude, solarDec);
delta = longitude - radToDeg(hourAngle);
timeDiff = 4 * delta;
timeUTC = 720 + timeDiff - eqTime;
return timeUTC;
}
This snippet calculating the following wrongly:
**** START CALCULATE SUN TIMES ****
GS - JD: 2459676.500000
GS - Time: 2022-4-7 9:38:47
GS - Latitude: 47.943272
GS - Longitude: 22.316900
**** Sun times calculated ****
GS - Sun Rise Info: 04:36
GS - Sun Set Info: 17:47
According to google for this latitude and longitude, the info should be:
**** Sun times FROM GOOGLE ****
Sun Rise Info: 05:55
Sun Set Info: 19:09
What do i do wrong? The logic is working fine but it has a big difference.
