ENG | Calculating sunrise and sunset

I asked AI how to compute sunrise and sunset and to generate code in Julia. I don’t use this language beyond better calculator, but it’s probably the most readable and actually suitable for this stuff. It generated functional code which had errors in order of minutes. Then I asked it why it’s inaccurate, it generated better code with error of minute or two. It was not that long - like 400 lines with lot of comments, nice output formatting, it included twilight times.

But then what to even use as a reference? timeanddate.com, Garmin smart watch? Why do they have one minute differences?

Known facts

• Earth travels on ellipsis with slight variation in distance from sun and angular speed - roughly 1° per day
• Earth rotates around is axis which has (almost) fixed orientation towards stars
• Earth axis is tilted by 23.44° relative to its orbital plane
• Earth rotation relative to stars (sideric day has 23h56m4.09s) and sun (24h) differs
• Earth rotates around sun on ellipsis going through nearest point every xxx days

Weird facts

• Sideric year: 365.256363 days (full rotation around sun) - but one extra rotation towards stars
• Tropical year: 365.242190 days (from spring equinox to spring equinox)
• Somehow difference is ~20 minutes yearly, 1 day in ~70 years and one year in ~26000 years -> polaris won’t be northern star forever, in 13000 years Orion can be summer constelation and earth will be closest to sun in July.
• Tilt of axis changes a bit

I assumed that sunrise/sunset could not be that hard cause earth moves on ellipsis so we need to solve position relative to sun, then it rotates around axis with period relative stars and we can simply solve position and at given time compute azimuth and elevation of sun at given location. Well, for some reason it’s not that simple.

Position on ellipsis can’t be directly solved so it seems to be approximated by 2nd or 3rd order polynom fit.

So where errors possibly come from?

• Rounding seconds to minutes
• Polynom fits
• Ignoring leap years
• Sun crossing horizon (center/top, athmopheric refraction correction - which, sadly, depend on current weather)

In the end I dropped AI generated Julia code almost entirely and use Excel sheet from NOAA as a starting point. Strength of Julia code was it’s structure.

Time conversion

We have various time references. Curiously, Excel has 1900-01-01 as day 1 for backward compatibility with Lotus 1-2-3 (some may remember Lotus AmiPro or Lotus Organizer from 90s). Except Lotus 1-2-3 had a bug, that considered 1900 as a leap year and Excel keeps this. That’s why day 0 in Excel is 1899-12-30.

Astronomers use an even weirder system: the Julian Day Number (JD), which counts days from January 1, 4713 BCE at noon—a date chosen because it’s when three ancient calendar cycles (28-year solar, 19-year lunar, 15-year Roman tax) all aligned at their starting point. It predates all recorded history, making it perfect for astronomical calculations without worrying about calendar reforms. And noon, cause whole night has the same date.

• 2400000.5 1858-11-17 00:00 Modified Julian day (shorter)
• 2415018.5 1899-12-30 00:00 (1-2-3, Excel epoch)
• 2451545.0 2000-01-01 12:00 (J2000.0 epoch)

Not to mention leap seconds and various different continous atomic times.

Code snippets

def excel_date_to_julian_day(excel_date, timezone_offset):
    # Excel's date system starts at 1900-01-01 (Excel day 1)
    # But Excel incorrectly assumes 1900 is a leap year, so day 60 is 1900-02-29 (which doesn't exist)
    # So we adjust for this bug (comes from Lotus-1-2-3)
    excel_epoch = datetime(1899, 12, 30)
    delta = timedelta(days=excel_date)  # Subtract 2 to correct for Excel's 1900 leap year bug
    date = excel_epoch + delta

    # Calculate Julian Day
    # Excel formula: JD = DATE + 2415018.5 + DAY_FRACTION - TIMEZONE/24
    # Here, DAY_FRACTION is 0.5 for noon, but you can adjust as needed
    jd = excel_date + 2415018.5 - timezone_offset / 24
    return jd

