Here comes the sun
When does the sun rise and set? I have a fuzzy idea about how this works, here I want to make it a little bit less fuzzy.
#collapse
import datetime
from suntime import Sun, SunTimeException
latitude = 51.21
longitude = 21.01
sun = Sun(latitude, longitude)
# Get today's sunrise and sunset in UTC
today_sr = sun.get_sunrise_time()
today_ss = sun.get_sunset_time()
print('Today ({}) at Warsaw the sun raised at {} and get down at {} UTC'.
format(datetime.datetime.today().strftime('%Y-%m-%d'), today_sr.strftime('%H:%M'), today_ss.strftime('%H:%M')))
# On a special date in your machine's local time zone
abd = datetime.date(2014, 10, 3)
abd_sr = sun.get_local_sunrise_time(abd)
abd_ss = sun.get_local_sunset_time(abd)
print('On {} the sun at Warsaw raised at {} and get down at {}.'.
format(abd, abd_sr.strftime('%H:%M'), abd_ss.strftime('%H:%M')))
# Error handling (no sunset or sunrise on given location)
latitude = 87.55
longitude = 0.1
sun = Sun(latitude, longitude)
try:
abd_sr = sun.get_local_sunrise_time(abd)
abd_ss = sun.get_local_sunset_time(abd)
print('On {} at somewhere in the north the sun raised at {} and get down at {}.'.
format(abd, abd_sr.strftime('%H:%M'), abd_ss.strftime('%H:%M')))
except SunTimeException as e:
print("Error: {0}.".format(e))
I currently live here: 58°04'00.8"N 11°42'10.7"E
This one can be converted to Decimal degrees.
def longitude_str_to_num(s:str):
"""
Convert a longitude or latitude in the following format : 58°04'00.8"N
to degrees
"""
s1,s_=s.split('°')
s2,s_=s_.split("'")
s3,s_=s_.split('"')
D = float(s1)
M = float(s2)
S = float(s3)
D_deg = D + M/60 + S/3600
if (s_=='W') or (s_=='S'):
D_deg*=-1
return D_deg
latitude = longitude_str_to_num('58°04\'00.8"N')
latitude
longitude = longitude_str_to_num('11°42\'10.7"E')
longitude
sun = Sun(latitude, longitude)
sunrise_time = sun.get_sunrise_time()
sunrise_time
sunset_time = sun.get_sunset_time()
sunset_time
Getting these times in local time zone is somewhat messy:
import pytz
timezone = pytz.timezone(r'Europe/Stockholm')
sun_rise_time_sweden = sunrise_time.replace(tzinfo=pytz.utc).astimezone(timezone)
sun_rise_time_sweden
sun_set_time_sweden = sunset_time.replace(tzinfo=pytz.utc).astimezone(timezone)
sun_set_time_sweden
day_time = sun_set_time_sweden-sun_rise_time_sweden
str(day_time)
The Sun rise equation is written: $$\cos \omega_{\circ}=-\tan \phi \times \tan \delta$$ where $\phi$ is the latitude and $\delta$ is the earth inclination angle to the sun, which changes over the seasons.
The declination-angle can be calculated: $$ \delta=-23.45^{\circ} \times \cos \left(\frac{360}{365} \times(d+10)\right) $$ where $d=1$ at january 1st.
Putting it all together:
import numpy as np
import matplotlib.pyplot as plt
def calculate_declination_angle(day):
angle_deg=360/365*(day+10)
delta_deg = -23.45*np.cos(np.deg2rad(angle_deg))
delta = np.deg2rad(delta_deg)
return delta
#collapse
days=np.linspace(1,360,360)
delta=calculate_declination_angle(day=days)
fig,ax=plt.subplots()
ax.plot(days,np.rad2deg(delta));
ax.set_xlabel('days')
ax.set_ylabel('$\delta$ inclination angle [deg]');
def sun_hour_angle(latitude:float, day:int):
"""
param: latitude [deg]
param: day, january 1 --> day=1
"""
delta = calculate_declination_angle(day=day)
phi=np.deg2rad(latitude)
omega0 = np.arccos(-np.tan(phi)*np.tan(delta))
return omega0
def sun_rise_time(latitude:float, day:int):
"""
param: latitude [deg]
param: day, january 1 --> day=1
"""
omega0 = sun_hour_angle(latitude=latitude, day=day)
rise_time=12-omega0/np.deg2rad(15)
return rise_time
def sun_set_time(latitude:float, day:int):
"""
param: latitude [deg]
param: day, january 1 --> day=1
"""
omega0 = sun_hour_angle(latitude=latitude, day=day)
set_time=12+omega0/np.deg2rad(15)
return set_time
#collapse
rise_time=sun_rise_time(latitude=latitude, day=days)
set_time=sun_set_time(latitude=latitude, day=days)
fig,ax=plt.subplots()
ax.plot(days,rise_time, label='rise', color='y');
ax.plot(days,set_time, label='set', color='r');
ax.legend()
ax.set_xlabel('days')
ax.set_ylabel('time [hour]');
now = datetime.datetime.now()
day = (now-datetime.datetime(year=now.year, month=1, day=1)).days
rise_time = sun_rise_time(latitude=latitude, day=day)
rise_time = datetime.timedelta(hours=rise_time)
'%s' % rise_time
... so we are missing with about 15 minutes compared to the previous result. I think that the given time here is by definition in local time, since we substract from local time 12 o'clock as in the code above.
Wikipedia also says that there is a more advanced Sun rise equation:
$$ \cos \omega_{\circ}=\frac{\sin a-\sin \phi \times \sin \delta}{\cos \phi \times \cos \delta} $$ where $a=−0.83°$
Will that fix the 15 minutes?
#collapse
from numpy import sin,cos
def sun_hour_angle2(latitude:float, day:int):
"""
param: latitude [deg]
param: day, january 1 --> day=1
"""
delta = calculate_declination_angle(day=day)
phi=np.deg2rad(latitude)
a = np.deg2rad(-0.83)
omega0 = np.arccos((sin(a)-sin(phi)*sin(delta))/(cos(phi)*cos(delta)))
return omega0
def sun_rise_time2(latitude:float, day:int):
"""
param: latitude [deg]
param: day, january 1 --> day=1
"""
omega0 = sun_hour_angle2(latitude=latitude, day=day)
rise_time=12-omega0/np.deg2rad(15)
return rise_time
def sun_set_time2(latitude:float, day:int):
"""
param: latitude [deg]
param: day, january 1 --> day=1
"""
omega0 = sun_hour_angle2(latitude=latitude, day=day)
set_time=12+omega0/np.deg2rad(15)
return set_time
rise_time = sun_rise_time2(latitude=latitude, day=day)
rise_time = datetime.timedelta(hours=rise_time)
'%s' % rise_time
...a little bit better but still not perfect. There was also a simplification in the calculation of the inclination angle, will that do it?
def calculate_declination_angle2(day):
angle_deg=360/365*(day+10)
delta = np.arcsin(np.sin(np.deg2rad(-23.45))*np.cos(np.deg2rad(angle_deg)))
return delta
#collapse
delta=calculate_declination_angle(day=days)
delta2=calculate_declination_angle2(day=days)
fig,ax=plt.subplots()
ax.plot(days,np.rad2deg(delta),label='inclination linear');
ax.plot(days,np.rad2deg(delta2),'--', label='inlination nonlinear');
ax.set_xlabel('days')
ax.set_ylabel('$\delta$ inclination angle [deg]');
ax.legend()
#collapse
def sun_hour_angle3(latitude:float, day:int):
"""
param: latitude [deg]
param: day, january 1 --> day=1
"""
delta = calculate_declination_angle2(day=day)
phi=np.deg2rad(latitude)
#a = np.deg2rad(-0.83)
a = np.deg2rad(-1.07)
omega0 = np.arccos((sin(a)-sin(phi)*sin(delta))/(cos(phi)*cos(delta)))
return omega0
def sun_rise_time3(latitude:float, day:int):
"""
param: latitude [deg]
param: day, january 1 --> day=1
"""
omega0 = sun_hour_angle3(latitude=latitude, day=day)
rise_time=12-omega0/np.deg2rad(15)
return rise_time
def sun_set_time3(latitude:float, day:int):
"""
param: latitude [deg]
param: day, january 1 --> day=1
"""
omega0 = sun_hour_angle3(latitude=latitude, day=day)
set_time=12+omega0/np.deg2rad(15)
return set_time
rise_time = sun_rise_time3(latitude=latitude, day=day)
rise_time = datetime.timedelta(hours=rise_time)
'%s' % rise_time
yet a little bit better, but there is still some time missing... It is however so typical me to jump to conclusion instead of reading the description thuroughly. The sun rise and sun set should be calculated as:
$$ \begin{array}{l} J_{\text {rise}}=J_{\text {transit}}-\frac{\omega_{\circ}}{360^{\circ}} \\ J_{\text {set}}=J_{\text {transit}}+\frac{\omega_{\circ}}{360^{\circ}} \end{array} $$where $\mathrm{J}_{\text {rise }}$ is the actual Julian date of sunrise; J set is the actual Julian date of sunset.
$$ J_{\text {transit}}=2451545.0+J^{\star}+0.0053 \sin M-0.0069 \sin (2 \lambda) $$where: J transit is the Julian date for the local true solar transit (or solar noon). 2451545.0 is noon of the equivalent Julian year reference. $0.0053 \sin M-0.0069 \sin (2 \lambda)$ is a simplified version of the equation of time. The coefficients are fractional day minutes.
$$ \lambda=(M+C+180+102.9372) \bmod 360 $$where: $\lambda$ is the ecliptic longitude.
$$ C=1.9148 \sin (M)+0.0200 \sin (2 M)+0.0003 \sin (3 M) $$where: $\mathrm{C}$ is the Equation of the center value needed to calculate lambda (see next equation). 1.9148 is the coefficient of the Equation of the Center for the planet the observer is on (in this case, Earth)
$$ M=\left(357.5291+0.98560028 \times J^{\star}\right) \bmod 360 $$$$ J^{\star}=n-\frac{l_{w}}{360^{\circ}} $$where: $J^{\star}$ is an approximation of mean solar time at $n$ expressed as a Julian day with the day fraction. $l_{\omega}$ is the longitude west (west is negative, east is positive) of the observer on the Earth;
$$ n=J_{\text {date}}-2451545.0+0.0008 $$where: $n$ is the number of days since Jan 1 st, 2000 12:00 . $J_{\text {date}}$ is the Julian date;
def calulate_n(time:datetime.datetime):
start = datetime.datetime(year=2000, month=1, day=1, hour=12)
n = (time-start).total_seconds() / (60*60*24)
return n
def calculate_julian_date(time:datetime.datetime):
n = calulate_n(time=time)
return n + 2451545.0 - 0.0008
calculate_julian_date(datetime.datetime.today())
def calculate_mean_solar_time(n,longitude):
J_star = n+longitude/360
return J_star
def calculate_j_transit(time:datetime.datetime, longitude:float):
n = calulate_n(time=time)
j_star = calculate_mean_solar_time(n=n, longitude=longitude)
m = np.mod((357.5291+0.98560028*j_star),360)
c = 1.9148*np.sin(np.deg2rad(m))+0.0200*np.sin(2*np.deg2rad(m))+0.0003*np.sin(3*np.deg2rad(m))
lamda = np.mod(m+c+180+102.9372, 360)
j_transit = 2451545.0+j_star+0.0053*np.sin(np.deg2rad(m))-0.0069*np.sin(2*np.deg2rad(lamda))
return j_transit
j_transit = calculate_j_transit(time=datetime.datetime.today(), longitude=longitude)
j_transit - calculate_julian_date(datetime.datetime.today())
#collapse
def time_of_day(latitude:float, day:int):
return sun_set_time3(latitude=latitude, day=day) - sun_rise_time3(latitude=latitude, day=day)
time=time_of_day(latitude=latitude, day=days)
fig,ax=plt.subplots()
ax.plot(days,time, 'b-', label='rise');
ax.set_xlabel('days')
ax.set_ylabel('Length of day [hour]');
ax.grid(True)
#collapse
from matplotlib import ticker, cm
N=400
latitudes = np.linspace(0.1,89.9,N)
days_ = np.linspace(0,365,N)
L = np.tile(latitudes,(len(days_),1))
D = np.tile(days_,(N,1)).T
time=time_of_day(latitude=L, day=D)
fig,ax=plt.subplots()
fig.set_size_inches(15,10)
cs = ax.contourf(D,L,time, cmap=cm.gray, levels=48);
cb = fig.colorbar(cs)
ax.set_xlabel('days')
ax.set_ylabel('Latitude [deg]');
latitude = longitude_str_to_num('58°04\'00.8"N')
ax.plot([0,365],[latitude,latitude],'b-');
ax.plot(day,latitude,'ro');
ax.set_title('Hours per day at various days and latitudes');
The blue line is where I normally experience the seasons and at this very moment I'm in the red dot, listening to this: