ESPHome  2022.11.3
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  time::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<time::ESPTime> sunrise(time::ESPTime date, num_t zenith) const { return event(true, date, zenith); }
231  optional<time::ESPTime> sunset(time::ESPTime date, num_t zenith) const { return event(false, date, zenith); }
232  optional<time::ESPTime> event(bool rise, time::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 time::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  time::ESPTime local_event_(time::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 time::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<time::ESPTime> Sun::calc_event_(bool rising, double zenith) {
291  SunAtLocation sun{location_};
292  auto now = this->time_->utcnow();
293  if (!now.is_valid())
294  return {};
295  // Calculate UT1 timestamp at 0h
296  auto today = now;
297  today.hour = today.minute = today.second = 0;
298  today.recalc_timestamp_utc();
299 
300  auto it = sun.event(rising, today, zenith);
301  if (it.has_value() && it->timestamp < now.timestamp) {
302  // We're calculating *next* sunrise/sunset, but calculated event
303  // is today, so try again tomorrow
304  time_t new_timestamp = today.timestamp + 24 * 60 * 60;
305  today = time::ESPTime::from_epoch_utc(new_timestamp);
306  it = sun.event(rising, today, zenith);
307  }
308  return it;
309 }
310 
311 optional<time::ESPTime> Sun::sunrise(double elevation) { return this->calc_event_(true, 90 - elevation); }
312 optional<time::ESPTime> Sun::sunset(double elevation) { return this->calc_event_(false, 90 - elevation); }
313 double Sun::elevation() { return this->calc_coords_().elevation; }
314 double Sun::azimuth() { return this->calc_coords_().azimuth; }
315 
316 } // namespace sun
317 } // namespace esphome
internal::HorizontalCoordinate calc_coords_()
Definition: sun.cpp:277
uint8_t month
month; january=1 [1-12]
double azimuth()
Definition: sun.cpp:314
time_t timestamp
unix epoch time (seconds since UTC Midnight January 1, 1970)
num_t degrees(num_t rad)
Definition: sun.cpp:27
num_t wmod(num_t x, num_t y)
Definition: sun.cpp:63
double elevation()
Definition: sun.cpp:313
static ESPTime from_epoch_local(time_t epoch)
Convert an UTC epoch timestamp to a local time ESPTime instance.
uint8_t minute
minutes after the hour [0-59]
optional< time::ESPTime > sunset(double elevation)
Definition: sun.cpp:312
num_t julian_day(time::ESPTime moment)
Definition: sun.cpp:40
uint8_t m
Definition: bl0939.h:20
optional< time::ESPTime > calc_event_(bool rising, double zenith)
Definition: sun.cpp:290
A more user-friendly version of struct tm from time.h.
num_t sq(num_t x)
Definition: sun.cpp:30
static ESPTime from_epoch_utc(time_t epoch)
Convert an UTC epoch timestamp to a UTC time ESPTime instance.
num_t radians(num_t deg)
Definition: sun.cpp:28
uint8_t second
seconds after the minute [0-60]
optional< time::ESPTime > sunrise(double elevation)
Definition: sun.cpp:311
Definition: a4988.cpp:4
num_t cb(num_t x)
Definition: sun.cpp:31
uint8_t day_of_month
day of the month [1-31]
num_t arcdeg(num_t deg, num_t minutes, num_t seconds)
Definition: sun.cpp:29
time::ESPTime dt
Definition: sun.h:29
uint8_t hour
hours since midnight [0-23]
num_t delta_t(time::ESPTime moment)
Definition: sun.cpp:57