Reuse the day 3 elevation notebook to get the solar elevation for Wallace and Hobbs Figure 4.24 -- data were taken on December 12, 1970 in Tucson, Arizona

The Wallace and Hobbs figure shows measurements of the radiance at solar zenith angles with sec(theta) ranging from 2 to 10 -- what time of day were these measurements made?

Step 1: find solar declination angle $\delta_s$

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 as dt
import pytz
from math import asin,sin,cos,pi
import numpy as np

deg2rad=pi/180.
rad2deg=1./deg2rad

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=dt.datetime(the_year,1,1,0,0,0)
    year_end=dt.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
    #now find the Julian day of the solstice
    #for fun, get the exact solstice from the pyephem package
    #http://shallowsky.com/blog/science/astro/calculating-solstice.html
    try:
        #
        # use the exact solstice for 1970 if the ephem module is installed
        #
        import ephem
        solstice=ephem.next_solstice('1970/6/1').datetime()
    except ImportError:
        #print("ephem package isn't available, using calendar solstice")
        solstice=dt.datetime(1970,6,21,0,0,0)
    #print('solstice set to: ',solstice)
    #print "ephem solstice is: ",solstice
    jan1=dt.datetime(1970,1,1,0,0,0)
    solstice_day=(solstice - jan1).days
    #print('solstice has {} days'.format(solstice_day))
    fraction=(the_day - solstice_day)/year_length
    deltas=phir*cos(2*pi*fraction)
    return deltas

Step 2: find the elevation using equation 2.6

use the timezone module pytz to convert Tucson time to UTC, correctly accounting for daylight savings time



In [0]:

    
def find_elevation(the_date):
    """find the solar elevation for Tucson 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
    """
    deltas=find_deltas(the_date)
    deltas=deltas*deg2rad
    #stull p. 32 gives lat and lon for Tucson
    phi=32.221*deg2rad #Tucson latitude deg N
    lambda_e = 110.926*deg2rad #Tucson longitude, deg W
    mountain=pytz.timezone('US/Mountain')
    utc=pytz.timezone('UTC')
    #
    # give the time a pacific timezone
    #
    tucson=mountain.localize(the_date)
    #convert to utc
    tucson=tucson.astimezone(utc)
    #
    #  turn minutes into fractions of an hour
    #
    t_utc=tucson.hour + tucson.minute/60.
    #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

Step 3: print cos(theta) and secant(theta) for the sun for Dec. 12, 1970 at 15 minute intervals



In [0]:

    
print("{:5s}  {:5s}  {:10s}  {:10s}".format('hour','zenith','cos(theta)','sec(theta)') )
print("="*40)
for hour in np.arange(6,19):
    for minute in np.arange(0,60,15):
        dec12=datetime.datetime(1970,12,12,hour,minute,0)   
        elevation=find_elevation(dec12) 
        #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

        #zenith angle is measured from straight up
        zenith=90. - elevation
        the_cos=cos(deg2rad*zenith)
        if the_cos < 1.e-5:
            the_cos=0.
            the_sec=0.
        else:
            the_sec=1./the_cos
        print("{:>02d}:{:>02d}   {:>5.2f}   {:>5.2f}    {:>5.2f}".format(hour,minute,zenith,the_cos,the_sec))



In [0]: