mup_astro.h

#pragma once
 
#include <math.h>
 
namespace ustd {
 
const double C_PI = 3.1415926535897932384626433;  
const double C_D2R = C_PI / 180.0;                
const double C_R2D = 180.0 / C_PI;                
const double C_AU = 149597870700.0;               
const double C_C = 299792458.0;                   
const double C_CAUD = C_C * 60 * 60 * 24 / C_AU;  
const double C_MJD = 2400000.5;                   
 
class Astro {
  public:
    double lat, lon, utcOffset;
 
#ifdef USTD_FEATURE_FILESYSTEM
    Astro() {
    }
#endif
    Astro(double lat, double lon, double utcOffset)
        : lat(lat), lon(lon), utcOffset(utcOffset) {
    }
 
    static long julianDayNumber(int year, uint8_t month, uint8_t day) {
        long JDN = (1461L * ((long)year + 4800L + ((long)month - 14L) / 12L)) / 4L +
                   (367L * ((long)month - 2L - 12L * (((long)month - 14L) / 12L))) / 12L -
                   (3L * (((long)year + 4900L + ((long)month - 14L) / 12L) / 100L)) / 4L +
                   (long)day - 32075L;
        return JDN;
    }
 
    static double fracDay(uint8_t hour, uint8_t min, double sec) {
        double dayFrac = ((double)hour + (double)min / 60.0 + (double)sec / 3600.0) / 24.0;
        return dayFrac;
    }
 
    static double julianDate(int year, uint8_t month, uint8_t day, uint8_t hour, uint8_t min,
                             double sec) {
        long JDN = julianDayNumber(year, month, day);
        double dfrac = fracDay(hour, min, sec) - 0.5;
        double JD = (double)JDN + dfrac;
        return JD;
    }
 
    static double modifiedJulianDate(int year, uint8_t month, uint8_t day, uint8_t hour,
                                     uint8_t min, double sec) {
        long JDN = julianDayNumber(year, month, day);
        double dfrac = fracDay(hour, min, sec) - 0.5;
        double JD = (double)(JDN - C_MJD) + dfrac;
        return JD;
    }
 
    static bool calculateSunRiseSet(int year, int month, int day, double lat, double lon,
                                    int localOffset, int daylightSavings, bool bRising,
                                    double *pSunTime) {
        // 1. first calculate the day of the year
        const double ZENITH = 90.0 + 50.0 / 60.0;
        double N1 = floor(275.0 * month / 9.0);
        double N2 = floor((month + 9.0) / 12.0);
        double N3 = (1.0 + floor((year - 4 * floor(year / 4.0) + 2.0) / 3.0));
        double N = N1 - (N2 * N3) + day - 30.0;
 
        // 2. convert the longitude to hour value and calculate an approximate time
        double lonHour = lon / 15.0;
        double t;
        if (bRising)
            t = N + ((6 - lonHour) / 24);  // if rising time is desired:
        else
            t = N + ((18 - lonHour) / 24);  // if setting time is desired:
 
        // 3. calculate the Sun's mean anomaly
        double M = (0.9856 * t) - 3.289;
 
        // 4. calculate the Sun's true longitude
        double L =
            fmod(M + (1.916 * sin(C_D2R * M)) + (0.020 * sin(2 * C_D2R * M)) + 282.634, 360.0);
 
        // 5a. calculate the Sun's right ascension
        double RA = fmod(C_R2D * atan(0.91764 * tan(C_D2R * L)), 360.0);
 
        // 5b. right ascension value needs to be in the same quadrant as L
        double Lquadrant = floor(L / 90.0) * 90.0;
        double RAquadrant = floor(RA / 90.0) * 90.0;
        RA = RA + (Lquadrant - RAquadrant);
 
        // 5c. right ascension value needs to be converted into hours
        RA = RA / 15.0;
 
        // 6. calculate the Sun's declination
        double sinDec = 0.39782 * sin(C_D2R * L);
        double cosDec = cos(asin(sinDec));
 
        // 7a. calculate the Sun's local hour angle
        // double cosH = (sin(C_D2R * ZENITH) - (sinDec * sin(C_D2R * lat))) /
        //              (cosDec * cos(C_D2R * lat));
        double cosH =
            (cos(C_D2R * ZENITH) - (sinDec * sin(lat * C_D2R))) / (cosDec * cos(C_D2R * lat));
 
        if (bRising) {
            if (cosH > 1.0) {  // the sun never rises on this location (on the specified date)
                *pSunTime = -1.0;
                return false;
            }
        } else {
            if (cosH < -1.0) {  // the sun never sets on this location (on the specified date)
                *pSunTime = -1.0;
                return false;
            }
        }
 
        // 7b. finish calculating H and convert into hours
        double H;
        if (bRising)
            H = 360.0 - C_R2D * acos(cosH);  //   if if rising time is desired:
        else
            H = C_R2D * acos(cosH);  //   if setting time is desired:
        H = H / 15.0;
 
        // 8. calculate local mean time of rising/setting
        double T = H + RA - (0.06571 * t) - 6.622;
 
        // 9. adjust back to UTC
        double UT = fmod(T - lonHour, 24.0);
 
        // 10. convert UT value to local time zone of latitude/longitude
        *pSunTime = UT + localOffset + daylightSavings;
        return true;
    }
 
    static int cmpHourMinuteTime(uint8_t h1, uint8_t m1, uint8_t h2, uint8_t m2) {
        if (h2 > h1)
            return 1;
        if (h2 < h1)
            return -1;
        if (m2 > m1)
            return 1;
        if (m2 < m1)
            return -1;
        return 0;
    }
 
    static int deltaHourMinuteTime(uint8_t h1, uint8_t m1, uint8_t h2, uint8_t m2) {
        if (m2 < m1) {
            m2 += 60;
            if (h2 > 0)
                --h2;
            else
                h2 = 23;
        }
        if (h2 < h1)
            h2 += 24;
        return (h2 - h1) * 60 + m2 - m1;
    }
 
    static bool inHourMinuteInterval(uint8_t test_hour, uint8_t test_minute, uint8_t start_hour,
                                     uint8_t start_minute, uint8_t end_hour, uint8_t end_minute) {
        if (cmpHourMinuteTime(start_hour, start_minute, end_hour, end_minute) > 0) {
            if (cmpHourMinuteTime(test_hour, test_minute, start_hour, start_minute) <= 0 &&
                cmpHourMinuteTime(test_hour, test_minute, end_hour, end_minute) >= 0)
                return true;
            else
                return false;
        } else {
            if (cmpHourMinuteTime(test_hour, test_minute, start_hour, start_minute) > 0 &&
                cmpHourMinuteTime(test_hour, test_minute, end_hour, end_minute) < 0)
                return false;
            else
                return true;
        }
        return false;
    }
 
    static bool parseHourMinuteString(String hourMinute, int *hour, int *minute) {
        int ind = hourMinute.indexOf(':');
        if (ind == -1)
            return false;
        int hr = hourMinute.substring(0, ind).toInt();
        int mn = hourMinute.substring(ind + 1).toInt();
        if (hr < 0 || hr > 23)
            return false;
        if (mn < 0 || mn > 59)
            return false;
        *hour = hr;
        *minute = mn;
        return true;
    }
};
 
}  // namespace ustd