Sunset/rise calculation

I'm using the suntime class that is based on: "Calculate sunrise and sunset times for a given GPS coordinate within PostgreSQL". I'm struggling with the following:
At the end of the function get_sun_timedelta, force range is applied to the UT. However, I checked the data from the class and compared it to the sunset / sunrise events from https://aa.usno.navy.mil/ for several locations over a year. I noticed that by force ranging UT, some events are missed (especially if 2 sunsets/rises happen in one day). This can be prevented by not using force range. BUT without force range, errors occur the first day of the year or the last. Somehow they are not calculated properly. I don't quite understand how to prevent the skipping of the first and/or last day of the year and accuratley calculate it without using force range which causes errors in other occasions. Besides that the class works great and is exactly what I need. Maybe anyone knows how to fix this? I attached the calculation between suntimes and USNO without forced range.

Addendum: Two sunsets can happen in a day if the time is measured in UTC and the days/nights are short: e.g. sunrise at 00:00:00, sunset at 19:00:00, sunrise at 23:50:00. Somehow the indentation is not applied properly when copy pasting the code here. But the code I posted is from the package "suntime". You can simply install it "pip install suntime" and run and reproduce my error case from there by changing the line after "UT = round(UT,2)" and enabling the force_range for all UT calculations or none. The default is, that force range is only applied when sunrises are calculated which still leads to the error I'm presenting here which happens for both sunsets and sunrises. I added a table showing the input I use when calling the function to calculate sunsets and the output I get compared to the USNO data representing the calculated events. I marked the events that are missed in the calculation with force range or without force range. I added images showing the difference over a year between calculation and USNO data with and without the force_range function when calculating UT.

WITH force_range in UT calculation:

Without force_range in UT calculation:

import math
import warnings
from datetime import datetime, timedelta, time, timezone
import pandas as pd
from collections import defaultdict
import numpy as np

TO_RAD = math.pi/180.0

class SunTimeException(Exception):

def __init__(self, message):
    super(SunTimeException, self).__init__(message)


class Sun:
"""
Approximated calculation of sunrise and sunset datetimes. Adapted from:
https://stackoverflow.com/questions/19615350/calculate-sunrise-and-sunset-times-for-a-given-gps-coordinate-within-postgresql
"""
def __init__(self, lat, lon):
    self._lat = lat
    self._lon = lon

    self.lngHour = self._lon / 15

def get_sunrise_time(self, at_date=datetime.now(), time_zone=timezone.utc):
    """
    :param at_date: Reference date. datetime.now() if not provided.
    :param time_zone: pytz object with .tzinfo() or None
    :return: sunrise datetime.
    :raises: SunTimeException when there is no sunrise and sunset on given location and date.
    """
    time_delta = self.get_UT(at_date, time_zone=time_zone, is_rise_time=True)
    if time_delta is None:
        raise SunTimeException('The sun never rises on this location (on the specified date)')
    elif isinstance(time_delta, str):
        return time_delta
    else:
        event = datetime.combine(at_date, time(tzinfo=time_zone)) + time_delta
        return event

def get_sunset_time(self, at_date=datetime.now(), time_zone=timezone.utc):
    """
    Calculate the sunset time for given date.
    :param at_date: Reference date. datetime.now() if not provided.
    :param time_zone: pytz object with .tzinfo() or None
    :return: sunset datetime.
    :raises: SunTimeException when there is no sunrise and sunset on given location and date.
    """
    time_delta = self.get_UT(at_date, time_zone=time_zone, is_rise_time=False)
    if time_delta is None:
        raise SunTimeException('The sun never rises on this location (on the specified date)')
    elif isinstance(time_delta, str):
        return time_delta
    else:
        event = datetime.combine(at_date, time(tzinfo=time_zone)) + time_delta
        return event
    