def date_to_excel_serial_date(date):
    # Excel epoch (1900-01-01) with the 1900 leap year bug correction
    excel_epoch = datetime(1899, 12, 30)
    delta = date - excel_epoch
    return delta.days

def date_to_julian_day(date, timezone_offset):
    # https://en.wikipedia.org/wiki/Julian_day#History yes, we have year 6738 in 2025
    excel_date = date_to_excel_serial_date(date)
    jd = excel_date_to_julian_day(excel_date, timezone_offset)    
    return jd

FIXME: date, make it function

# Výpočet Julian Century (G)
julian_century      = (julian_day(datetime(2025, 12, 24), timezone_offset) + time_past_midnight - 2451545)/36525    

Because nothing is nice and constant, there are some weird approximation by polynoms of 2nd or 3rd degree. In the end we need sun_declination (angle from equator) and equation_of_time which is offset of solar noon (sundial time vs mechanical clock).

These depend only on time. Note that this equation gives different result for morning, noon, evening and midnight and it also ignores timezones. So it’s valid for every place on earth at any given moment.

Curiously equation of time can differ by tens seconds per day (winter days are roughly 30s longer than average)

Getting equation of time and declination

def precalc(julian_century):
    # Geometrická střední délka slunce (I)
    geom_mean_long_sun = (280.46646 + julian_century * (36000.76983 + julian_century * 0.0003032)) % 360

    # Geometrická střední anomálie slunce (J)
    geom_mean_anom_sun = 357.52911 + julian_century * (35999.05029 - 0.0001537 * julian_century)

    # Excentricita zemské dráhy (K)
    eccent_earth_orbit = 0.016708634 - julian_century * (0.000042037 + 0.0000001267 * julian_century)

    # Rovnice středu (L)
    sun_eq_of_ctr = (
          math.sin(math.radians(    geom_mean_anom_sun)) * (1.914602 - julian_century * (0.004817 + 0.000014 * julian_century))
        + math.sin(math.radians(2 * geom_mean_anom_sun)) * (0.019993 - 0.000101 * julian_century)
        + math.sin(math.radians(3 * geom_mean_anom_sun)) * 0.000289
    )

    # Pravá délka slunce (M)
    sun_true_long = geom_mean_long_sun + sun_eq_of_ctr

    # Pravá anomálie slunce (N)
    sun_true_anom = geom_mean_anom_sun + sun_eq_of_ctr

    # Vzdálenost slunce od Země (O)
    sun_rad_vector = (1.000001018 * (1 - eccent_earth_orbit**2)) / (1 + eccent_earth_orbit * math.cos(math.radians(sun_true_anom)))

    # Zdánlivá délka slunce (P)
    sun_app_long = sun_true_long - 0.00569 - 0.00478 * math.sin(math.radians(125.04 - 1934.136 * julian_century))

    # Střední sklon ekliptiky (Q)
    mean_obliq_ecliptic = 23 + (26 + ((21.448 - julian_century * (46.815 + julian_century * (0.00059 - julian_century * 0.001813)))) / 60) / 60

    # Korekce sklonu ekliptiky (R)
    obliq_corr = mean_obliq_ecliptic + 0.00256 * math.cos(math.radians(125.04 - 1934.136 * julian_century))

    # Rektascenze slunce (S)
    # https://en.wikipedia.org/wiki/Right_ascension
    sun_rt_ascen = math.degrees(
        math.atan2(
            math.cos(math.radians(sun_app_long)),
            math.cos(math.radians(obliq_corr)) * math.sin(math.radians(sun_app_long)),
        )
    )

    # Deklinace slunce (T)
    sun_declin = math.degrees(math.asin(math.sin(math.radians(obliq_corr)) * math.sin(math.radians(sun_app_long))))

    # Pomocná proměnná (U)
    var_y = math.tan(math.radians(obliq_corr / 2)) ** 2

    # Rovnice času (V)
    eq_of_time = 4 * math.degrees(
        var_y * math.sin(2 * math.radians(geom_mean_long_sun))
        - 2 * eccent_earth_orbit * math.sin(math.radians(geom_mean_anom_sun))
        + 4 * eccent_earth_orbit * var_y * math.sin(math.radians(geom_mean_anom_sun)) * math.cos(2 * math.radians(geom_mean_long_sun))
        - 0.5 * var_y**2 * math.sin(4 * math.radians(geom_mean_long_sun))
        - 1.25 * eccent_earth_orbit**2 * math.sin(2 * math.radians(geom_mean_anom_sun))
    )

    return (sun_declin, eq_of_time)

