19 static const char *
const TAG =
"sun";
26 static const num_t PI = 3.141592653589793;
46 d += moment.
hour / 24.0;
47 d += moment.
minute / (24.0 * 60);
48 d += moment.
second / (24.0 * 60 * 60);
54 int b = 2 - a + a / 4;
55 return ((
int) (365.25 * (y + 4716))) + ((
int) (30.6001 * (m + 1))) + d + b - 1524.5;
59 int t = moment.
year - 2000;
60 return 62.92 + t * (0.32217 + t * 0.005589);
64 num_t res = fmod(x, y);
75 return jd() +
delta_t(dt) / (60 * 60 * 24);
83 SunAtTime(
num_t jde) : jde(jde), t((jde - 2451545) / 36525.0) {}
85 num_t mean_obliquity()
const {
95 num_t omega = 125.05 - 1934.136 * t;
99 num_t true_obliquity()
const {
102 num_t epsilon = mean_obliquity() + delta_epsilon;
106 num_t mean_longitude()
const {
108 num_t l0 = 280.46646 + 36000.76983 * t + 0.0003032 *
sq(t);
109 return wmod(l0, 360);
112 num_t eccentricity()
const {
114 num_t e = 0.016708634 - 0.000042037 * t - 0.0000001267 *
sq(t);
118 num_t mean_anomaly()
const {
120 num_t m = 357.52911 + 35999.05029 * t - 0.0001537 *
sq(t);
124 num_t equation_of_center()
const {
127 num_t c = ((1.914602 - 0.004817 * t - 0.000014 *
sq(t)) * sin(m_rad) + (0.019993 - 0.000101 * t) * sin(2 * m_rad) +
128 0.000289 * sin(3 * m_rad));
132 num_t true_longitude()
const {
134 num_t x = mean_longitude() + equation_of_center();
138 num_t true_anomaly()
const {
140 num_t x = mean_anomaly() + equation_of_center();
144 num_t apparent_longitude()
const {
146 num_t x = true_longitude() - 0.00569 - 0.00478 * sin(
radians(omega()));
154 num_t right_ascension_rad = atan2(cos(epsilon_rad) * sin(app_lon_rad), cos(app_lon_rad));
155 num_t declination_rad = asin(sin(epsilon_rad) * sin(app_lon_rad));
159 num_t equation_of_time()
const {
163 num_t l2 = 2 * mean_longitude();
165 num_t e = eccentricity();
168 num_t sin_m = sin(m_rad);
169 num_t eot = (y * sin(l2_rad) - 2 * e * sin_m + 4 * e * y * sin_m * cos(l2_rad) - 1 / 2.0 *
sq(y) * sin(2 * l2_rad) -
170 5 / 4.0 *
sq(e) * sin(2 * m_rad));
176 ESP_LOGV(TAG,
"jde: %f", jde);
177 ESP_LOGV(TAG,
"T: %f", t);
178 ESP_LOGV(TAG,
"L_0: %f", mean_longitude());
179 ESP_LOGV(TAG,
"M: %f", mean_anomaly());
180 ESP_LOGV(TAG,
"e: %f", eccentricity());
181 ESP_LOGV(TAG,
"C: %f", equation_of_center());
182 ESP_LOGV(TAG,
"Odot: %f", true_longitude());
183 ESP_LOGV(TAG,
"Omega: %f", omega());
184 ESP_LOGV(TAG,
"lambda: %f", apparent_longitude());
185 ESP_LOGV(TAG,
"epsilon_0: %f", mean_obliquity());
186 ESP_LOGV(TAG,
"epsilon: %f", true_obliquity());
187 ESP_LOGV(TAG,
"v: %f", true_anomaly());
188 auto eq = equatorial_coordinate();
189 ESP_LOGV(TAG,
"right_ascension: %f", eq.right_ascension);
190 ESP_LOGV(TAG,
"declination: %f", eq.declination);
194 struct SunAtLocation {
197 num_t greenwich_sidereal_time(
Moment moment)
const {
205 num_t t = (jd0 - 2451545) / 36525;
207 num_t gmst = (+280.46061837 + 360.98564736629 * (jd - 2451545) + 0.000387933 *
sq(t) - (1 / 38710000.0) *
cb(t));
208 return wmod(gmst, 360);
212 auto eq = SunAtTime(moment.
jde()).equatorial_coordinate();
213 num_t gmst = greenwich_sidereal_time(moment);
215 num_t nutation_corr = 0;
217 num_t ra = eq.right_ascension;
219 alpha =
wmod(alpha, 360);
224 num_t sin_elevation = (+sin_lat * sin(eq.declination_rad()) + cos_lat * cos(eq.declination_rad()) * cos(alpha_rad));
225 num_t elevation_rad = asin(sin_elevation);
226 num_t azimuth_rad = atan2(sin(alpha_rad), cos(alpha_rad) * sin_lat - tan(eq.declination_rad()) * cos_lat);
236 auto m = local_event_(date, 12);
238 num_t new_h = 0, old_h;
241 auto x = local_hour_angle_(jde + old_h / 86400, rise, zenith);
245 }
while (std::abs(new_h - old_h) >= 15);
246 time_t new_timestamp =
m.timestamp + (time_t) new_h;
252 auto pos = SunAtTime(jde).equatorial_coordinate();
253 num_t dec_rad = pos.declination_rad();
255 num_t num = cos(
radians(zenith)) - (sin(dec_rad) * sin(lat_rad));
256 num_t denom = cos(dec_rad) * cos(lat_rad);
257 num_t cos_h = num / denom;
258 if (cos_h > 1 || cos_h < -1)
271 num_t eot = SunAtTime(jd).equation_of_time() * 240;
272 time_t new_timestamp = (time_t) (date.
timestamp + added_d * 86400 - eot);
278 SunAtLocation sun{location_};
280 if (!
m.dt.is_valid())
288 return sun.true_coordinate(
m);
291 SunAtLocation sun{location_};
296 today.
hour = today.minute = today.second = 0;
297 today.recalc_timestamp_utc();
299 auto it = sun.event(rising, today, zenith);
300 if (it.has_value() && it->timestamp < date.
timestamp) {
303 time_t new_timestamp = today.timestamp + 24 * 60 * 60;
305 it = sun.event(rising, today, zenith);
internal::HorizontalCoordinate calc_coords_()
static ESPTime from_epoch_utc(time_t epoch)
Convert an UTC epoch timestamp to a UTC time ESPTime instance.
num_t right_ascension_rad() const
num_t declination_rad() const
num_t delta_t(ESPTime moment)
A more user-friendly version of struct tm from time.h.
optional< ESPTime > calc_event_(bool rising, double zenith)
num_t wmod(num_t x, num_t y)
optional< ESPTime > sunrise(double elevation)
optional< ESPTime > sunset(double elevation)
static ESPTime from_epoch_local(time_t epoch)
Convert an UTC epoch timestamp to a local time ESPTime instance.
num_t julian_day(ESPTime moment)
uint8_t second
seconds after the minute [0-60]
time_t timestamp
unix epoch time (seconds since UTC Midnight January 1, 1970)
uint8_t minute
minutes after the hour [0-59]
num_t elevation_rad() const
bool is_valid() const
Check if this ESPTime is valid (all fields in range and year is greater than 2018) ...
Implementation of SPI Controller mode.
num_t longitude_rad() const
uint8_t month
month; january=1 [1-12]
num_t arcdeg(num_t deg, num_t minutes, num_t seconds)
num_t azimuth_rad() const
uint8_t hour
hours since midnight [0-23]
uint8_t day_of_month
day of the month [1-31]
num_t latitude_rad() const