diff --git a/esphome/components/sun/__init__.py b/esphome/components/sun/__init__.py
index 5241f1bb55..b24d62f32f 100644
--- a/esphome/components/sun/__init__.py
+++ b/esphome/components/sun/__init__.py
@@ -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),
diff --git a/esphome/components/sun/sun.cpp b/esphome/components/sun/sun.cpp
index 6744a418c5..0722e66e67 100644
--- a/esphome/components/sun/sun.cpp
+++ b/esphome/components/sun/sun.cpp
@@ -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
diff --git a/esphome/components/sun/sun.h b/esphome/components/sun/sun.h
index 501d122da0..6a8364a5f0 100644
--- a/esphome/components/sun/sun.h
+++ b/esphome/components/sun/sun.h
@@ -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; }
diff --git a/script/ci-custom.py b/script/ci-custom.py
index 3106cabeff..4ec7c664a4 100755
--- a/script/ci-custom.py
+++ b/script/ci-custom.py
@@ -405,6 +405,7 @@ ARDUINO_FORBIDDEN_RE = r"[^\w\d](" + r"|".join(ARDUINO_FORBIDDEN) + r")\(.*"
     include=cpp_include,
     exclude=[
         "esphome/components/mqtt/custom_mqtt_device.h",
+        "esphome/components/sun/sun.cpp",
         "esphome/core/esphal.*",
     ],
 )