From 02e26c2caa6efb60de284046761fde422f43603f Mon Sep 17 00:00:00 2001 From: jpb10 <4648556+jpb10@users.noreply.github.com> Date: Mon, 16 Jan 2023 13:57:44 -0500 Subject: [PATCH] Update SolarCalculator.h/.cpp --- src/SolarCalculator.cpp | 146 ++++++++++++---------------------------- src/SolarCalculator.h | 18 ++--- 2 files changed, 47 insertions(+), 117 deletions(-) diff --git a/src/SolarCalculator.cpp b/src/SolarCalculator.cpp index 264e443..2e002d5 100644 --- a/src/SolarCalculator.cpp +++ b/src/SolarCalculator.cpp @@ -4,10 +4,12 @@ // This library provides functions to calculate the times of sunrise, sunset, solar noon, twilight (dawn and dusk), // Sun's apparent position in the sky, equation of time, etc. // -// Most formulae are taken from Astronomical Algorithms by Jean Meeus and optimized for 8-bit AVR platform. +// Most formulae are taken from Astronomical Algorithms by Jean Meeus and optimized for Arduino. //====================================================================================================================== -#ifndef ARDUINO +#ifdef ARDUINO +#include +#else #include #endif @@ -21,10 +23,11 @@ JulianDay::JulianDay(unsigned long utc) m = (utc % 86400) / 86400.0; } -JulianDay::JulianDay(int year, int month, int day, int hour, int minute, int second) +// Valid from 1901 to 2099, Van Flandern & Pulkkinen (1979) +JulianDay::JulianDay(int Y, int M, int D, int hour, int minute, int second) { - JD = calcJulianDay(year, month, day); - m = fractionalDay(hour, minute, second); + JD = 367.0 * Y - static_cast(7 * (Y + (M + 9) / 12) / 4) + static_cast(275 * M / 9) + D + 1721013.5; + m = (hour + minute / 60.0 + second / 3600.0) / 24.0; } //====================================================================================================================== @@ -58,38 +61,6 @@ double wrapTo180(double angle) return angle - 180; // [-180, 180) } -// Interpolation of three tabular values, valid for n between -1 and +1 -double interpolateCoordinates(double n, double y1, double y2, double y3) -{ - if (fabs(y2 - y1) > 180) // if coordinate is discontinuous - { // add or subtract 360 degrees - if (y1 < 0) y1 += 360; - else if (y2 < 0) y1 -= 360; - } - else if (fabs(y3 - y2) > 180) - { - if (y3 < 0) y3 += 360; - else if (y2 < 0) y3 -= 360; - } - - double a = y2 - y1; - double b = y3 - y2; - double c = b - a; - return y2 + n * (a + b + n * c) / 2; -} - -double fractionalDay(int hour, int minute, int second) -{ - return (hour + minute / 60.0 + second / 3600.0) / 24; -} - -// Valid from 1901 to 2099, Van Flandern & Pulkkinen (1979) -double calcJulianDay(int year, int month, int day) -{ - return 367.0 * year - static_cast(7 * (year + (month + 9) / 12) / 4) + static_cast(275 * month / 9) + - day + 1721013.5; -} - double calcJulianCent(JulianDay jd) { return (jd.JD - 2451545 + jd.m) / 36525; @@ -136,52 +107,30 @@ void calcSolarCoordinates(double T, double &ra, double &dec) double calcGrMeanSiderealTime(JulianDay jd) { - double GMST = wrapTo360(100.46061837 + 0.98564736629 * (jd.JD - 2451545)); - return wrapTo360(GMST + 360.985647 * jd.m); // in degrees + double GMST0 = wrapTo360(100.46061837 + 0.98564736629 * (jd.JD - 2451545)); + return wrapTo360(GMST0 + 360.985647 * jd.m); // in degrees } void equatorial2horizontal(double H, double dec, double lat, double &az, double &el) { - az = degrees(atan2(sin(radians(H)), cos(radians(H)) * sin(radians(lat)) - - tan(radians(dec)) * cos(radians(lat)))); - el = degrees(asin(sin(radians(lat)) * sin(radians(dec)) + - cos(radians(lat)) * cos(radians(dec)) * cos(radians(H)))); -} - -// Approximate atmospheric refraction correction, in degrees -double calcRefraction(double elev) -{ - if (elev < -0.575) - return -20.774 / tan(radians(elev)) / 3600; // Zimmerman (1981) - else - return 1.02 / tan(radians(elev + 10.3 / (elev + 5.11))) / 60; // Sæmundsson (1986) + az = degrees(atan2(sin(radians(H)), cos(radians(H)) * sin(radians(lat)) - tan(radians(dec)) * cos(radians(lat)))); + el = degrees(asin(sin(radians(lat)) * sin(radians(dec)) + cos(radians(lat)) * cos(radians(dec)) * cos(radians(H)))); } -// Equation of ephemeris time by Smart (1978) -double equationOfTimeSmart(double T) +// Hour angle at sunrise or sunset, returns NaN if circumpolar +double calcHourAngleRiseSet(double dec, double lat, double h0) { - double L0 = calcGeomMeanLongSun(T); - double M = calcGeomMeanAnomalySun(T); - - // Using numerical values for the eccentricity and obliquity in 2025 - return 2.465 * sin(2 * radians(L0)) - 1.913 * sin(radians(M)) + - 0.165 * sin(radians(M)) * cos(2 * radians(L0)) - - 0.053 * sin(4 * radians(L0)) - 0.02 * sin(2 * radians(M)); // in degrees + return degrees(acos((sin(radians(h0)) - sin(radians(lat)) * sin(radians(dec))) / + (cos(radians(lat)) * cos(radians(dec))))); } -// Simple polynomial expressions for delta T (ΔT), in seconds of time -double calcDeltaT(double year) +// Approximate atmospheric refraction correction, in degrees +double calcRefraction(double el) { - if (year > 1997) - { - double t = year - 2015; - return 67.62 + t * (0.3645 + 0.0039755 * t); // Fred Espenak (2014) - } - else // y > 948 - { - double u = (year - 2000) / 100; - return 64.69 + u * (80.59 + 23.604 * u); // fitted to historical data, very approximate before 1900 - } + if (el < -0.575) + return -20.774 / tan(radians(el)) / 3600; // Zimmerman (1981) + else + return 1.02 / tan(radians(el + 10.3 / (el + 5.11))) / 60; // Sæmundsson (1986) } //====================================================================================================================== @@ -194,7 +143,12 @@ double calcDeltaT(double year) void calcEquationOfTime(JulianDay jd, double &E) { double T = calcJulianCent(jd); - E = 4 * equationOfTimeSmart(T); + double L0 = calcGeomMeanLongSun(T); + + double ra, dec; + calcSolarCoordinates(T, ra, dec); + + E = 4 * wrapTo180(L0 - 0.00569 - ra); } // Sun's geocentric (as seen from the center of the Earth) equatorial coordinates, in degrees and AUs @@ -220,49 +174,33 @@ void calcHorizontalCoordinates(JulianDay jd, double latitude, double longitude, equatorial2horizontal(H, dec, latitude, azimuth, elevation); azimuth += 180; // measured from the North -// elevation -= 8.794 * cos(radians(elevation)) / 3600; // parallax in altitude, always < 0.0025 degrees elevation += calcRefraction(elevation); } -// Helper function -void calcRiseSetTimes(double (&m)[3], JulianDay jd, double latitude, double longitude, double h0) -{ - double T = calcJulianCent(jd); - double GMST = calcGrMeanSiderealTime(jd); - - double ra, dec; - calcSolarCoordinates(T, ra, dec); - - // Local hour angle at sunrise or sunset (±NaN if body is circumpolar) - double H0 = degrees(acos((sin(radians(h0)) - sin(radians(latitude)) * sin(radians(dec))) / - (cos(radians(latitude)) * cos(radians(dec))))); - - m[0] = jd.m + wrapTo180(ra - longitude - GMST) / 360; - m[1] = m[0] - H0 / 360; - m[2] = m[0] + H0 / 360; -} - // Find the times of sunrise, transit, and sunset, in hours void calcSunriseSunset(JulianDay jd, double latitude, double longitude, double &transit, double &sunrise, double &sunset, double altitude, int iterations) { - double m[3], times[3]; + double m[3]; m[0] = 0.5 - longitude / 360; for (int i = 0; i <= iterations; i++) for (int event = 0; event < 3; event++) { jd.m = m[event]; - calcRiseSetTimes(times, jd, latitude, longitude, altitude); - m[event] = times[event]; - - // First iteration, approximate rise and set times - if (i == 0) - { - m[1] = times[1]; - m[2] = times[2]; - break; - } + double T = calcJulianCent(jd); + double GMST = calcGrMeanSiderealTime(jd); + + double ra, dec; + calcSolarCoordinates(T, ra, dec); + + double m0 = jd.m + wrapTo180(ra - longitude - GMST) / 360; + double d0 = calcHourAngleRiseSet(dec, latitude, altitude) / 360; + + if (event == 0) m[0] = m0; + if (event == 1 || i == 0) m[1] = m0 - d0; + if (event == 2 || i == 0) m[2] = m0 + d0; + if (i == 0) break; } transit = m[0] * 24; diff --git a/src/SolarCalculator.h b/src/SolarCalculator.h index b501b1a..0ca039f 100644 --- a/src/SolarCalculator.h +++ b/src/SolarCalculator.h @@ -4,16 +4,12 @@ // This library provides functions to calculate the times of sunrise, sunset, solar noon, twilight (dawn and dusk), // Sun's apparent position in the sky, equation of time, etc. // -// Most formulae are taken from Astronomical Algorithms by Jean Meeus and optimized for 8-bit AVR platform. +// Most formulae are taken from Astronomical Algorithms by Jean Meeus and optimized for Arduino. //====================================================================================================================== #ifndef SOLARCALCULATOR_H #define SOLARCALCULATOR_H -#ifdef ARDUINO -#include -#endif - //namespace solarcalculator { constexpr double SUNRISESET_STD_ALTITUDE = -0.8333; @@ -39,11 +35,8 @@ struct JulianDay // Utilities double wrapTo360(double angle); double wrapTo180(double angle); -double interpolateCoordinates(double n, double y1, double y2, double y3); -// Julian day -double fractionalDay(int hour, int minute, int second); -double calcJulianDay(int year, int month, int day); +// Julian centuries double calcJulianCent(JulianDay jd); // Solar coordinates @@ -54,14 +47,13 @@ double calcSunRadVector(double T); double calcMeanObliquityOfEcliptic(double T); void calcSolarCoordinates(double T, double &ra, double &dec); -// Sidereal time at Greenwich, solar time, and ΔT +// Sidereal time at Greenwich double calcGrMeanSiderealTime(JulianDay jd); -double equationOfTimeSmart(double T); -double calcDeltaT(double year); // Sun's position in the sky void equatorial2horizontal(double H, double dec, double lat, double &az, double &el); -double calcRefraction(double elev); +double calcHourAngleRiseSet(double dec, double lat, double h0); +double calcRefraction(double el); //====================================================================================================================== // Solar calculator