ESPHome  2024.10.2
sun.cpp
Go to the documentation of this file.
1 #include "sun.h"
2 #include "esphome/core/log.h"
3 
4 /*
5 The formulas/algorithms in this module are based on the book
6 "Astronomical algorithms" by Jean Meeus (2nd edition)
7 
8 The target accuracy of this implementation is ~1min for sunrise/sunset calculations,
9 and 6 arcminutes for elevation/azimuth. As such, some of the advanced correction factors
10 like exact nutation are not included. But in some testing the accuracy appears to be within range
11 for random spots around the globe.
12 */
13 
14 namespace esphome {
15 namespace sun {
16 
17 using namespace esphome::sun::internal;
18 
19 static const char *const TAG = "sun";
20 
21 #undef PI
22 #undef degrees
23 #undef radians
24 #undef sq
25 
26 static const num_t PI = 3.141592653589793;
27 inline num_t degrees(num_t rad) { return rad * 180 / PI; }
28 inline num_t radians(num_t deg) { return deg * PI / 180; }
29 inline num_t arcdeg(num_t deg, num_t minutes, num_t seconds) { return deg + minutes / 60 + seconds / 3600; }
30 inline num_t sq(num_t x) { return x * x; }
31 inline num_t cb(num_t x) { return x * x * x; }
32 
35 num_t EquatorialCoordinate::right_ascension_rad() const { return radians(right_ascension); }
36 num_t EquatorialCoordinate::declination_rad() const { return radians(declination); }
37 num_t HorizontalCoordinate::elevation_rad() const { return radians(elevation); }
38 num_t HorizontalCoordinate::azimuth_rad() const { return radians(azimuth); }
39 
41  // p. 59
42  // UT -> JD, TT -> JDE
43  int y = moment.year;
44  int m = moment.month;
45  num_t d = moment.day_of_month;
46  d += moment.hour / 24.0;
47  d += moment.minute / (24.0 * 60);
48  d += moment.second / (24.0 * 60 * 60);
49  if (m <= 2) {
50  y -= 1;
51  m += 12;
52  }
53  int a = y / 100;
54  int b = 2 - a + a / 4;
55  return ((int) (365.25 * (y + 4716))) + ((int) (30.6001 * (m + 1))) + d + b - 1524.5;
56 }
58  // approximation for 2005-2050 from NASA (https://eclipse.gsfc.nasa.gov/SEhelp/deltatpoly2004.html)
59  int t = moment.year - 2000;
60  return 62.92 + t * (0.32217 + t * 0.005589);
61 }
62 // Perform a fractional module operation where the result will always be positive (wrapping around)
64  num_t res = fmod(x, y);
65  if (res < 0)
66  res += y;
67  return res;
68 }
69 
70 num_t internal::Moment::jd() const { return julian_day(dt); }
71 
73  // dt is in UT1, but JDE is based on TT
74  // so add deltaT factor
75  return jd() + delta_t(dt) / (60 * 60 * 24);
76 }
77 
78 struct SunAtTime {
79  num_t jde;
80  num_t t;
81 
82  // eq 25.1, p. 163; julian centuries from the epoch J2000.0
83  SunAtTime(num_t jde) : jde(jde), t((jde - 2451545) / 36525.0) {}
84 
85  num_t mean_obliquity() const {
86  // eq. 22.2, p. 147; mean obliquity of the ecliptic
87  num_t epsilon_0 = (+arcdeg(23, 26, 21.448) - arcdeg(0, 0, 46.8150) * t - arcdeg(0, 0, 0.00059) * sq(t) +
88  arcdeg(0, 0, 0.001813) * cb(t));
89  return epsilon_0;
90  }
91 
92  num_t omega() const {
93  // eq. 25.8, p. 165; correction factor for obliquity of the ecliptic
94  // in degrees
95  num_t omega = 125.05 - 1934.136 * t;
96  return omega;
97  }
98 
99  num_t true_obliquity() const {
100  // eq. 25.8, p. 165; correction factor for obliquity of the ecliptic
101  num_t delta_epsilon = 0.00256 * cos(radians(omega()));
102  num_t epsilon = mean_obliquity() + delta_epsilon;
103  return epsilon;
104  }
105 
106  num_t mean_longitude() const {
107  // eq 25.2, p. 163; geometric mean longitude = mean equinox of the date in degrees
108  num_t l0 = 280.46646 + 36000.76983 * t + 0.0003032 * sq(t);
109  return wmod(l0, 360);
110  }
111 
112  num_t eccentricity() const {
113  // eq 25.4, p. 163; eccentricity of earth's orbit
114  num_t e = 0.016708634 - 0.000042037 * t - 0.0000001267 * sq(t);
115  return e;
116  }
117 
118  num_t mean_anomaly() const {
119  // eq 25.3, p. 163; mean anomaly of the sun in degrees
120  num_t m = 357.52911 + 35999.05029 * t - 0.0001537 * sq(t);
121  return wmod(m, 360);
122  }
123 
124  num_t equation_of_center() const {
125  // p. 164; sun's equation of the center c in degrees
126  num_t m_rad = radians(mean_anomaly());
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));
129  return wmod(c, 360);
130  }
131 
132  num_t true_longitude() const {
133  // p. 164; sun's true longitude in degrees
134  num_t x = mean_longitude() + equation_of_center();
135  return wmod(x, 360);
136  }
137 
138  num_t true_anomaly() const {
139  // p. 164; sun's true anomaly in degrees
140  num_t x = mean_anomaly() + equation_of_center();
141  return wmod(x, 360);
142  }
143 
144  num_t apparent_longitude() const {
145  // p. 164; sun's apparent longitude = true equinox in degrees
146  num_t x = true_longitude() - 0.00569 - 0.00478 * sin(radians(omega()));
147  return wmod(x, 360);
148  }
149 
150  EquatorialCoordinate equatorial_coordinate() const {
151  num_t epsilon_rad = radians(true_obliquity());
152  // eq. 25.6; p. 165; sun's right ascension alpha
153  num_t app_lon_rad = radians(apparent_longitude());
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));
156  return EquatorialCoordinate{degrees(right_ascension_rad), degrees(declination_rad)};
157  }
158 
159  num_t equation_of_time() const {
160  // chapter 28, p. 185
161  num_t epsilon_half = radians(true_obliquity() / 2);
162  num_t y = sq(tan(epsilon_half));
163  num_t l2 = 2 * mean_longitude();
164  num_t l2_rad = radians(l2);
165  num_t e = eccentricity();
166  num_t m = mean_anomaly();
167  num_t m_rad = radians(m);
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));
171  return degrees(eot);
172  }
173 
174  void debug() const {
175  // debug output like in example 25.a, p. 165
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);
191  }
192 };
193 
194 struct SunAtLocation {
195  GeoLocation location;
196 
197  num_t greenwich_sidereal_time(Moment moment) const {
198  // Return the greenwich mean sidereal time for this instant in degrees
199  // see chapter 12, p. 87
200  num_t jd = moment.jd();
201  // eq 12.1, p.87; jd for 0h UT of this date
202  ESPTime moment_0h = moment.dt;
203  moment_0h.hour = moment_0h.minute = moment_0h.second = 0;
204  num_t jd0 = Moment{moment_0h}.jd();
205  num_t t = (jd0 - 2451545) / 36525;
206  // eq. 12.4, p.88
207  num_t gmst = (+280.46061837 + 360.98564736629 * (jd - 2451545) + 0.000387933 * sq(t) - (1 / 38710000.0) * cb(t));
208  return wmod(gmst, 360);
209  }
210 
211  HorizontalCoordinate true_coordinate(Moment moment) const {
212  auto eq = SunAtTime(moment.jde()).equatorial_coordinate();
213  num_t gmst = greenwich_sidereal_time(moment);
214  // do not apply any nutation correction (not important for our target accuracy)
215  num_t nutation_corr = 0;
216 
217  num_t ra = eq.right_ascension;
218  num_t alpha = gmst + nutation_corr + location.longitude - ra;
219  alpha = wmod(alpha, 360);
220  num_t alpha_rad = radians(alpha);
221 
222  num_t sin_lat = sin(location.latitude_rad());
223  num_t cos_lat = cos(location.latitude_rad());
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);
227  return HorizontalCoordinate{degrees(elevation_rad), degrees(azimuth_rad) + 180};
228  }
229 
230  optional<ESPTime> sunrise(ESPTime date, num_t zenith) const { return event(true, date, zenith); }
231  optional<ESPTime> sunset(ESPTime date, num_t zenith) const { return event(false, date, zenith); }
232  optional<ESPTime> event(bool rise, ESPTime date, num_t zenith) const {
233  // couldn't get the method described in chapter 15 to work,
234  // so instead this is based on the algorithm in time4j
235  // https://github.com/MenoData/Time4J/blob/master/base/src/main/java/net/time4j/calendar/astro/StdSolarCalculator.java
236  auto m = local_event_(date, 12); // noon
237  num_t jde = julian_day(m);
238  num_t new_h = 0, old_h;
239  do {
240  old_h = new_h;
241  auto x = local_hour_angle_(jde + old_h / 86400, rise, zenith);
242  if (!x.has_value())
243  return {};
244  new_h = *x;
245  } while (std::abs(new_h - old_h) >= 15);
246  time_t new_timestamp = m.timestamp + (time_t) new_h;
247  return ESPTime::from_epoch_local(new_timestamp);
248  }
249 
250  protected:
251  optional<num_t> local_hour_angle_(num_t jde, bool rise, num_t zenith) const {
252  auto pos = SunAtTime(jde).equatorial_coordinate();
253  num_t dec_rad = pos.declination_rad();
254  num_t lat_rad = location.latitude_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)
259  return {};
260  num_t hour_angle = degrees(acos(cos_h)) * 240;
261  if (rise)
262  hour_angle *= -1;
263  return hour_angle;
264  }
265 
266  ESPTime local_event_(ESPTime date, int hour) const {
267  // input date should be in UTC, and hour/minute/second fields 0
268  num_t added_d = hour / 24.0 - location.longitude / 360;
269  num_t jd = julian_day(date) + added_d;
270 
271  num_t eot = SunAtTime(jd).equation_of_time() * 240;
272  time_t new_timestamp = (time_t) (date.timestamp + added_d * 86400 - eot);
273  return ESPTime::from_epoch_utc(new_timestamp);
274  }
275 };
276 
278  SunAtLocation sun{location_};
279  Moment m{time_->utcnow()};
280  if (!m.dt.is_valid())
281  return HorizontalCoordinate{NAN, NAN};
282 
283  // uncomment to print some debug output
284  /*
285  SunAtTime st{m.jde()};
286  st.debug();
287  */
288  return sun.true_coordinate(m);
289 }
290 optional<ESPTime> Sun::calc_event_(ESPTime date, bool rising, double zenith) {
291  SunAtLocation sun{location_};
292  if (!date.is_valid())
293  return {};
294  // Calculate UT1 timestamp at 0h
295  auto today = date;
296  today.hour = today.minute = today.second = 0;
297  today.recalc_timestamp_utc();
298 
299  auto it = sun.event(rising, today, zenith);
300  if (it.has_value() && it->timestamp < date.timestamp) {
301  // We're calculating *next* sunrise/sunset, but calculated event
302  // is today, so try again tomorrow
303  time_t new_timestamp = today.timestamp + 24 * 60 * 60;
304  today = ESPTime::from_epoch_utc(new_timestamp);
305  it = sun.event(rising, today, zenith);
306  }
307  return it;
308 }
309 optional<ESPTime> Sun::calc_event_(bool rising, double zenith) {
310  auto it = Sun::calc_event_(this->time_->utcnow(), rising, zenith);
311  return it;
312 }
313 
314 optional<ESPTime> Sun::sunrise(double elevation) { return this->calc_event_(true, 90 - elevation); }
315 optional<ESPTime> Sun::sunset(double elevation) { return this->calc_event_(false, 90 - elevation); }
316 optional<ESPTime> Sun::sunrise(ESPTime date, double elevation) { return this->calc_event_(date, true, 90 - elevation); }
317 optional<ESPTime> Sun::sunset(ESPTime date, double elevation) { return this->calc_event_(date, false, 90 - elevation); }
318 double Sun::elevation() { return this->calc_coords_().elevation; }
319 double Sun::azimuth() { return this->calc_coords_().azimuth; }
320 
321 } // namespace sun
322 } // namespace esphome
internal::HorizontalCoordinate calc_coords_()
Definition: sun.cpp:277
static ESPTime from_epoch_utc(time_t epoch)
Convert an UTC epoch timestamp to a UTC time ESPTime instance.
Definition: time.h:94
double azimuth()
Definition: sun.cpp:319
num_t delta_t(ESPTime moment)
Definition: sun.cpp:57
uint16_t x
Definition: tt21100.cpp:17
A more user-friendly version of struct tm from time.h.
Definition: time.h:17
num_t degrees(num_t rad)
Definition: sun.cpp:27
optional< ESPTime > calc_event_(bool rising, double zenith)
Definition: sun.cpp:309
num_t wmod(num_t x, num_t y)
Definition: sun.cpp:63
double elevation()
Definition: sun.cpp:318
uint16_t y
Definition: tt21100.cpp:18
optional< ESPTime > sunrise(double elevation)
Definition: sun.cpp:314
const char *const TAG
Definition: spi.cpp:8
optional< ESPTime > sunset(double elevation)
Definition: sun.cpp:315
static ESPTime from_epoch_local(time_t epoch)
Convert an UTC epoch timestamp to a local time ESPTime instance.
Definition: time.h:85
num_t sq(num_t x)
Definition: sun.cpp:30
num_t julian_day(ESPTime moment)
Definition: sun.cpp:40
uint8_t second
seconds after the minute [0-60]
Definition: time.h:21
time_t timestamp
unix epoch time (seconds since UTC Midnight January 1, 1970)
Definition: time.h:39
uint8_t hour
num_t radians(num_t deg)
Definition: sun.cpp:28
uint8_t minute
minutes after the hour [0-59]
Definition: time.h:23
bool is_valid() const
Check if this ESPTime is valid (all fields in range and year is greater than 2018) ...
Definition: time.h:61
uint16_t year
year
Definition: time.h:35
Implementation of SPI Controller mode.
Definition: a01nyub.cpp:7
num_t cb(num_t x)
Definition: sun.cpp:31
uint8_t month
month; january=1 [1-12]
Definition: time.h:33
num_t arcdeg(num_t deg, num_t minutes, num_t seconds)
Definition: sun.cpp:29
uint8_t m
Definition: bl0906.h:208
uint8_t hour
hours since midnight [0-23]
Definition: time.h:25
uint8_t day_of_month
day of the month [1-31]
Definition: time.h:29