this function uses the python datetime module to correctly find the length of the year for normal or leap years
In [0]:
#make all division floating point
from __future__ import division
import datetime
import pytz
from math import asin,sin,cos,pi
def find_deltas(the_date):
"""given a python datetime object (local time)
find the solar declination angle in degrees
using Stull equation 2.5
input: datetime object with or without timezone
output: solar declination angle in degrees
"""
the_year=the_date.year
#
# find the length of the year in days by subtracting
# two datetimes exactly 1 year apart -- jan 1, 0 hours, 0 minutes, 0 seconds
#
year_start=datetime.datetime(the_year,1,1,0,0,0)
year_end=datetime.datetime(the_year+1,1,1,0,0,0)
year_length=(year_end - year_start).days
phir=23.44 #axis tilt in degrees from stull
#number of days since the new year
the_day=(the_date - year_start).days
solstice=172 #stull
fraction=(the_day - solstice)/year_length
deltas=phir*cos(2*pi*fraction)
return deltas
use the timezone module pytz to convert local vancouver time to UTC, correctly accounting for daylight savings time
In [0]:
def find_elevation(the_date):
"""find the solar elevation for Vancouver in degrees given a python
datetime object without a timezone representing
vancouver local time, using Stull eqn. 2.6
input: datetime object (no timezone)
output: solar elevation in degrees
"""
deg2rad=pi/180.
rad2deg=1./deg2rad
deltas=find_deltas(the_date)
deltas=deltas*deg2rad
#stull p. 32 gives lat and lon for vancouver
phi=49.25*deg2rad #vancouver latitude deg N
lambda_e = 123.1*deg2rad #vancouver longitude, deg W
pacific=pytz.timezone('US/Pacific')
utc=pytz.timezone('UTC')
#
# give the time a pacific timezone
#
vancouver=pacific.localize(the_date)
#convert to utc
vancouver=vancouver.astimezone(utc)
t_utc=vancouver.hour
#stull eqn 2.6
sin_psi=sin(phi)*sin(deltas) - cos(phi)*cos(deltas)*cos(2*pi*t_utc/24. - lambda_e)
elevation=asin(sin_psi)*rad2deg
#write 0 if under the horizon
if elevation < 0:
elevation=0.
return elevation
In [0]:
print("note that on 6/22 the highest elevation occurs at 1 pm, because of DST")
print("{:5s} {:5s} {:5s} {:5s}".format('hour','12/22','3/23','6/22') )
print("="*30)
for hour in range(3,23):
dec22=datetime.datetime(2014,12,22,hour,0,0)
mar23=datetime.datetime(2014,3,23,hour,0,0)
jun22=datetime.datetime(2014,6,22,hour,0,0)
elevation=[find_elevation(item) for item in [dec22,mar23,jun22]]
#right justify the columns and make them all 5 characters wide
# see http://clouds.eos.ubc.ca/~phil/djpine_python/Book/_build/html/chap4/chap4_io.html?highlight=format
# see http://clouds.eos.ubc.ca/~phil/djpine_python/Book/_build/html/chap7/chap7_funcs.html#user-defined-functions
# for an explanation of *elevation
print("{:>5d} {:>5.2f} {:>5.2f} {:>5.2f}".format(hour,*elevation))