import math
from math import floor
import time
from datetime import date, timedelta
import datetime

#date = datetime.date(2015,5,16)
# set date and get LMT/clockTime (in hours)

#Date = datetime.datetime(2015,9,1,12,00)
date = datetime.datetime.now()
summerCorrection = 1 # BST or GMT


[input data]
    date = datetime.datetime.now()
    summerCorrection = 1 # BST or GMT
    longitude = 0 
    latitude = 51*(math.pi/180)    
    zoneMeridian = 0

clockTime = (int(date.strftime("%M")))/60 + int(date.strftime("%H"))+summerCorrection
# calculate declination angle and equation of time (based on time of year)
dayOfYear = date.timetuple().tm_yday
dayAngle = (dayOfYear*(360/365.25))*(math.pi/180)
declinationAngle = math.asin(0.3978*math.sin(dayAngle-1.4+0.0355*math.sin(dayAngle-0.0489)))
EOT_h = -0.128*math.sin((dayAngle-0.04887))-0.165*math.sin((2*dayAngle+0.3438))
# calculate local apparent time (solarTime) and equivilant solar hour angle
# account for position of longitude within Time Zone
solarTime = clockTime + ((longitude - zoneMeridian)/15) + EOT_h - summerCorrection 
hourAngle = (15*(solarTime-12))*(math.pi/180)
# set latitude and calculate solar altitude and azimuth
solarAltitude = math.asin((math.sin(latitude)*math.sin(declinationAngle))+((math.cos(latitude)*math.cos(declinationAngle)*math.cos(hourAngle))))
solarAzimuth = math.asin((math.sin(hourAngle)*math.cos(declinationAngle)/(math.cos(solarAltitude))))
#sunsetHourAngle = math.acos(-math.tan(latitude)*math.tan(declinationAngle))
#sunriseLAT = 12 - (sunsetHourAngle/(math.radians(15)))
#sunsetLAT = 12 + (sunsetHourAngle/(math.radians(15)))
#daylength = 2*(sunsetHourAngle/(math.radians(15)))
# daily extraerrestial irradiance
#solarConstant = 1367
#solarDistanceCorrection = 1.0 + 0.03344*math.cos(dayAngle-math.radians(2.80))
#dailyExtraterrestialIrradation = ((1/3600)/math.pi)*solarDistanceCorrection*solarConstant*((sunsetHourAngle*math.sin(latitude)*math.sin(declinationAngle))+(math.cos(latitude)*math.cos(declinationAngle)*math.sin(sunsetHourAngle)))


print('dayOfYear = %d' % dayOfYear)
print('hourAngle = %.3f' % ((hourAngle)*(180/math.pi)))
print('dayAngle = %.3f' % ((dayAngle)*(180/math.pi)))
print('EOT_h = %.3f' % (EOT_h*60))
print('summerCorrection = %.3f' % summerCorrection)
print('latitude = %.2f' % ((latitude)*(180/math.pi)))
print('longitude = %.3f' % longitude)
print('zoneMeridian = %.3f' % zoneMeridian)
print('solarTime = %.3f' % solarTime)
print('clockTime = %.3f' % clockTime)
print('declinationAngle = %.2f' % ((declinationAngle)*(180/math.pi)))
print('solarAltitude = %.3f' % ((solarAltitude)*(180/math.pi)))
print('solarAzimuth = %.3f' % ((solarAzimuth)*(180/math.pi)))
print('sunsetHourAngle = %.2f' % ((sunsetHourAngle)*(180/math.pi)))
print('sunriseLAT = %.2f' % sunriseLAT)
print('sunsetLAT = %.2f' % sunsetLAT)
print('daylength = %.2f' % daylength)
#print('solarConstant = %.2f' % solarConstant)
#print('solarDistanceCorrection = %.2f' % solarDistanceCorrection)
#print('dailyExtraterrestialIrradation = %.4f' % dailyExtraterrestialIrradation)
#print('ExtraterrestialIrradiance = %.4f' % (dailyExtraterrestialIrradation*3600))