def get_UT(self, at_date, time_zone, is_rise_time=True, zenith=90.8):
    """
    Calculate sunrise or sunset date.
    :param at_date: Reference date
    :param time_zone: pytz object with .tzinfo() or None
    :param is_rise_time: True if you want to calculate sunrise time.
    :param zenith: Sun reference zenith
    :return: timedelta showing hour, minute, and second of sunrise or sunset
    """

    # If not set get local timezone from datetime
    if time_zone is None:
        time_zone = datetime.now().tzinfo

    # 1. first get the day of the year
    N = at_date.timetuple().tm_yday
    # if at_date.date() == datetime(at_date.year, 1,1).date():
    #     N = 0 #TODO: BAU DAS MIT DEN 3 OPTIONEN HIER EIN IN DER KLASSE. TESTE DAS.
    # 2. convert the longitude to hour value and calculate an approximate time
    if is_rise_time:
        t = N + ((6 - self.lngHour) / 24)
    else:   # sunset
        t = N + ((18 - self.lngHour) / 24)
    
    # 3a. calculate the Sun's mean anomaly
    M = (0.9856 * t) - 3.289

    # 3b. calculate the Sun's true longitude
    L = M + (1.916 * math.sin(TO_RAD*M)) + (0.020 * math.sin(TO_RAD * 2 * M)) + 282.634
    L = self._force_range(L, 360)   # NOTE: L adjusted into the range [0,360)

    # 4a. calculate the Sun's declination
    sinDec = 0.39782 * math.sin(TO_RAD*L)
    cosDec = math.cos(math.asin(sinDec))

    # 4b. calculate the Sun's local hour angle
    cosH = (math.cos(TO_RAD*zenith) - (sinDec * math.sin(TO_RAD*self._lat))) / (cosDec * math.cos(TO_RAD*self._lat))

    if cosH > 1:
        return "constant_eclipse"     #TODO: -->eclipse # The sun never rises on this location (on the specified date)
    if cosH < -1:
        return "constant_sun"     #TODO: -->daylight The sun never sets on this location (on the specified date)

    # 4c. finish calculating H and convert into hours
    if is_rise_time:
        H = 360 - (1/TO_RAD) * math.acos(cosH)
    else:   # setting
        H = (1/TO_RAD) * math.acos(cosH)
    H = H / 15

    # 5a. calculate the Sun's right ascension
    RA = (1/TO_RAD) * math.atan(0.91764 * math.tan(TO_RAD*L))
    RA = self._force_range(RA, 360)     # NOTE: RA adjusted into the range [0,360)

    # 5b. right ascension value needs to be in the same quadrant as L
    Lquadrant = (math.floor(L/90)) * 90
    RAquadrant = (math.floor(RA/90)) * 90
    RA = RA + (Lquadrant - RAquadrant)

    # 5c. right ascension value needs to be converted into hours
    RA = RA / 15

    # 6. calculate local mean time of rising/setting
    T = H + RA - (0.06571 * t) - 6.622

    # 7a. adjust back to UTC
    UT = T - self.lngHour

    if time_zone:
        # 7b. adjust back to local time
        UT += time_zone.utcoffset(at_date).total_seconds() / 3600

    # 7c. rounding and impose range bounds
    UT = round(UT, 2)
    # NOTE: UT with force range : miss events, without: problems for 01.01/31.12
    # UT = self._force_range(UT, 24)
    # NOTE: original code, force range only for rise time
    # if is_rise_time:
    #     UT = self._force_range(UT, 24)
    
    # 8. return timedelta
    return timedelta(hours=UT)
    #return UT

@staticmethod
def _force_range(v, max):
    # force v to be >= 0 and < max
    if v < 0:
        return v + max
    elif v >= max:
        return v - max
    return v

I'll treat this as a debugging problem, and not discuss the virtues of the algorithm itself or the lack thereof. Still, relevant to the solution of the problem is the fact that the algorithm is limited in the fact that it computes a single value for the right ascension and for the declination of the Sun in an interval of 24h; it only changes with at_date at 24h/0h UTC. This (and the other corners cut) make the approximation very poor when and where the daily variations of sunrise or sunset times are very high, like for the first (or last) sunrises (or sunsets) in the polar areas, at the edge of polar nights or polar days.

