ESPHome  2024.5.1
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
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
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
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 {
152  // eq. 25.6; p. 165; sun's right ascension alpha
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();
165  num_t e = eccentricity();
166  num_t m = mean_anomaly();
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);
221
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();
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
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
uint8_t m
Definition: bl0939.h:20
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
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
This is a workaround until we can figure out a way to get the tflite-micro idf component code availab...
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 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