1#
2# Calculate solar position in equatorial coordinates (ascension, declination)
3# then convert to horizontal coordinates (altitude, azimuth).
4# Parameters L, g, lambda, and eps are appoximated using values chosen
5# for accuracy in the era on either side of Jan 2000.
6#
7# Input:
8#	Date - string with format "dd-mm-YYYY HH:MM"
9#	Latitude - in degrees
10#	Longitude - in degrees
11#
12# Output:
13#	Function definitions
14#		Azimuth(t) = Solar azimuth at time t
15#		Altitude(t) = Solar altitude at time t
16#	where t is time in seconds relative to local solar noon on Date.
17#		sunrise, sunset (string containing local standard time)
18#		sunlight (string containing time between sunrise, sunset)
19#
20# Local values:
21	Minute = 60.
22	Hour = 3600.
23	Day = 86400.
24	TimeFormat = "%d-%m-%Y %H:%M"
25	Phi = Latitude
26	set angle degrees
27#
28#     n	= local time in days since noon on 1 Jan 2000
29#     L	= mean solar longitude in ecliptic coordinates
30#  	= 280.460° + 0.9856474° * n
31#     g	= orbital anomaly
32#  	= 357.528° + 0.9856003° * n
33# lamba	= solar longitude
34#      	= L + 1.915° * sin(g) + 0.020° * sin(2g)
35#   eps	= obliquity of the ecliptic
36#	= 23.439° - 0.0000004° * n
37
38n = (strptime(TimeFormat, Date) - strptime(TimeFormat, "01-01-2000 12:00")) / Day
39
40L = n * 0.9856474 + 280.460
41L = L - 360 * int(L / 360.)
42if (L < 0) { L += 360. }
43
44g = n * 0.9856003 + 357.528
45g = g - 360 * int(g / 360.)
46if (g < 0) { g += 360. }
47
48lambda = L + 1.915 * sin(g) + 0.020 * sin(2*g)
49
50eps = 23.439 - n * 0.0000004
51
52# Equatorial coordinates
53#     RAsc  = right ascension
54#     Dec = declination
55
56RAsc = atan2( cos(eps) * sin(lambda), cos(lambda) )
57Dec = asin( sin(eps) * sin(lambda) )
58
59# Length of day
60daylength = 2.0 * acos( -sin(Dec)*sin(Phi) / cos(Phi) )
61# Correction for displacement of this location from center of timezone
62corr = Longitude - floor(Longitude/15.0) * 15.0 - 7.5
63rise = (daylength/2.0 - corr) * Day/360.
64set  = (daylength/2.0 + corr) * Day/360.
65
66sunlight = strftime("%tH h %tM m", daylength * Day/360.)
67sunrise = strftime("%tH:%tM", (Day/2)-rise)
68sunset  = strftime("%tH:%tM", (Day/2)+set)
69
70# Horizontal coordinates
71
72h(t) = 360. * t / Day
73Altitude(t) = asin( sin(Dec) * sin(Phi) + cos(Dec) * cos(Phi) * cos(h(t)) )
74
75cosAzi(t)  = ( sin(Dec) * sin(Phi) - cos(Dec) * cos(h(t)) * sin(Phi) ) \
76           / cos(Altitude(t))
77sinAzi(t)  = ( -cos(Dec) * sin(h(t)) ) \
78           / cos(Altitude(t))
79Azimuth(t) = atan2( sinAzi(t), cosAzi(t) )
80
81