Finding the position of the sun

sun_position.py

def sun_position ( date, hour, longitude, latitude ):
    """
    Calculates the position of the sun given a position and time.
    Basically, all you need to know is here: 
    http://answers.google.com/answers/threadview/id/782886.html
    """
    import time
    from math import sin, cos, degrees, radians, acos

    doy = int ( time.strftime( "%j", time.strptime( date, "%Y-%m-%d" ) ) )
    latitude = np.deg2rad ( latitude )
    longitude = np.deg2rad ( longitude )
    ( hh, mm )  = hour.split(":")
    h = float(hh)
    h = h + float(mm)/60.
    ## this is the fractional year
    gdeg = (360/365.25)*(doy+ h/24)
    g = radians(gdeg)

    ## SOLAR DECLINATION
    ## Follows the calculation of the declination of the sun, again you an use a table or a formula
    ##"Table of the Declination of the Sun": http://www.wsanford.com/~wsanford/exo/sundials/DEC_Sun.html
    ## Or use the following formula:
    D = 0.396372-22.91327*cos(g)+4.02543*sin(g)-0.387205*cos(2*g) + \
    0.051967*sin(2*g)-0.154527*cos(3*g) + 0.084798*sin(3*g)


    ##Now calculate the TIME CORRECTION for solar angle:
    TC = 0.004297+0.107029*cos(g)-1.837877*sin(g)-0.837378*cos(2*g) - \
    2.340475*sin(2*g)


    ## Now we can calculate the Solar Hour Angle (SHA)
    SHA = (h-12)*15 + longitude + TC
    if SHA > 180:
        SHA = SHA - 360
    elif SHA< -180:
        SHA = SHA + 360

    ##Now we can calculate the Sun Zenith Angle (SZA):
    cosSZA =  sin(radians(latitude))*sin(radians(D)) + \
            cos(radians(latitude))*cos(radians(D))*cos(radians(SHA))
    SZA = degrees(acos(cosSZA))
    altitude = 90 - SZA

    ##To finish we will calculate the Azimuth Angle (AZ):
    cosAZ = (sin(radians(D)) - sin(radians(latitude)) * cos(radians(SZA))) /\
            (cos(radians(latitude))*sin(radians(SZA)))
    AZ = degrees(acos(cosAZ))


    return (SZA, AZ)