REM Based on a program by Roger W. Sinnott
REM lat = latitude (degrees north)
REM long = longitude (degrees east)
REM Z = time zone (hours offset from GMT)
REM Y = year
REM M = month (1-12)
REM D = day (1-31)
REM tr = sunrise time (hours), -1 indicates no sunrise
REM ts = sunset time (hours), -1 indicates no sunset
REM ar = sunrise azimuth (degrees)
REM as = sunset azimuth (degrees)
DEF PROCsunriseset(lat, long, Z, Y, M, D, RETURN tr, RETURN ts, RETURN ar, RETURN as)
LOCAL A0, A1, C, D0, D1, D2, DD, H0, H2, HD, lst, N, S, T, T0, Z0
T = FN_mjd(D,M,Y) - (FN_mjd(1,1,2000) + 0.5) : REM Days from noon on 1st Jan 2000
T0 = T / 36525 : REM Centuries from noon on 1st Jan 2000
Z0 = -Z / 24 : REM Time zone offset in days
lst = 24110.5 + 8640184.812999999 * T0 + 86636.6 * Z0 + 240 * long : REM seconds
lst = RAD((lst MOD 86400) / 240) : REM Local sidereal time (radians)
PROCsun(T + Z0, A0, D0) : REM Sun's position 'today'
PROCsun(T + Z0 + 1, A1, D1) : REM Sun's position 'tomorrow'
IF A1 < A0 THEN A1 += PI * 2
H0 = lst - A0
HD = 15 * RAD(1.0027379) - (A1 - A0) / 24
DD = (D1 - D0) / 24
S = SIN(RAD(lat))
C = COS(RAD(lat))
rise = -1
set = -1
FOR N = 0 TO 23
D2 = D0 + DD : H2 = H0 + HD
PROCtestevent(C, S, N, H0, H2, D0, D2, tr, ts, ar, as)
D0 = D2 : H0 = H2
NEXT
ENDPROC
REM Test an hour window for an event
REM H0,H2 = start and end Hour angle
REM D0,D2 = start and end Declination
DEF PROCtestevent(C, S, N, H0, H2, D0, D2, RETURN tr, RETURN ts, RETURN ar, RETURN as)
LOCAL A, B, D, E, H1, D1, R0, R1, R2, U, V, W, Z
Z = COS(RAD(90.833)) : REM Zenith distance
H1 = (H2 + H0) / 2 : REM Hour angle at half-hour
D1 = (D2 + D0) / 2 : REM declination at half-hour
R0 = S * SIN(D0) + C * COS(D0) * COS(H0) - Z
R1 = S * SIN(D1) + C * COS(D1) * COS(H1) - Z
R2 = S * SIN(D2) + C * COS(D2) * COS(H2) - Z
IF SGN(R0) = SGN(R2) THEN ENDPROC : REM No sunrise or sunset this hour
A = 2 * R2 - 4 * R1 + 2 * R0
B = 4 * R1 - 3 * R0 - R2
D = B * B - 4 * A * R0 : IF D < 0 THEN ENDPROC : REM No solution to quadratic
D = SQR(D)
E = (0 - B + D) / (2 * A)
IF E > 1 OR E < 0 THEN E = (0 - B - D) / (2 * A)
U = H0 + E * (H2 - H0)
V = 0 - COS(D1) * SIN(U)
W = C * SIN(D1) - S * COS(D1) * COS(U)
Z = DEG(ATN(V / W))
IF W < 0 THEN Z += 180
IF Z < 0 THEN Z += 360
IF Z > 360 THEN Z -= 360
IF R0 < 0 THEN tr = N + E : ar = Z ELSE ts = N + E : as = Z
ENDPROC
REM Sun's position in the sky (Van Flandern & Pulkkinen, 1979)
REM D = Delta days from noon on 1st January 2000
DEF PROCsun(D, RETURN RA, RETURN Dec)
LOCAL L, G, T, U, V, W
T = D / 36525 + 1 : REM Number of centuries since 1900
L = 0.779072 + 0.00273790931 * D
G = 0.993126 + 0.00273777850 * D
L = (L - INT(L)) * 2 * PI : REM Mean latitude (radians)
G = (G - INT(G)) * 2 * PI : REM Mean anomaly (radians)
U = 1 - 0.03349 * COS(G) - 0.00014 * COS(2 * L) + 0.00008 * COS(L)
V = (0.39785 - 0.00021 * T) * SIN(L) - 0.01 * SIN(L - G) + 0.00333 * SIN(L + G)
W = (0.03211 - 0.00008 * T) * SIN(G) - 0.04129 * SIN(2 * L) + \
\ 0.00104 * SIN(2 * L - G) - 0.00035 * SIN(2 * L + G) - 0.0001
REM Compute Sun's Right Ascension and Declination (radians)
RA = L + ASN(W / SQR(U - V * V))
Dec = ASN(V / SQR(U))
ENDPROC