#!/usr/bin/env python # coding: utf-8 # # 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. # # - toc: true # - branch: master # - badges: true # - comments: true # - categories: [math] # - image: https://exitoina.uol.com.br/media/_versions/beatlessss_widelg.jpg # - hide: false # - search_exclude: true # - metadata_key1: metadata_value1 # - metadata_key2: metadata_value2 # ## Using suntime package # Eventhough this could be cheating, it's a good starting point to get a feel... # This example is taken straight from the documentation: # In[1]: #hide import warnings warnings.filterwarnings("ignore") # In[2]: #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](https://en.wikipedia.org/wiki/Decimal_degrees). # In[3]: 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 # In[4]: latitude = longitude_str_to_num('58°04\'00.8"N') latitude # In[5]: longitude = longitude_str_to_num('11°42\'10.7"E') longitude # In[6]: sun = Sun(latitude, longitude) # In[7]: sunrise_time = sun.get_sunrise_time() sunrise_time # In[8]: sunset_time = sun.get_sunset_time() sunset_time # Getting these times in local time zone is somewhat messy: # In[9]: import pytz timezone = pytz.timezone(r'Europe/Stockholm') sun_rise_time_sweden = sunrise_time.replace(tzinfo=pytz.utc).astimezone(timezone) sun_rise_time_sweden # In[10]: sun_set_time_sweden = sunset_time.replace(tzinfo=pytz.utc).astimezone(timezone) sun_set_time_sweden # In[11]: day_time = sun_set_time_sweden-sun_rise_time_sweden str(day_time) # ## Own implementation # The [Sun rise equation](https://en.wikipedia.org/wiki/Sunrise_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](https://www.pveducation.org/pvcdrom/properties-of-sunlight/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: # In[12]: 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 # In[13]: #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]'); # In[14]: 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 # In[15]: #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]'); # In[16]: 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](https://en.wikipedia.org/wiki/Sunrise_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? # In[108]: #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 # In[109]: 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? # In[110]: 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 # In[111]: #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() # In[114]: #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 # In[115]: 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; # In[132]: 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 # In[142]: calculate_julian_date(datetime.datetime.today()) # In[150]: 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 # In[151]: j_transit = calculate_j_transit(time=datetime.datetime.today(), longitude=longitude) j_transit - calculate_julian_date(datetime.datetime.today()) # In[43]: #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) # In[50]: #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: # # # In[ ]: