long latitude       = 47.94327164;
long longitude      = 22.31690025;
long sunSetHour     = 0;
long sunSetMinute   = 0;
long sunRiseHour    = 0;
long sunRiseMinute  = 0;
boolean dst         = false;
boolean gotGeoData  = true;
#define SUN_TIMES_DEBUG true

struct t_Config {
  int year  = 2022;
  int month = 8;
  int day   = 23;
  int hour  = 16;
  int min   = 6;
  int sec   = 30;
};
t_Config timeConf;

long lastTimeUpdate   = 0;
long lastSunTimeCalc  = 0;

void setup() {
  Serial.begin(115200);
  updateTime();
  calcSunTimes(timeConf.year,timeConf.month,timeConf.day,timeConf.hour,timeConf.min,timeConf.sec);
}

void loop() {
  updateTime();
  if( millis() - lastSunTimeCalc >= 5000 ){
    lastSunTimeCalc = millis();
    calcSunTimes(timeConf.year,timeConf.month,timeConf.day,timeConf.hour,timeConf.min,timeConf.sec);
  }
}


void calcSunTimes( int year, int month, int day, int hour, int min, int sec ) {
  // First we need to check if we already got latitude data.

  if ( !gotGeoData ) {
    return;
  }

  time_t seconds;
  time_t tseconds;
  struct tm *ptm = NULL;
  struct tm tm;

  float JD = calcJD( year, month, day );

  tm.tm_year = year -1900;
  tm.tm_mon  = month -1;
  tm.tm_mday = day;
  tm.tm_hour = hour;
  tm.tm_min  = min;
  tm.tm_sec  = sec;

  seconds = mktime(&tm);
  int delta;
  ptm = gmtime(&seconds);
  delta = ptm->tm_hour;

  // Is winter time.
  int offsetToAdd = 1;
  if ( dst ) {
    offsetToAdd = 0;
  }

  // Calc sunRise
  tseconds = seconds;
  seconds = seconds + calcSunriseUTC(JD, latitude, -longitude) * 60;
  seconds = seconds - delta * 3600;
  ptm = gmtime(&seconds);

  sunRiseHour   = ptm->tm_hour + offsetToAdd;
  sunRiseMinute = ptm->tm_min;

  // Calc sunSet
  seconds = tseconds;
  seconds += calcSunsetUTC(JD, latitude, -longitude) * 60;
  seconds = seconds - delta * 3600;

  ptm = gmtime(&seconds);
  sunSetHour   = ptm->tm_hour + offsetToAdd;
  sunSetMinute = ptm->tm_min;

#if SUN_TIMES_DEBUG
  Serial.println("*********** ***********");
  Serial.printf("[SUN_TIMES] - Millis: %lu\n", millis() );
  Serial.printf("[SUN_TIMES] - Time: %d - %02d/%02d %02d:%02d:%02d\n",
      year,month,day,hour,min,sec);
  Serial.printf("[SUN_TIMES] - SunRise: %d:%d\n", sunRiseHour, sunRiseMinute);
  Serial.printf("[SUN_TIMES] - SunSet: %d:%d\n", sunSetHour, sunSetMinute);
  Serial.println("*********** ***********");
#endif

}

void updateTime() {
  if( millis() - lastTimeUpdate >= 1000 ){
    lastTimeUpdate = millis();
    timeConf.sec++;
    if ( timeConf.sec > 60 ) {
      timeConf.sec = 0;
      timeConf.min++;
      if ( timeConf.min > 60 ) {
        timeConf.min = 0;
        timeConf.hour++;
        if (timeConf.hour > 24 ) {
          timeConf.hour = 0;
          timeConf.day++;
        }
      }
    }
  }
}

/**************************  CALCULATIONS BELOW **************************/

/* Convert degree angle to radians */
double degToRad(double angleDeg) {
  return (PI * angleDeg / 180.0);
}

double radToDeg(double angleRad) {
  return (180.0 * angleRad / PI);
}

double 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 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 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 calcEccentricityEarthOrbit(double t) {
  double e = 0.016708634 - t * (0.000042037 + 0.0000001267 * t);
  return e;  // unitless
}

double calcGeomMeanAnomalySun(double t) {
  double M = 357.52911 + t * (35999.05029 - 0.0001537 * t);
  return M;  // in degrees
}

double 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 calcTimeJulianCent(double jd) {
  double T = (jd - 2451545.0) / 36525.0;
  return T;
}

double calcSunTrueLong(double t) {
  double l0 = calcGeomMeanLongSun(t);
  double c = calcSunEqOfCenter(t);
  double O = l0 + c;
  return O;  // in degrees
}

double 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 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 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 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 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 calcJDFromJulianCent(double t) {
  double JD = t * 36525.0 + 2451545.0;
  return JD;
}

double 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 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 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;
}