Function range(x)
tpi = 2 * 3.1415927
b = x / tpi
a = tpi * (b - Int(b))
If (a < 0) Then a = tpi + a
range = a
End Function
Function jday(y, m, d)
J = 367 * y - Int(7 * (y + Int((m + 9) / 12)) / 4)
J = J - Int(3 * (Int((y + Int((m - 9) / 7)) / 100) + 1) / 4)
J = J + Int(275 * m / 9) + d + 1721029
jday = J
End Function
Function Lambda(jday)
rads = 3.14159 / 180
L = range(280.461 * rads + 0.9856474 * rads * jday)
g = range(357.528 * rads + 0.9856003 * rads * jday)
Lambda = range(L + 1.915 * rads * Sin(g) + 0.02 * rads * Sin(2 * g))
End Function
Function oblique(jday)
rads = 3.14159 / 180
oblique = 23.439 * rads - 0.0000004 * rads * jday
End Function
Function delta(jd)
sdelta = Sin(oblique(jd)) * Sin(Lambda(jd))
delta = Atn(sdelta / Sqr(1 - sdelta * sdelta))
End Function
Function alpha(jd)
Pi = 3.14159
x = Cos(Lambda(jd))
y = Cos(oblique(jd)) * Sin(Lambda(jd))
alpha = Atn(y / x)
If ((x < 0) And (y < 0)) Then alpha = alpha - Pi
If ((x < 0) And (y > 0)) Then alpha = alpha + Pi
End Function
Function timeq(jd)
rads = 3.14159 / 180
L = range(280.461 * rads + 0.9856474 * rads * jd)
korj = 720 * (L - alpha(jd)) / 3.14159
timeq = 1440 - korj
End Function
Function hrmnss(dectime)
hours = Int(dectime)
minutes = 60 * (dectime - hours)
seconds = Int(60 * (minutes - Int(minutes)))
minutes = Int(minutes)
hrmnss = TimeSerial(hours, minutes, seconds)
End Function
Function RiseSet(f0, longit, timeq)
Pi = 3.1415927
RiseSet = 12 + 12 * f0 / Pi - longit / 15 + timeq / 60
End Function
Function f0(declin, latit)
Pi = 3.1415927
refract = Pi * 17 / 30 / 180
Sundia = Pi * 4 / 15 / 180
f0 = 0
x1 = Tan(declin + 0.5 * Sundia + refract) * Tan(Pi * latit / 180)
If (x1 >= 1) Then f0 = Pi
If (x1 > -1) And (x1 < 1) Then f0 = Pi / 2 + Atn(x1 / Sqr(1 - x1 * x1))
End Function
Function dif_jd(jour, mois, année)
dif_jd = DateDiff("d", "01/01/2000", DateSerial(année, mois, jour))
End Function
Function jd_2k()
jd_2k = jday(2000, 1, 1)
End Function
Function lever(latitude, longitude, jour, mois, année)
equation_t = timeq(dif_jd(jour, mois, année))
lever = hrmnss(RiseSet(-f0(delta(dif_jd(jour, mois, année)), latitude), longitude, equation_t))
End Function
Function coucher(latitude, longitude, jour, mois, année)
equation_t = timeq(dif_jd(jour, mois, année))
coucher = hrmnss(RiseSet(f0(delta(dif_jd(jour, mois, année)), latitude), longitude, equation_t))
End Function