#

# 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) )