We can precompute sun declination and equation of time, because this part is the most difficult calculation.

Solar position at time & Time of solar elevation

Then there are relatively simple functions that depend on latitude, longitude and can give answer to how early before noon sun reaches certain elevation over/below horizon (to calculate twilight, sunrise, golden hour) and second function is for solar position.

def hour_angle_for_elevation(latitude, declination, elevation):
    lat_rad  = math.radians(latitude)
    dec_rad  = math.radians(declination)
    elev_rad = math.radians(elevation)
    
    # Spherical astronomy formula
    cos_h = (math.sin(elev_rad) - math.sin(lat_rad) * math.sin(dec_rad)) / (math.cos(lat_rad) * math.cos(dec_rad))
    
    # Check if sun reaches this elevation today
    if cos_h < -1.0:
        return None  # Sun always above this elevation (polar summer)
    elif cos_h > 1.0:
        return None  # Sun never reaches this elevation (polar winter)

    return math.degrees(math.acos(cos_h))


def solar_position(latitude, sun_declin, eq_of_time, local_hour):
    """
    Calculate solar azimuth and elevation at given time

    Parameters:
    - latitude: observer latitude (degrees)
    - doy: day of year
    - solar_hour: time in solar hours (0-24)

    Returns:
    - (azimuth, elevation) in degrees
    - azimuth: 0° = North, 90° = East, 180° = South, 270° = West
    - elevation: angle above horizon (negative = below horizon)
    """

    # Hour angle from solar noon (15° per hour)
    solar_hour = local_hour + longitude/15.0 - timezone_offset
    hour_angle = 15.0 * (solar_hour - 12.0 + eq_of_time/60.0)
    
    lat_rad = math.radians(latitude)
    dec_rad = math.radians(sun_declin)
    h_rad = math.radians(hour_angle)
    
    # Solar elevation (altitude)
    sin_elev = math.sin(lat_rad) * math.sin(dec_rad) + \
               math.cos(lat_rad) * math.cos(dec_rad) * math.cos(h_rad)
    elevation = math.degrees(math.asin(sin_elev))
    
    # Solar azimuth
    cos_azim = (math.sin(dec_rad) - math.sin(lat_rad) * sin_elev) / \
               (math.cos(lat_rad) * math.cos(math.asin(sin_elev)))
    # Clamp to prevent numerical errors
    #cos_azim = math.clamp(cos_azim, -1.0, 1.0)
    if cos_azim > 1:
        cos_azim = 1
    if cos_azim < -1:
        cos_azim = -1
    azimuth = math.degrees(math.acos(cos_azim))
    
    # Adjust azimuth for afternoon (hour angle > 0)
    if hour_angle > 0:
        azimuth = 360.0 - azimuth

    return (azimuth, elevation)

Test code

Test code may look like this:

# Brno, CZ, Náměstí svobody
latitude  = 49.19528 
longitude = 16.60778

timezone_offset = 1   # UTC+1, CET

time_past_midnight = 0.5

# Precompute array of sun_declin, eq_of_time
start_date = datetime(2025, 12, 1)
end_date   = datetime(2025, 12, 31)

