mirror of
				https://github.com/esphome/esphome.git
				synced 2025-10-30 22:53:59 +00:00 
			
		
		
		
	Rewrite sun component calculations (#1661)
This commit is contained in:
		| @@ -1,3 +1,5 @@ | ||||
| import re | ||||
|  | ||||
| import esphome.codegen as cg | ||||
| import esphome.config_validation as cv | ||||
| from esphome import automation | ||||
| @@ -22,7 +24,7 @@ CONF_ON_SUNSET = "on_sunset" | ||||
|  | ||||
| # Default sun elevation is a bit below horizon because sunset | ||||
| # means time when the entire sun disk is below the horizon | ||||
| DEFAULT_ELEVATION = -0.883 | ||||
| DEFAULT_ELEVATION = -0.83333 | ||||
|  | ||||
| ELEVATION_MAP = { | ||||
|     "sunrise": 0.0, | ||||
| @@ -45,12 +47,54 @@ def elevation(value): | ||||
|     return cv.float_range(min=-180, max=180)(value) | ||||
|  | ||||
|  | ||||
| # Parses sexagesimal values like 22°57′7″S | ||||
| LAT_LON_REGEX = re.compile( | ||||
|     r"([+\-])?\s*" | ||||
|     r"(?:([0-9]+)\s*°)?\s*" | ||||
|     r"(?:([0-9]+)\s*[′\'])?\s*" | ||||
|     r'(?:([0-9]+)\s*[″"])?\s*' | ||||
|     r"([NESW])?" | ||||
| ) | ||||
|  | ||||
|  | ||||
| def parse_latlon(value): | ||||
|     if isinstance(value, str) and value.endswith("°"): | ||||
|         # strip trailing degree character | ||||
|         value = value[:-1] | ||||
|     try: | ||||
|         return cv.float_(value) | ||||
|     except cv.Invalid: | ||||
|         pass | ||||
|  | ||||
|     value = cv.string_strict(value) | ||||
|     m = LAT_LON_REGEX.match(value) | ||||
|  | ||||
|     if m is None: | ||||
|         raise cv.Invalid("Invalid format for latitude/longitude") | ||||
|     sign = m.group(1) | ||||
|     deg = m.group(2) | ||||
|     minute = m.group(3) | ||||
|     second = m.group(4) | ||||
|     d = m.group(5) | ||||
|  | ||||
|     val = float(deg or 0) + float(minute or 0) / 60 + float(second or 0) / 3600 | ||||
|     if sign == "-": | ||||
|         val *= -1 | ||||
|     if d and d in "SW": | ||||
|         val *= -1 | ||||
|     return val | ||||
|  | ||||
|  | ||||
| CONFIG_SCHEMA = cv.Schema( | ||||
|     { | ||||
|         cv.GenerateID(): cv.declare_id(Sun), | ||||
|         cv.GenerateID(CONF_TIME_ID): cv.use_id(time.RealTimeClock), | ||||
|         cv.Required(CONF_LATITUDE): cv.float_range(min=-90, max=90), | ||||
|         cv.Required(CONF_LONGITUDE): cv.float_range(min=-180, max=180), | ||||
|         cv.Required(CONF_LATITUDE): cv.All( | ||||
|             parse_latlon, cv.float_range(min=-90, max=90) | ||||
|         ), | ||||
|         cv.Required(CONF_LONGITUDE): cv.All( | ||||
|             parse_latlon, cv.float_range(min=-180, max=180) | ||||
|         ), | ||||
|         cv.Optional(CONF_ON_SUNRISE): automation.validate_automation( | ||||
|             { | ||||
|                 cv.GenerateID(CONF_TRIGGER_ID): cv.declare_id(SunTrigger), | ||||
|   | ||||
| @@ -1,176 +1,319 @@ | ||||
| #include "sun.h" | ||||
| #include "esphome/core/log.h" | ||||
|  | ||||
| /* | ||||
| The formulas/algorithms in this module are based on the book | ||||
| "Astronomical algorithms" by Jean Meeus (2nd edition) | ||||
|  | ||||
| The target accuracy of this implementation is ~1min for sunrise/sunset calculations, | ||||
| and 6 arcminutes for elevation/azimuth. As such, some of the advanced correction factors | ||||
| like exact nutation are not included. But in some testing the accuracy appears to be within range | ||||
| for random spots around the globe. | ||||
| */ | ||||
|  | ||||
| namespace esphome { | ||||
| namespace sun { | ||||
|  | ||||
| using namespace esphome::sun::internal; | ||||
|  | ||||
| static const char *TAG = "sun"; | ||||
|  | ||||
| #undef PI | ||||
| #undef degrees | ||||
| #undef radians | ||||
| #undef sq | ||||
|  | ||||
| /* Usually, ESPHome uses single-precision floating point values | ||||
|  * because those tend to be accurate enough and are more efficient. | ||||
|  * | ||||
|  * However, some of the data in this class has to be quite accurate, so double is | ||||
|  * used everywhere. | ||||
|  */ | ||||
| static const double PI = 3.141592653589793; | ||||
| static const double TAU = 6.283185307179586; | ||||
| static const double TO_RADIANS = PI / 180.0; | ||||
| static const double TO_DEGREES = 180.0 / PI; | ||||
| static const double EARTH_TILT = 23.44 * TO_RADIANS; | ||||
| static const num_t PI = 3.141592653589793; | ||||
| inline num_t degrees(num_t rad) { return rad * 180 / PI; } | ||||
| inline num_t radians(num_t deg) { return deg * PI / 180; } | ||||
| inline num_t arcdeg(num_t deg, num_t minutes, num_t seconds) { return deg + minutes / 60 + seconds / 3600; } | ||||
| inline num_t sq(num_t x) { return x * x; } | ||||
| inline num_t cb(num_t x) { return x * x * x; } | ||||
|  | ||||
| optional<time::ESPTime> Sun::sunrise(double elevation) { | ||||
|   auto time = this->time_->now(); | ||||
|   if (!time.is_valid()) | ||||
|     return {}; | ||||
|   double sun_time = this->sun_time_for_elevation_(time.day_of_year, elevation, true); | ||||
|   if (isnan(sun_time)) | ||||
|     return {}; | ||||
|   uint32_t epoch = this->calc_epoch_(time, sun_time); | ||||
|   return time::ESPTime::from_epoch_local(epoch); | ||||
| } | ||||
| optional<time::ESPTime> Sun::sunset(double elevation) { | ||||
|   auto time = this->time_->now(); | ||||
|   if (!time.is_valid()) | ||||
|     return {}; | ||||
|   double sun_time = this->sun_time_for_elevation_(time.day_of_year, elevation, false); | ||||
|   if (isnan(sun_time)) | ||||
|     return {}; | ||||
|   uint32_t epoch = this->calc_epoch_(time, sun_time); | ||||
|   return time::ESPTime::from_epoch_local(epoch); | ||||
| } | ||||
| double Sun::elevation() { | ||||
|   auto time = this->current_sun_time_(); | ||||
|   if (isnan(time)) | ||||
|     return NAN; | ||||
|   return this->elevation_(time); | ||||
| } | ||||
| double Sun::azimuth() { | ||||
|   auto time = this->current_sun_time_(); | ||||
|   if (isnan(time)) | ||||
|     return NAN; | ||||
|   return this->azimuth_(time); | ||||
| } | ||||
| // like clamp, but with doubles | ||||
| double clampd(double val, double min, double max) { | ||||
|   if (val < min) | ||||
|     return min; | ||||
|   if (val > max) | ||||
|     return max; | ||||
|   return val; | ||||
| } | ||||
| double Sun::sun_declination_(double sun_time) { | ||||
|   double n = sun_time - 1.0; | ||||
|   // maximum declination | ||||
|   const double tot = -sin(EARTH_TILT); | ||||
| num_t GeoLocation::latitude_rad() const { return radians(latitude); } | ||||
| num_t GeoLocation::longitude_rad() const { return radians(longitude); } | ||||
| num_t EquatorialCoordinate::right_ascension_rad() const { return radians(right_ascension); } | ||||
| num_t EquatorialCoordinate::declination_rad() const { return radians(declination); } | ||||
| num_t HorizontalCoordinate::elevation_rad() const { return radians(elevation); } | ||||
| num_t HorizontalCoordinate::azimuth_rad() const { return radians(azimuth); } | ||||
|  | ||||
|   // eccentricity of the earth's orbit (ellipse) | ||||
|   double eccentricity = 0.0167; | ||||
|  | ||||
|   // days since perihelion (January 3rd) | ||||
|   double days_since_perihelion = n - 2; | ||||
|   // days since december solstice (december 22) | ||||
|   double days_since_december_solstice = n + 10; | ||||
|   const double c = TAU / 365.24; | ||||
|   double v = cos(c * days_since_december_solstice + 2 * eccentricity * sin(c * days_since_perihelion)); | ||||
|   // Make sure value is in range (double error may lead to results slightly larger than 1) | ||||
|   double x = clampd(tot * v, -1.0, 1.0); | ||||
|   return asin(x); | ||||
| num_t julian_day(time::ESPTime moment) { | ||||
|   // p. 59 | ||||
|   // UT -> JD, TT -> JDE | ||||
|   int y = moment.year; | ||||
|   int m = moment.month; | ||||
|   num_t d = moment.day_of_month; | ||||
|   d += moment.hour / 24.0; | ||||
|   d += moment.minute / (24.0 * 60); | ||||
|   d += moment.second / (24.0 * 60 * 60); | ||||
|   if (m <= 2) { | ||||
|     y -= 1; | ||||
|     m += 12; | ||||
|   } | ||||
|   int a = y / 100; | ||||
|   int b = 2 - a + a / 4; | ||||
|   return ((int) (365.25 * (y + 4716))) + ((int) (30.6001 * (m + 1))) + d + b - 1524.5; | ||||
| } | ||||
| double Sun::elevation_ratio_(double sun_time) { | ||||
|   double decl = this->sun_declination_(sun_time); | ||||
|   double hangle = this->hour_angle_(sun_time); | ||||
|   double a = sin(this->latitude_rad_()) * sin(decl); | ||||
|   double b = cos(this->latitude_rad_()) * cos(decl) * cos(hangle); | ||||
|   double val = clampd(a + b, -1.0, 1.0); | ||||
|   return val; | ||||
| num_t delta_t(time::ESPTime moment) { | ||||
|   // approximation for 2005-2050 from NASA (https://eclipse.gsfc.nasa.gov/SEhelp/deltatpoly2004.html) | ||||
|   int t = moment.year - 2000; | ||||
|   return 62.92 + t * (0.32217 + t * 0.005589); | ||||
| } | ||||
| double Sun::latitude_rad_() { return this->latitude_ * TO_RADIANS; } | ||||
| double Sun::hour_angle_(double sun_time) { | ||||
|   double time_of_day = fmod(sun_time, 1.0) * 24.0; | ||||
|   return -PI * (time_of_day - 12) / 12; | ||||
| // Perform a fractional module operation where the result will always be positive (wrapping around) | ||||
| num_t wmod(num_t x, num_t y) { | ||||
|   num_t res = fmod(x, y); | ||||
|   if (res < 0) | ||||
|     res += y; | ||||
|   return res; | ||||
| } | ||||
| double Sun::elevation_(double sun_time) { return this->elevation_rad_(sun_time) * TO_DEGREES; } | ||||
| double Sun::elevation_rad_(double sun_time) { return asin(this->elevation_ratio_(sun_time)); } | ||||
| double Sun::zenith_rad_(double sun_time) { return acos(this->elevation_ratio_(sun_time)); } | ||||
| double Sun::azimuth_rad_(double sun_time) { | ||||
|   double hangle = -this->hour_angle_(sun_time); | ||||
|   double decl = this->sun_declination_(sun_time); | ||||
|   double zen = this->zenith_rad_(sun_time); | ||||
|   double nom = cos(zen) * sin(this->latitude_rad_()) - sin(decl); | ||||
|   double denom = sin(zen) * cos(this->latitude_rad_()); | ||||
|   double v = clampd(nom / denom, -1.0, 1.0); | ||||
|   double az = PI - acos(v); | ||||
|   if (hangle > 0) | ||||
|     az = -az; | ||||
|   if (az < 0) | ||||
|     az += TAU; | ||||
|   return az; | ||||
|  | ||||
| num_t internal::Moment::jd() const { return julian_day(dt); } | ||||
|  | ||||
| num_t internal::Moment::jde() const { | ||||
|   // dt is in UT1, but JDE is based on TT | ||||
|   // so add deltaT factor | ||||
|   return jd() + delta_t(dt) / (60 * 60 * 24); | ||||
| } | ||||
| double Sun::azimuth_(double sun_time) { return this->azimuth_rad_(sun_time) * TO_DEGREES; } | ||||
| double Sun::calc_sun_time_(const time::ESPTime &time) { | ||||
|   // Time as seen at 0° longitude | ||||
|   if (!time.is_valid()) | ||||
|     return NAN; | ||||
|  | ||||
|   double base = (time.day_of_year + time.hour / 24.0 + time.minute / 24.0 / 60.0 + time.second / 24.0 / 60.0 / 60.0); | ||||
|   // Add longitude correction | ||||
|   double add = this->longitude_ / 360.0; | ||||
|   return base + add; | ||||
| } | ||||
| uint32_t Sun::calc_epoch_(time::ESPTime base, double sun_time) { | ||||
|   sun_time -= this->longitude_ / 360.0; | ||||
|   base.day_of_year = uint32_t(floor(sun_time)); | ||||
| struct SunAtTime { | ||||
|   num_t jde; | ||||
|   num_t t; | ||||
|  | ||||
|   sun_time = (sun_time - base.day_of_year) * 24.0; | ||||
|   base.hour = uint32_t(floor(sun_time)); | ||||
|  | ||||
|   sun_time = (sun_time - base.hour) * 60.0; | ||||
|   base.minute = uint32_t(floor(sun_time)); | ||||
|  | ||||
|   sun_time = (sun_time - base.minute) * 60.0; | ||||
|   base.second = uint32_t(floor(sun_time)); | ||||
|  | ||||
|   base.recalc_timestamp_utc(true); | ||||
|   return base.timestamp; | ||||
| } | ||||
| double Sun::sun_time_for_elevation_(int32_t day_of_year, double elevation, bool rising) { | ||||
|   // Use binary search, newton's method would be better but binary search already | ||||
|   // converges quite well (19 cycles) and much simpler. Function is guaranteed to be | ||||
|   // monotonous. | ||||
|   double lo, hi; | ||||
|   if (rising) { | ||||
|     lo = day_of_year + 0.0; | ||||
|     hi = day_of_year + 0.5; | ||||
|   } else { | ||||
|     lo = day_of_year + 1.0; | ||||
|     hi = day_of_year + 0.5; | ||||
|   SunAtTime(num_t jde) : jde(jde) { | ||||
|     // eq 25.1, p. 163; julian centuries from the epoch J2000.0 | ||||
|     t = (jde - 2451545) / 36525.0; | ||||
|   } | ||||
|  | ||||
|   double min_elevation = this->elevation_(lo); | ||||
|   double max_elevation = this->elevation_(hi); | ||||
|   if (elevation < min_elevation || elevation > max_elevation) | ||||
|     return NAN; | ||||
|  | ||||
|   // Accuracy: 0.1s | ||||
|   const double accuracy = 1.0 / (24.0 * 60.0 * 60.0 * 10.0); | ||||
|  | ||||
|   while (fabs(hi - lo) > accuracy) { | ||||
|     double mid = (lo + hi) / 2.0; | ||||
|     double value = this->elevation_(mid) - elevation; | ||||
|     if (value < 0) { | ||||
|       lo = mid; | ||||
|     } else if (value > 0) { | ||||
|       hi = mid; | ||||
|     } else { | ||||
|       lo = hi = mid; | ||||
|       break; | ||||
|     } | ||||
|   num_t mean_obliquity() const { | ||||
|     // eq. 22.2, p. 147; mean obliquity of the ecliptic | ||||
|     num_t epsilon_0 = (+arcdeg(23, 26, 21.448) - arcdeg(0, 0, 46.8150) * t - arcdeg(0, 0, 0.00059) * sq(t) + | ||||
|                        arcdeg(0, 0, 0.001813) * cb(t)); | ||||
|     return epsilon_0; | ||||
|   } | ||||
|  | ||||
|   return (lo + hi) / 2.0; | ||||
|   num_t omega() const { | ||||
|     // eq. 25.8, p. 165; correction factor for obliquity of the ecliptic | ||||
|     // in degrees | ||||
|     num_t omega = 125.05 - 1934.136 * t; | ||||
|     return omega; | ||||
|   } | ||||
|  | ||||
|   num_t true_obliquity() const { | ||||
|     // eq. 25.8, p. 165; correction factor for obliquity of the ecliptic | ||||
|     num_t delta_epsilon = 0.00256 * cos(radians(omega())); | ||||
|     num_t epsilon = mean_obliquity() + delta_epsilon; | ||||
|     return epsilon; | ||||
|   } | ||||
|  | ||||
|   num_t mean_longitude() const { | ||||
|     // eq 25.2, p. 163; geometric mean longitude = mean equinox of the date in degrees | ||||
|     num_t l0 = 280.46646 + 36000.76983 * t + 0.0003032 * sq(t); | ||||
|     return wmod(l0, 360); | ||||
|   } | ||||
|  | ||||
|   num_t eccentricity() const { | ||||
|     // eq 25.4, p. 163; eccentricity of earth's orbit | ||||
|     num_t e = 0.016708634 - 0.000042037 * t - 0.0000001267 * sq(t); | ||||
|     return e; | ||||
|   } | ||||
|  | ||||
|   num_t mean_anomaly() const { | ||||
|     // eq 25.3, p. 163; mean anomaly of the sun in degrees | ||||
|     num_t m = 357.52911 + 35999.05029 * t - 0.0001537 * sq(t); | ||||
|     return wmod(m, 360); | ||||
|   } | ||||
|  | ||||
|   num_t equation_of_center() const { | ||||
|     // p. 164; sun's equation of the center c in degrees | ||||
|     num_t m_rad = radians(mean_anomaly()); | ||||
|     num_t c = ((1.914602 - 0.004817 * t - 0.000014 * sq(t)) * sin(m_rad) + (0.019993 - 0.000101 * t) * sin(2 * m_rad) + | ||||
|                0.000289 * sin(3 * m_rad)); | ||||
|     return wmod(c, 360); | ||||
|   } | ||||
|  | ||||
|   num_t true_longitude() const { | ||||
|     // p. 164; sun's true longitude in degrees | ||||
|     num_t x = mean_longitude() + equation_of_center(); | ||||
|     return wmod(x, 360); | ||||
|   } | ||||
|  | ||||
|   num_t true_anomaly() const { | ||||
|     // p. 164; sun's true anomaly in degrees | ||||
|     num_t x = mean_anomaly() + equation_of_center(); | ||||
|     return wmod(x, 360); | ||||
|   } | ||||
|  | ||||
|   num_t apparent_longitude() const { | ||||
|     // p. 164; sun's apparent longitude = true equinox in degrees | ||||
|     num_t x = true_longitude() - 0.00569 - 0.00478 * sin(radians(omega())); | ||||
|     return wmod(x, 360); | ||||
|   } | ||||
|  | ||||
|   EquatorialCoordinate equatorial_coordinate() const { | ||||
|     num_t epsilon_rad = radians(true_obliquity()); | ||||
|     // eq. 25.6; p. 165; sun's right ascension alpha | ||||
|     num_t app_lon_rad = radians(apparent_longitude()); | ||||
|     num_t right_ascension_rad = atan2(cos(epsilon_rad) * sin(app_lon_rad), cos(app_lon_rad)); | ||||
|     num_t declination_rad = asin(sin(epsilon_rad) * sin(app_lon_rad)); | ||||
|     return EquatorialCoordinate{degrees(right_ascension_rad), degrees(declination_rad)}; | ||||
|   } | ||||
|  | ||||
|   num_t equation_of_time() const { | ||||
|     // chapter 28, p. 185 | ||||
|     num_t epsilon_half = radians(true_obliquity() / 2); | ||||
|     num_t y = sq(tan(epsilon_half)); | ||||
|     num_t l2 = 2 * mean_longitude(); | ||||
|     num_t l2_rad = radians(l2); | ||||
|     num_t e = eccentricity(); | ||||
|     num_t m = mean_anomaly(); | ||||
|     num_t m_rad = radians(m); | ||||
|     num_t sin_m = sin(m_rad); | ||||
|     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) - | ||||
|                  5 / 4.0 * sq(e) * sin(2 * m_rad)); | ||||
|     return degrees(eot); | ||||
|   } | ||||
|  | ||||
|   void debug() const { | ||||
|     // debug output like in example 25.a, p. 165 | ||||
|     ESP_LOGV(TAG, "jde: %f", jde); | ||||
|     ESP_LOGV(TAG, "T: %f", t); | ||||
|     ESP_LOGV(TAG, "L_0: %f", mean_longitude()); | ||||
|     ESP_LOGV(TAG, "M: %f", mean_anomaly()); | ||||
|     ESP_LOGV(TAG, "e: %f", eccentricity()); | ||||
|     ESP_LOGV(TAG, "C: %f", equation_of_center()); | ||||
|     ESP_LOGV(TAG, "Odot: %f", true_longitude()); | ||||
|     ESP_LOGV(TAG, "Omega: %f", omega()); | ||||
|     ESP_LOGV(TAG, "lambda: %f", apparent_longitude()); | ||||
|     ESP_LOGV(TAG, "epsilon_0: %f", mean_obliquity()); | ||||
|     ESP_LOGV(TAG, "epsilon: %f", true_obliquity()); | ||||
|     ESP_LOGV(TAG, "v: %f", true_anomaly()); | ||||
|     auto eq = equatorial_coordinate(); | ||||
|     ESP_LOGV(TAG, "right_ascension: %f", eq.right_ascension); | ||||
|     ESP_LOGV(TAG, "declination: %f", eq.declination); | ||||
|   } | ||||
| }; | ||||
|  | ||||
| struct SunAtLocation { | ||||
|   GeoLocation location; | ||||
|  | ||||
|   num_t greenwich_sidereal_time(Moment moment) const { | ||||
|     // Return the greenwich mean sidereal time for this instant in degrees | ||||
|     // see chapter 12, p. 87 | ||||
|     num_t jd = moment.jd(); | ||||
|     // eq 12.1, p.87; jd for 0h UT of this date | ||||
|     time::ESPTime moment_0h = moment.dt; | ||||
|     moment_0h.hour = moment_0h.minute = moment_0h.second = 0; | ||||
|     num_t jd0 = Moment{moment_0h}.jd(); | ||||
|     num_t t = (jd0 - 2451545) / 36525; | ||||
|     // eq. 12.4, p.88 | ||||
|     num_t gmst = (+280.46061837 + 360.98564736629 * (jd - 2451545) + 0.000387933 * sq(t) - (1 / 38710000.0) * cb(t)); | ||||
|     return wmod(gmst, 360); | ||||
|   } | ||||
|  | ||||
|   HorizontalCoordinate true_coordinate(Moment moment) const { | ||||
|     auto eq = SunAtTime(moment.jde()).equatorial_coordinate(); | ||||
|     num_t gmst = greenwich_sidereal_time(moment); | ||||
|     // do not apply any nutation correction (not important for our target accuracy) | ||||
|     num_t nutation_corr = 0; | ||||
|  | ||||
|     num_t ra = eq.right_ascension; | ||||
|     num_t alpha = gmst + nutation_corr + location.longitude - ra; | ||||
|     alpha = wmod(alpha, 360); | ||||
|     num_t alpha_rad = radians(alpha); | ||||
|  | ||||
|     num_t sin_lat = sin(location.latitude_rad()); | ||||
|     num_t cos_lat = cos(location.latitude_rad()); | ||||
|     num_t sin_elevation = (+sin_lat * sin(eq.declination_rad()) + cos_lat * cos(eq.declination_rad()) * cos(alpha_rad)); | ||||
|     num_t elevation_rad = asin(sin_elevation); | ||||
|     num_t azimuth_rad = atan2(sin(alpha_rad), cos(alpha_rad) * sin_lat - tan(eq.declination_rad()) * cos_lat); | ||||
|     return HorizontalCoordinate{degrees(elevation_rad), degrees(azimuth_rad) + 180}; | ||||
|   } | ||||
|  | ||||
|   optional<time::ESPTime> sunrise(time::ESPTime date, num_t zenith) const { return event(true, date, zenith); } | ||||
|   optional<time::ESPTime> sunset(time::ESPTime date, num_t zenith) const { return event(false, date, zenith); } | ||||
|   optional<time::ESPTime> event(bool rise, time::ESPTime date, num_t zenith) const { | ||||
|     // couldn't get the method described in chapter 15 to work, | ||||
|     // so instead this is based on the algorithm in time4j | ||||
|     // https://github.com/MenoData/Time4J/blob/master/base/src/main/java/net/time4j/calendar/astro/StdSolarCalculator.java | ||||
|     auto m = local_event_(date, 12);  // noon | ||||
|     num_t jde = julian_day(m); | ||||
|     num_t new_h = 0, old_h; | ||||
|     do { | ||||
|       old_h = new_h; | ||||
|       auto x = local_hour_angle_(jde + old_h / 86400, rise, zenith); | ||||
|       if (!x.has_value()) | ||||
|         return {}; | ||||
|       new_h = *x; | ||||
|     } while (std::abs(new_h - old_h) >= 15); | ||||
|     time_t new_timestamp = m.timestamp + (time_t) new_h; | ||||
|     return time::ESPTime::from_epoch_local(new_timestamp); | ||||
|   } | ||||
|  | ||||
|  protected: | ||||
|   optional<num_t> local_hour_angle_(num_t jde, bool rise, num_t zenith) const { | ||||
|     auto pos = SunAtTime(jde).equatorial_coordinate(); | ||||
|     num_t dec_rad = pos.declination_rad(); | ||||
|     num_t lat_rad = location.latitude_rad(); | ||||
|     num_t num = cos(radians(zenith)) - (sin(dec_rad) * sin(lat_rad)); | ||||
|     num_t denom = cos(dec_rad) * cos(lat_rad); | ||||
|     num_t cos_h = num / denom; | ||||
|     if (cos_h > 1 || cos_h < -1) | ||||
|       return {}; | ||||
|     num_t hour_angle = degrees(acos(cos_h)) * 240; | ||||
|     if (rise) | ||||
|       hour_angle *= -1; | ||||
|     return hour_angle; | ||||
|   } | ||||
|  | ||||
|   time::ESPTime local_event_(time::ESPTime date, int hour) const { | ||||
|     // input date should be in UTC, and hour/minute/second fields 0 | ||||
|     num_t added_d = hour / 24.0 - location.longitude / 360; | ||||
|     num_t jd = julian_day(date) + added_d; | ||||
|  | ||||
|     num_t eot = SunAtTime(jd).equation_of_time() * 240; | ||||
|     time_t new_timestamp = (time_t)(date.timestamp + added_d * 86400 - eot); | ||||
|     return time::ESPTime::from_epoch_utc(new_timestamp); | ||||
|   } | ||||
| }; | ||||
|  | ||||
| HorizontalCoordinate Sun::calc_coords_() { | ||||
|   SunAtLocation sun{location_}; | ||||
|   Moment m{time_->utcnow()}; | ||||
|   if (!m.dt.is_valid()) | ||||
|     return HorizontalCoordinate{NAN, NAN}; | ||||
|  | ||||
|   // uncomment to print some debug output | ||||
|   /* | ||||
|   SunAtTime st{m.jde()}; | ||||
|   st.debug(); | ||||
|   */ | ||||
|   return sun.true_coordinate(m); | ||||
| } | ||||
| optional<time::ESPTime> Sun::calc_event_(bool rising, double zenith) { | ||||
|   SunAtLocation sun{location_}; | ||||
|   auto now = this->time_->utcnow(); | ||||
|   if (!now.is_valid()) | ||||
|     return {}; | ||||
|   // Calculate UT1 timestamp at 0h | ||||
|   auto today = now; | ||||
|   today.hour = today.minute = today.second = 0; | ||||
|   today.recalc_timestamp_utc(); | ||||
|  | ||||
|   auto it = sun.event(rising, today, zenith); | ||||
|   if (it.has_value() && it->timestamp < now.timestamp) { | ||||
|     // We're calculating *next* sunrise/sunset, but calculated event | ||||
|     // is today, so try again tomorrow | ||||
|     time_t new_timestamp = today.timestamp + 24 * 60 * 60; | ||||
|     today = time::ESPTime::from_epoch_utc(new_timestamp); | ||||
|     it = sun.event(rising, today, zenith); | ||||
|   } | ||||
|   return it; | ||||
| } | ||||
|  | ||||
| optional<time::ESPTime> Sun::sunrise(double elevation) { return this->calc_event_(true, 90 - elevation); } | ||||
| optional<time::ESPTime> Sun::sunset(double elevation) { return this->calc_event_(false, 90 - elevation); } | ||||
| double Sun::elevation() { return this->calc_coords_().elevation; } | ||||
| double Sun::azimuth() { return this->calc_coords_().azimuth; } | ||||
|  | ||||
| }  // namespace sun | ||||
| }  // namespace esphome | ||||
|   | ||||
| @@ -8,85 +8,72 @@ | ||||
| namespace esphome { | ||||
| namespace sun { | ||||
|  | ||||
| namespace internal { | ||||
|  | ||||
| /* Usually, ESPHome uses single-precision floating point values | ||||
|  * because those tend to be accurate enough and are more efficient. | ||||
|  * | ||||
|  * However, some of the data in this class has to be quite accurate, so double is | ||||
|  * used everywhere. | ||||
|  */ | ||||
| using num_t = double; | ||||
| struct GeoLocation { | ||||
|   num_t latitude; | ||||
|   num_t longitude; | ||||
|  | ||||
|   num_t latitude_rad() const; | ||||
|   num_t longitude_rad() const; | ||||
| }; | ||||
|  | ||||
| struct Moment { | ||||
|   time::ESPTime dt; | ||||
|  | ||||
|   num_t jd() const; | ||||
|   num_t jde() const; | ||||
| }; | ||||
|  | ||||
| struct EquatorialCoordinate { | ||||
|   num_t right_ascension; | ||||
|   num_t declination; | ||||
|  | ||||
|   num_t right_ascension_rad() const; | ||||
|   num_t declination_rad() const; | ||||
| }; | ||||
|  | ||||
| struct HorizontalCoordinate { | ||||
|   num_t elevation; | ||||
|   num_t azimuth; | ||||
|  | ||||
|   num_t elevation_rad() const; | ||||
|   num_t azimuth_rad() const; | ||||
| }; | ||||
|  | ||||
| }  // namespace internal | ||||
|  | ||||
| class Sun { | ||||
|  public: | ||||
|   void set_time(time::RealTimeClock *time) { time_ = time; } | ||||
|   time::RealTimeClock *get_time() const { return time_; } | ||||
|   void set_latitude(double latitude) { latitude_ = latitude; } | ||||
|   void set_longitude(double longitude) { longitude_ = longitude; } | ||||
|   void set_latitude(double latitude) { location_.latitude = latitude; } | ||||
|   void set_longitude(double longitude) { location_.longitude = longitude; } | ||||
|  | ||||
|   optional<time::ESPTime> sunrise(double elevation = 0.0); | ||||
|   optional<time::ESPTime> sunset(double elevation = 0.0); | ||||
|   optional<time::ESPTime> sunrise(double elevation); | ||||
|   optional<time::ESPTime> sunset(double elevation); | ||||
|  | ||||
|   double elevation(); | ||||
|   double azimuth(); | ||||
|  | ||||
|  protected: | ||||
|   double current_sun_time_() { return this->calc_sun_time_(this->time_->utcnow()); } | ||||
|  | ||||
|   /** Calculate the declination of the sun in rad. | ||||
|    * | ||||
|    * See https://en.wikipedia.org/wiki/Position_of_the_Sun#Declination_of_the_Sun_as_seen_from_Earth | ||||
|    * | ||||
|    * Accuracy: ±0.2° | ||||
|    * | ||||
|    * @param sun_time The day of the year, 1 means January 1st. See calc_sun_time_. | ||||
|    * @return Sun declination in degrees | ||||
|    */ | ||||
|   double sun_declination_(double sun_time); | ||||
|  | ||||
|   double elevation_ratio_(double sun_time); | ||||
|  | ||||
|   /** Calculate the hour angle based on the sun time of day in hours. | ||||
|    * | ||||
|    * Positive in morning, 0 at noon, negative in afternoon. | ||||
|    * | ||||
|    * @param sun_time Sun time, see calc_sun_time_. | ||||
|    * @return Hour angle in rad. | ||||
|    */ | ||||
|   double hour_angle_(double sun_time); | ||||
|  | ||||
|   double elevation_(double sun_time); | ||||
|  | ||||
|   double elevation_rad_(double sun_time); | ||||
|  | ||||
|   double zenith_rad_(double sun_time); | ||||
|  | ||||
|   double azimuth_rad_(double sun_time); | ||||
|  | ||||
|   double azimuth_(double sun_time); | ||||
|  | ||||
|   /** Return the sun time given by the time_ object. | ||||
|    * | ||||
|    * Sun time is defined as doubleing point day of year. | ||||
|    * Integer part encodes the day of the year (1=January 1st) | ||||
|    * Decimal part encodes time of day (1/24 = 1 hour) | ||||
|    */ | ||||
|   double calc_sun_time_(const time::ESPTime &time); | ||||
|  | ||||
|   uint32_t calc_epoch_(time::ESPTime base, double sun_time); | ||||
|  | ||||
|   /** Calculate the sun time of day | ||||
|    * | ||||
|    * @param day_of_year | ||||
|    * @param elevation | ||||
|    * @param rising | ||||
|    * @return | ||||
|    */ | ||||
|   double sun_time_for_elevation_(int32_t day_of_year, double elevation, bool rising); | ||||
|  | ||||
|   double latitude_rad_(); | ||||
|   internal::HorizontalCoordinate calc_coords_(); | ||||
|   optional<time::ESPTime> calc_event_(bool rising, double zenith); | ||||
|  | ||||
|   time::RealTimeClock *time_; | ||||
|   /// Latitude in degrees, range: -90 to 90. | ||||
|   double latitude_; | ||||
|   /// Longitude in degrees, range: -180 to 180. | ||||
|   double longitude_; | ||||
|   internal::GeoLocation location_; | ||||
| }; | ||||
|  | ||||
| class SunTrigger : public Trigger<>, public PollingComponent, public Parented<Sun> { | ||||
|  public: | ||||
|   SunTrigger() : PollingComponent(1000) {} | ||||
|   SunTrigger() : PollingComponent(60000) {} | ||||
|  | ||||
|   void set_sunrise(bool sunrise) { sunrise_ = sunrise; } | ||||
|   void set_elevation(double elevation) { elevation_ = elevation; } | ||||
|   | ||||
		Reference in New Issue
	
	Block a user