Day 3 exercise: use Stull eqn 2.6 to find elevation for Dec. 22, Mar. 23 and Jun. 22 2014

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


    

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))