current_date = start_date
julian_century = date_to_julian_century(start_date, timezone_offset, time_past_midnight)
array = []
while current_date <= end_date:
    sun_declin, eq_of_time = precalc(julian_century)
    array.append((sun_declin, eq_of_time))
    current_date   += timedelta(days=1)
    julian_century += 1/DAYS_PER_CENTURY


for index, (sun_declin, eq_of_time) in enumerate(array):
    sun_declin, eq_of_time = array[index]
    current_date = start_date + timedelta(days=index)
    ha_sunrise = hour_angle_for_elevation(latitude, sun_declin, -0.8333)

    # Sluneční poledne (X)
    solar_noon = (720 - 4 * longitude - eq_of_time + timezone_offset * 60) / 1440

    # Čas východu slunce (Y)
    sunrise_time = solar_noon - ha_sunrise * 4 / 1440

    # Čas západu slunce (Z)
    sunset_time = solar_noon + ha_sunrise * 4 / 1440
    print(f"{current_date:%Y-%m-%d}\t{fraction_to_time(sunrise_time)}\t{fraction_to_time(solar_noon)}\t{fraction_to_time(sunset_time)}")


date_to_julian_century(datetime(2025, 12, 24), timezone_offset, time_past_midnight)
for hour in range (0, 23):
    sun_declin, eq_of_time = precalc(julian_century)
    azimuth, elevation = solar_position(latitude, sun_declin, eq_of_time, hour)
    print(f"{hour:02d}\t{azimuth}\t{elevation}")

Example output

Sunrise, noon, and sunset times in December 2025, Brno, Czech Republic

Date	Sunrise	Noon	Sunset
2025-12-01	07:27:21	11:42:40	15:57:59
2025-12-02	07:28:37	11:43:03	15:57:28
2025-12-03	07:29:52	11:43:27	15:57:01
2025-12-04	07:31:05	11:43:51	15:56:37
2025-12-05	07:32:15	11:44:16	15:56:16
2025-12-06	07:33:24	11:44:41	15:55:57
2025-12-07	07:34:31	11:45:07	15:55:42
2025-12-08	07:35:36	11:45:33	15:55:30
2025-12-09	07:36:38	11:46:00	15:55:21
2025-12-10	07:37:39	11:46:27	15:55:15
2025-12-11	07:38:37	11:46:54	15:55:12
2025-12-12	07:39:32	11:47:22	15:55:12
2025-12-13	07:40:25	11:47:50	15:55:16
2025-12-14	07:41:16	11:48:19	15:55:22
2025-12-15	07:42:04	11:48:48	15:55:31
2025-12-16	07:42:50	11:49:17	15:55:44
2025-12-17	07:43:33	11:49:46	15:56:00
2025-12-18	07:44:13	11:50:15	15:56:18
2025-12-19	07:44:50	11:50:45	15:56:40
2025-12-20	07:45:24	11:51:15	15:57:05
2025-12-21	07:45:56	11:51:44	15:57:33
2025-12-22	07:46:25	11:52:14	15:58:03
2025-12-23	07:46:51	11:52:44	15:58:37
2025-12-24	07:47:14	11:53:14	15:59:14
2025-12-25	07:47:34	11:53:43	15:59:53
2025-12-26	07:47:51	11:54:13	16:00:35
2025-12-27	07:48:04	11:54:42	16:01:20
2025-12-28	07:48:15	11:55:12	16:02:08
2025-12-29	07:48:23	11:55:41	16:02:58
2025-12-30	07:48:28	11:56:10	16:03:51
2025-12-31	07:48:30	11:56:38	16:04:47

Sun azimuth and elevation on 2025-12-24