It's quite clear that force-ranging UT to 24h will result in loss of information. However, what can and should be force-ranged is the "local mean time", T. By listing the values computed for T during a year, we'll see that on the March equinox, when the Sun's right ascension passes from 23.9h to 0h, T (computed as H + RA - ....) has a false jump of -24h that is directly impacting the UT, introducing false features into the results:

        # .................................       
        # 6. calculate local mean time of rising/setting
        T = H + RA - (0.06571 * t) - 6.622
        T = self._force_range(T, 24)

Without force-ranging UT we'll get UT > 24, that will be correctly interpreted (see the values of t) as the sunset time for the next day (the date in the result will be the input date +1day).

The question is then how to work with the initial moment at_date in two scenarios: (1) when tabulating all the rise and set times in a larger interval, like a year, and (2) when attempting to find the rise and set times or the Sun for a certain day.

With these, the solution for tabulating the sunset/sunrise times for a long interval (e.g., a year) could be:

• in the input dates, start with one day prior to the intended start date, and end with one day after the end date
• save data only based on the resulting datetime, i.e., get the date from the result and completely disregard the input date,
• discard those values that are outside the target interval

Since the functions get_sunrise_time and get_sunset_time take a timezone second argument, I'll assume the first argument, at_time is a datetime value with timezone 0 (tzinfo=timezone.utc).

You haven't posted your "main" program, so the following example will use simple lists and save data as strings for a "visual" comparison with the reference values:

sun = Sun(-30, -100)
tz = timezone(timedelta(hours=0))

sun_table = []
current_date_iso = None

date0 = datetime(2024, 1, 1, tzinfo=timezone.utc)

date = date0
n_days = 0
# Initialize the list and get as a bonus the number of days in the year 2024
while date.year == date0.year:
    sun_table.append([date.date().isoformat(), [], []]) # date, sunrise(s), sunset(s)
    date += timedelta(days=1)
    n_days += 1

# start one day prior to the intended start day
date = date0 - timedelta(days=1)
idx1 = 0
idx2 = 0
for _ in range(n_days + 2):  # go one more day after the intended end
    sunrise = sun.get_sunrise_time(date, tz)
    if isinstance(sunrise, datetime) and sunrise.year == date0.year: # discard results outside the target interval
        date_iso = sunrise.date().isoformat()
        while sun_table[idx1][0] != date_iso:
            idx1 += 1
        # or, less efficiently:
        # idx1 = next(i for i,v in enumerate(sun_table) if v[0] == date_iso)
        sun_table[idx1][1].append(sunrise.time().isoformat())

    sunset = sun.get_sunset_time(date, tz)
    if isinstance(sunset, datetime) and sunset.year == date0.year:
        date_iso = sunset.date().isoformat()
        while sun_table[idx2][0] != date_iso:
            idx2 += 1
        sun_table[idx2][2].append(sunset.time().isoformat())
    date += timedelta(days=1)

print(*sun_table, sep='\n')

playground

I don't know exactly how your plots were calculated. However, what is still possible to happen is this: for the day when the sunset changes from 00:00 to 23:59 (the day with two sunsets), because of the imprecision of the algorithm, you might get, let's say, 23:58 instead of 23:59 and 23:59 instead of 00:00, so the date when this shift 0->24 is calculated will be +1day from the real date. If we plot the actual difference for each day, this results in a huge spike of almost 24 h, but the actual error of the algorithm is just 1 min.

I also created plots of the errors, but reducing them modulo 24h. I got results of the following type (I can't post the code as it's an ugly mess of local scripts and not even in python). I used them to compare the results of various choices for the algorithm and validate the ones I described above (i.e, about force-ranging T and not force-ranging UT).

If you want to find the sunset or sunrise times for a certain day, you'll have to do exactly the same thing: compute the values for 3 days: the previous one, the targeted one and the next day, and collect those results that are in the target day.

def get_sunsets_for_date(date0, tz):
    sunsets = []
    date = date0 - timedelta(days=1)
    for _ in range(3):
        sunset = sun.get_sunset_time(date, tz)
        if isinstance(sunset, datetime) and sunset.date() == date0.date():
            sunsets.append(sunset.time().isoformat())
        date += timedelta(days=1)
    return sunsets

playground