#
# Calculate solar position in equitorial coordinates (ascension, declination)
# then convert to horizontal coordinates (altitude, azimuth).
# Parameters L, g, lambda, and eps are appoximated using values chosen
# for accuracy in the era on either side of Jan 2000.
#
# Input:
#	Date - string with format "dd-mm-YYYY HH:MM"
#	Latitude - in degrees
#	Longitude - in degrees
#
# Output:
#	Function definitions
#		Azimuth(t) = Solar azimuth at time t
#		Altitude(t) = Solar altitude at time t
#	where t is time in seconds relative to local solar noon on Date.
#		sunrise, sunset (string containing local standard time)
#		sunlight (string containing time between sunrise, sunset)
#
# Local values:
	Minute = 60.
	Hour = 3600.
	Day = 86400.
	TimeFormat = "%d-%m-%Y %H:%M"
	Phi = Latitude
	set angle degrees
#
#     n	= local time in days since noon on 1 Jan 2000
#     L	= mean solar longitude in ecliptic coordinates
#  	= 280.460° + 0.9856474° * n
#     g	= orbital anomaly
#  	= 357.528° + 0.9856003° * n
# lamba	= solar longitude
#      	= L + 1.915° * sin(g) + 0.020° * sin(2g)
#   eps	= obliquity of the ecliptic
#	= 23.439° - 0.0000004° * n

n = (strptime(TimeFormat, Date) - strptime(TimeFormat, "01-01-2000 12:00")) / Day

L = n * 0.9856474 + 280.460
L = L - 360 * int(L / 360.)
if (L < 0) { L += 360. }

g = n * 0.9856003 + 357.528
g = g - 360 * int(g / 360.)
if (g < 0) { g += 360. }

lambda = L + 1.915 * sin(g) + 0.020 * sin(2*g)

eps = 23.439 - n * 0.0000004

# Equatorial coordinates
#     RAsc  = right ascension
#     Dec = declination

RAsc = atan2( cos(eps) * sin(lambda), cos(lambda) )
Dec = asin( sin(eps) * sin(lambda) )

# Length of day
daylength = 2.0 * acos( -sin(Dec)*sin(Phi) / cos(Phi) )
# Correction for displacement of this location from center of timezone
corr = Longitude - floor(Longitude/15.0) * 15.0 - 7.5
rise = (daylength/2.0 - corr) * Day/360.
set  = (daylength/2.0 + corr) * Day/360.

sunlight = strftime("%tH h %tM m", daylength * Day/360.)
sunrise = strftime("%tH:%tM", (Day/2)-rise)
sunset  = strftime("%tH:%tM", (Day/2)+set)

# Horizontal coordinates

h(t) = 360. * t / Day
Altitude(t) = asin( sin(Dec) * sin(Phi) + cos(Dec) * cos(Phi) * cos(h(t)) )

cosAzi(t)  = ( sin(Dec) * sin(Phi) - cos(Dec) * cos(h(t)) * sin(Phi) ) \
           / cos(Altitude(t))
sinAzi(t)  = ( -cos(Dec) * sin(h(t)) ) \
           / cos(Altitude(t))
Azimuth(t) = atan2( sinAzi(t), cosAzi(t) )