/* https://www.esrl.noaa.gov/gmd/grad/solcalc/
   https://www.esrl.noaa.gov/gmd/grad/solcalc/azel.html

  The calculations in the NOAA Sunrise/Sunset and Solar Position Calculators are based on equations from Astronomical Algorithms, by Jean Meeus.
  The sunrise and sunset results are theoretically accurate to within a minute for locations between +/- 72° latitude,
  and within 10 minutes outside of those latitudes.
  However, due to variations in atmospheric composition, temperature, pressure, and conditions, observed values may vary from calculations

  The approximations used in these programs are very good for years between 1800 and 2100.
  Results should still be sufficiently accurate for the range from -1000 to 3000.
  Outside of this range, results may be given, but the potential for error is higher.
*/

/*
  The solar altitude is the angle between the direction of the Sun and the horizontal plane.
  It ranges from 0° above the horizon to 90° at the zenith.

  The solar azimuth is the angle created between the vertical plane passing through both the Sun and the observer,
  and the North-South vertical plane. This angle is 0° at North, 90° at East, 180° at South, and 270° at West.
*/

/* Formulas extracted from Excel files at https://www.esrl.noaa.gov/gmd/grad/solcalc/calcdetails.html */

/* Verification with https://www.suncalc.org/ - of course double on 4 bytes creates approximation issues on UNO */

// Decimal Degrees
double latitude = 46.174838; // MAISON ST DIDIER
double longitude = 4.827561; // MAISON ST DIDIER

// ***********************************************************************
// IDEALLY SET THE RTC TO UTC TIME AND SET deltaTimeZone TO 0
// This way, you don't need to change the program for daylight saving time.
// ***********************************************************************
int8_t deltaTimeZone = +2;  // FRANCE UTC+1 during standard time and UTC+2 during daylight saving time

#include <RTClib.h>

RTC_DS1307 rtc;
const char* dayNames[] = {"Dimanche", "Lundi", "Mardi", "Mercredi", "Jeudi", "Vendredi", "Samedi"};

void fatalError(const __FlashStringHelper *message)
{
  Serial.print(F("ERROR: ")); Serial.println(message);
  while (true); // We halt here
}

void printWithTwoDigits(uint32_t number)
{
  if (number < 10) Serial.write('0');
  Serial.print(number);
}

void printFractionalTimeOfDay(double t)
{
  uint32_t secondsSinceMidnight = t * (3600ul * 24ul);
  Serial.print(secondsSinceMidnight / 3600UL);
  Serial.write(':');
  printWithTwoDigits((secondsSinceMidnight % 3600UL) / 60);
}


void printMinuteDuration(double nbMinutes)
{
  byte h = nbMinutes / 60;
  byte m = (nbMinutes / 60 - h) * 60;
  byte s = (((nbMinutes / 60 - h) * 60) - m) * 60;

  Serial.print(h);
  Serial.write(':');
  printWithTwoDigits(m);
  // the seconds but we are not that precise anyway
  // Serial.write(':');
  // printWithTwoDigits(s);
}

inline double RADIANS(double a)
{
  return a * (M_PI / 180.0);
}

inline double DEGREES(double a)
{
  return a * (180.0 / M_PI);
}