Hour	Azimuth	Elevation
0	1.50	-63.78
1	30.97	-61.00
2	53.82	-54.36
3	70.66	-45.69
4	83.91	-36.14
5	95.33	-26.35
6	106.00	-16.73
7	116.61	-7.61
8	127.67	0.68
9	139.56	7.77
10	152.45	13.25
11	166.29	16.71
12	180.70	17.82
13	195.08	16.48
14	208.84	12.81
15	221.63	7.15
16	233.42	-0.07
17	244.42	-8.46
18	255.02	-17.64
19	265.72	-27.29
20	277.26	-37.07
21	290.77	-46.57
22	308.07	-55.11

MicroPython code using precomputed table

Code in MicroPython, may look like this and possibly run on ESP32 or Raspberry Pi Pico.

#!/usr/bin/micropython

import struct, math, time

declination = []
equation_of_time = []

def date_to_index(year, month, day, ref_year=2025, ref_month=12, ref_day=1):
    """Convert a date to table index (days since reference date)"""
    date_ts = time.mktime((year, month, day, 0, 0, 0, 0, 0))
    ref_ts = time.mktime((ref_year, ref_month, ref_day, 0, 0, 0, 0, 0))
    return (date_ts - ref_ts) // 86400

def index_to_date(index, ref_year=2025, ref_month=12, ref_day=1):
    """Convert table index to (year, month, day) tuple"""
    ref_ts = time.mktime((ref_year, ref_month, ref_day, 0, 0, 0, 0, 0, 0))
    target_ts = ref_ts + (index * 86400)
    date_tuple = time.localtime(target_ts)
    return (date_tuple[0], date_tuple[1], date_tuple[2])

def index_to_date_str(index, ref_year=2025, ref_month=12, ref_day=1):
    """Convert table index to date string YYYY-MM-DD"""
    year, month, day = index_to_date(index, ref_year, ref_month, ref_day)
    return f"{year}-{month:02d}-{day:02d}"

def fraction_to_time(fraction):
    hours   = int(fraction * 24)
    minutes = int((fraction * 24 - hours) * 60)
    seconds = int((fraction * 86400) % 60)
    return f"{hours:02d}:{minutes:02d}:{seconds:02d}"

def hour_angle_for_elevation(latitude, declination, elevation):
    lat_rad  = math.radians(latitude)
    dec_rad  = math.radians(declination)
    elev_rad = math.radians(elevation)
    
    # Spherical astronomy formula
    cos_h = (math.sin(elev_rad) - math.sin(lat_rad) * math.sin(dec_rad)) / (math.cos(lat_rad) * math.cos(dec_rad))
    
    # Check if sun reaches this elevation today
    if cos_h < -1.0:
        return None  # Sun always above this elevation (polar summer)
    elif cos_h > 1.0:
        return None  # Sun never reaches this elevation (polar winter)

    return math.degrees(math.acos(cos_h))

## Main program #############################################3

with open('sun_data.bin', 'rb') as f:
    count = struct.unpack('H', f.read(2))[0]
    for _ in range(count):
        decl, eot = struct.unpack('ff', f.read(8))
        declination.append(decl)
        equation_of_time.append(eot)

latitude  = 49.19528 
longitude = 16.60778
timezone_offset = 1   # UTC+1, CET

for i in range(0, 35):
    sun_declin = declination[i]
    eq_of_time = equation_of_time[i]
    ha_sunrise = hour_angle_for_elevation(latitude, sun_declin, -0.8333)
    solar_noon = (720 - 4 * longitude - eq_of_time + timezone_offset * 60) / 1440
    sunrise_time = solar_noon - ha_sunrise * 4 / 1440
    sunset_time  = solar_noon + ha_sunrise * 4 / 1440
    print(f"| {index_to_date_str(i)} | {fraction_to_time(sunrise_time)} | {fraction_to_time(solar_noon)} | {fraction_to_time(sunset_time)} |")

Quick Summary

• Sunrise/sunset calculation is surprisingly complex due to Earth’s elliptical orbit and axial tilt
• Accuracy within 1 minute requires polynomial approximations and careful time conversions
• Can pre-compute declination/equation-of-time tables for efficient embedded use
• Link to code snippets on GitHub