Skip to content

Commit

Permalink
Update SolarCalculator.h/.cpp
Browse files Browse the repository at this point in the history
  • Loading branch information
@jpb10
jpb10 committed Jan 16, 2023
1 parent 45fa4e3 commit 02e26c2
Show file tree
Hide file tree
Showing 2 changed files with 47 additions and 117 deletions.
146 changes: 42 additions & 104 deletions src/SolarCalculator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 <Arduino.h>
#else
#include <cmath>
#endif

Expand All @@ -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<int>(7 * (Y + (M + 9) / 12) / 4) + static_cast<int>(275 * M / 9) + D + 1721013.5;
m = (hour + minute / 60.0 + second / 3600.0) / 24.0;
}

//======================================================================================================================
Expand Down Expand Up @@ -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<int>(7 * (year + (month + 9) / 12) / 4) + static_cast<int>(275 * month / 9) +
day + 1721013.5;
}

double calcJulianCent(JulianDay jd)
{
return (jd.JD - 2451545 + jd.m) / 36525;
Expand Down Expand Up @@ -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)
}

//======================================================================================================================
Expand All @@ -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
Expand All @@ -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;
Expand Down
18 changes: 5 additions & 13 deletions src/SolarCalculator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 <Arduino.h>
#endif

//namespace solarcalculator {

constexpr double SUNRISESET_STD_ALTITUDE = -0.8333;
Expand All @@ -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
Expand All @@ -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
Expand Down

0 comments on commit 02e26c2

Please sign in to comment.