[Python] 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 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. Calculated data compared to USNO data

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):
"""
aram at_date: Reference date. datetime.now() if not provided.
aram 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.
aram at_date: Reference date. datetime.now() if not provided.
aram 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.
aram at_date: Reference date
aram time_zone: pytz object with .tzinfo() or None
aram is_rise_time: True if you want to calculate sunrise time.
aram 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

Continue reading...