def CalculateSunPosition(date, observerPoint):

    if debug == True:
        arcpy.AddMessage("Determining sun position for date: " + str(date))

    twoPi = math.pi * 2
    lat = observerPoint.Y
    lon = observerPoint.X

    if debug == True:
        arcpy.AddMessage("lat: " + str(lat))
        arcpy.AddMessage("lon: " + str(lon))

    # Determine the the day in the year (i.e. 123rd, 345th, 366th, etc...) on which the requested date falls
    timetuple = date.timetuple()
    tempHour = timetuple.tm_hour + (timetuple.tm_min / 60) + (timetuple.tm_sec/3600)
    julianDate = date_to_julian_day(date)
    time = julianDate - 2451545.0

    # Ecliptic coordinates

    # Mean Longitude
    mnlong = 280.460 + .9856474 * time
    mnlong = mnlong % 360
    if mnlong < 0:
        mnlong += 360

    # Mean Anomaly
    mnanom = 357.528 + .9856003 * time
    mnanom = mnanom % 360
    if mnanom < 0:
        mnanom += 360
    mnanom = math.radians(mnanom)

    # Ecliptic longitude and obliquity of ecliptic
    eclong = mnlong + 1.915 * math.sin(mnanom) + 0.020 * math.sin(2*mnanom)
    eclong = eclong % 360
    if eclong<0:
        eclong += 360
    oblqec = 23.439 - 0.0000004 * time
    eclong = math.radians(eclong)
    oblqec = math.radians(oblqec)

    # Celestial Coordinates
    # Right ascension and declination
    num = math.cos(oblqec) * math.sin(eclong)
    den = math.cos(eclong)
    ra = math.atan(num / den)
    if den < 0:
        ra += math.pi

    if den >= 0 and num < 0:
        ra += twoPi
    dec = math.asin(math.sin(oblqec) * math.sin(eclong))

    # Local coordinates
    # Greenwich mean sidereal time
    gmst = 6.697375 + .0657098242 * time + tempHour
    gmst = gmst % 24
    if gmst < 0:
        gmst += 24

    # Local mean sidereal time
    lmst = gmst + lon / 15
    lmst = lmst % 24
    if lmst < 0:
        lmst += 24
    lmst = lmst * 15.0
    lmst = math.radians(lmst)

    # Hour angle
    ha = lmst - ra
    if ha < (-math.pi):
        ha += twoPi
    if ha > math.pi:
        ha -= twoPi

    # Latitude to Radians
    latRad = math.radians(lat)

    # Solar zenith angle amd Azimuth
    zenithAngle = math.acos(math.sin(latRad) * math.sin(dec) + math.cos(latRad) * math.cos(dec) * math.cos(ha))
    az = math.acos(((math.sin(latRad) * math.cos(zenithAngle)) - math.sin(dec)) / (math.cos(latRad) * math.sin(zenithAngle)))

    # Elevation
    el = math.asin(math.sin(dec) * math.sin(latRad) + math.cos(dec) * math.cos(latRad) * math.cos(ha))

    if debug == True:
        arcpy.AddMessage("dec: " + str(dec))
        arcpy.AddMessage("latRad: " + str(latRad))
        arcpy.AddMessage("ha: " + str(ha))

    # For logic and names, see Spencer, J.W. 1989 Solar Energy.  42(4):353
    cosAzPos = ( 0 <= (math.sin(dec) - math.sin(el) * math.sin(lat)) )
    sinAzNeg = ( math.sin(az) < 0 )

    if debug == True:
        arcpy.AddMessage("Calculated Elevation Radians: " + str(el))
    if debug == True:
        arcpy.AddMessage("Calculated Azimuth Radians: " + str(az))

    el = math.degrees(el)
    az = math.degrees(az)

    if ha > 0:
        az = az + 180
    else:
        az = 540 - az
    az = az % 360

    if debug == True:
        arcpy.AddMessage("Calculated Elevation Degrees: " + str(el))
    if debug == True:
        arcpy.AddMessage("Calculated Azimuth Degrees: " + str(az))

    sp = SunPosition(date,observerPoint,az,el)
    return sp

# Check out any necessary licenses
dataset/ETHPy150Open Esri/solutions-geoprocessing-toolbox/visibility/toolboxes/scripts/spa.py/CalculateSunPosition/original