void calculateAstronomy(
  double latitude,         /* Latitude  (+ to N) */
  double longitude,        /* Longitude (+ to E) */
  int8_t deltaTimeZone,    /* Time Zone (+ to E) */
  uint32_t numberOfDays,  /* Number of days since January 1, 1900 */
  uint32_t secondsElapsed /* Number of seconds since midnight */
)
{
  double localHour = ((double) secondsElapsed) / (24.0 * 3600.0); // Fraction of a day
  double julianDay = 2415018.5 + (double) numberOfDays + (double) localHour - ((double) deltaTimeZone) / 24.0;
  double julianCentury = (julianDay - 2451545L) / 36525L;
  double meanSolarLongitude = fmod(280.46646 + julianCentury * (36000.76983 + julianCentury * 0.0003032), 360);
  double meanSolarAnomaly = 357.52911 + julianCentury * (35999.05029 - 0.0001537 * julianCentury);
  double eccentricityEarthOrbit = 0.016708634 - julianCentury * (0.000042037 + 0.0000001267 * julianCentury);
  double sunEquationOfCenter = sin(RADIANS(meanSolarAnomaly)) * (1.914602 - julianCentury * (0.004817 + 0.000014 * julianCentury)) + sin(RADIANS(2 * meanSolarAnomaly)) * (0.019993 - 0.000101 * julianCentury) + sin(RADIANS(3 * meanSolarAnomaly)) * 0.000289;
  double trueSolarLongitude = meanSolarLongitude + sunEquationOfCenter;
  double apparentSolarLongitude = trueSolarLongitude - 0.00569 - 0.00478 * sin(RADIANS(125.04 - 1934.136 * julianCentury));
  double meanObliquityEcliptic = 23 + (26 + ((21.448 - julianCentury * (46.815 + julianCentury * (0.00059 - julianCentury * 0.001813)))) / 60) / 60;
  double obliquityCorrection = meanObliquityEcliptic + 0.00256 * cos(RADIANS(125.04 - 1934.136 * julianCentury));
  double solarDeclination = DEGREES(asin(sin(RADIANS(obliquityCorrection)) * sin(RADIANS(apparentSolarLongitude))));
  double varY = tan(RADIANS(obliquityCorrection / 2)) * tan(RADIANS(obliquityCorrection / 2));
  double equationOfTime = 4 * DEGREES(varY * sin(2 * RADIANS(meanSolarLongitude)) - 2 * eccentricityEarthOrbit * sin(RADIANS(meanSolarAnomaly)) + 4 * eccentricityEarthOrbit * varY * sin(RADIANS(meanSolarAnomaly)) * cos(2 * RADIANS(meanSolarLongitude)) - 0.5 * varY * varY * sin(4 * RADIANS(meanSolarLongitude)) - 1.25 * eccentricityEarthOrbit * eccentricityEarthOrbit * sin(2 * RADIANS(meanSolarAnomaly)));
  double trueSolarTime = fmod(localHour * 1440 + equationOfTime + 4 * longitude - 60 * deltaTimeZone, 1440);
  double hourAngle = (trueSolarTime / 4 < 0) ? (trueSolarTime / 4 + 180) : (trueSolarTime / 4 - 180);
  double solarZenithAngle = DEGREES(acos(sin(RADIANS(latitude)) * sin(RADIANS(solarDeclination)) + cos(RADIANS(latitude)) * cos(RADIANS(solarDeclination)) * cos(RADIANS(hourAngle))));
  double solarElevationAngle = 90 - solarZenithAngle;

  double approximateAtmosphericRefractionAngle ;
  if (solarElevationAngle > 85) approximateAtmosphericRefractionAngle = 0;
  else if (solarElevationAngle > 5)
    approximateAtmosphericRefractionAngle = 58.1 / tan(RADIANS(solarElevationAngle)) - 0.07 / pow(tan(RADIANS(solarElevationAngle)), 3) + 0.000086 / pow(tan(RADIANS(solarElevationAngle)), 5);
  else if (solarElevationAngle > -0.575) approximateAtmosphericRefractionAngle = 1735 + solarElevationAngle * (-518.2 + solarElevationAngle * (103.4 + solarElevationAngle * (-12.79 + solarElevationAngle * 0.711)));
  else approximateAtmosphericRefractionAngle = -20.772 / tan(RADIANS(solarElevationAngle));
  approximateAtmosphericRefractionAngle /= 3600.0;

  double finalSolarElevation = solarElevationAngle + approximateAtmosphericRefractionAngle;
  double finalAzimuthAngle = (hourAngle > 0) ?
                             fmod(DEGREES(acos(((sin(RADIANS(latitude)) * cos(RADIANS(solarZenithAngle))) - sin(RADIANS(solarDeclination))) / (cos(RADIANS(latitude)) * sin(RADIANS(solarZenithAngle))))) + 180, 360) :
                             fmod(540 - DEGREES(acos(((sin(RADIANS(latitude)) * cos(RADIANS(solarZenithAngle))) - sin(RADIANS(solarDeclination))) / (cos(RADIANS(latitude)) * sin(RADIANS(solarZenithAngle))))), 360);



  double sunTrueAnomaly = meanSolarAnomaly + sunEquationOfCenter;
  double sunRadiusVector = (1.000001018 * (1 - eccentricityEarthOrbit * eccentricityEarthOrbit)) / (1 + eccentricityEarthOrbit * cos(RADIANS(sunTrueAnomaly)));
  double sunRightAscension = DEGREES(atan2(cos(RADIANS(apparentSolarLongitude)), cos(RADIANS(obliquityCorrection)) * sin(RADIANS(apparentSolarLongitude))));
  double hourAngleSunrise = DEGREES(acos(cos(RADIANS(90.833)) / (cos(RADIANS(latitude)) * cos(RADIANS(solarDeclination))) - tan(RADIANS(latitude)) * tan(RADIANS(solarDeclination))));
  double solarNoon = (720 - 4 * longitude - equationOfTime + deltaTimeZone * 60) / 1440;
  double sunriseTime = solarNoon - hourAngleSunrise * 4 / 1440;
  double sunsetTime = solarNoon + hourAngleSunrise * 4 / 1440;
  double daylightDuration = 8 * hourAngleSunrise;

  Serial.print(secondsElapsed / 3600UL);
  Serial.write(':');
  printWithTwoDigits((secondsElapsed % 3600UL) / 60);

  Serial.write('\t');
  Serial.print(finalAzimuthAngle, 1);

  Serial.print("°\t");
  Serial.print(finalSolarElevation, 1);

  Serial.write('°\t');
  printFractionalTimeOfDay(sunriseTime);

  Serial.write('\t');
  printFractionalTimeOfDay(sunsetTime);

  Serial.write('\t');
  printMinuteDuration(daylightDuration);

  Serial.print("\n");
}

void setup() {
  uint32_t numberOfDays;        // Number of days since January 1, 1900
  Serial.begin(115200);

  if (! rtc.begin()) fatalError(F("RTC Clock not found"));
  rtc.adjust(DateTime(F(__DATE__), F(__TIME__))); // Do this ONCE to set the RTC to the current time if your RTC is not already set, then comment it out.
  // rtc.adjust(DateTime(2020, 1, 2, 10, 30, 0)); // Or use this to set a specific date and time, e.g., January 2nd at 10:30 AM
  // rtc.adjust(DateTime(2020, 1, 3, 21, 45, 0)); // Or use this to set a specific date and time, e.g., January 3rd at 9:45 PM

  DateTime now = rtc.now(); // Current time
  DateTime firstJanuary2020(2020, 1, 1, 0, 0, 0); // Reference point: January 1, 2020, at 0:00 UTC
  TimeSpan sinceJanuary2020 = now - firstJanuary2020 ; // Time difference between now and January 1, 1900
  numberOfDays = 43829UL + sinceJanuary2020.days() ; // + 43829 days to count from January 1, 1900

  Serial.println(F("\nSun information for : "));
  Serial.print(dayNames[now.dayOfTheWeek()]);
  Serial.write(' ');
  Serial.print(now.day());
  Serial.write('/');
  Serial.print(now.month());
  Serial.write('/');
  Serial.print(now.year());
  Serial.write('\t');
  Serial.print(now.hour());
  Serial.write(':');
  printWithTwoDigits(now.minute());
  Serial.write(':');
  printWithTwoDigits(now.second());
  Serial.write('\n');

  Serial.print(F("Latitude : ")); Serial.println(latitude, 6);
  Serial.print(F("Longitude : ")); Serial.println(longitude, 6);

  Serial.println("\nHeure\tAzim\tElev\tLevé\tCouché\tDuréeJour");

  for (uint32_t moment = 22500; moment <= 79200UL; moment += 900UL) { // from 6:15 AM to 10:00 PM
    calculateAstronomy(latitude, longitude, deltaTimeZone, numberOfDays, moment); // every 15 minutes
  }

}

void loop() {}