'Astro-functions
'Keith@xylem.demon.co.uk
'http://www.xylem.demon.co.uk/kepler/
'This text file contains the current set of User defined functions
'for Microsoft Excel 5, 95 and 97 versions. The functions are written
'in Visual Basic for Applications, and may be useable from other
'Office97 applications.
'RA and time UT and all angles measured in degrees, all instants
'defined by the decimal fraction of days after J2000.0
'Value of daily change in Pluto's Longitude corrected 19990227
Public Const pi As Double = 3.14159265358979
Public Const tpi As Double = 6.28318530717958
Public Const degs As Double = 57.2957795130823
Public Const rads As Double = 1.74532925199433E-02
' Trig functions working in degrees
Function DegSin(x As Double) As Double
DegSin = Sin(rads * x)
End Function
Function DegCos(x As Double) As Double
DegCos = Cos(rads * x)
End Function
Function DegTan(x As Double) As Double
DegTan = Tan(rads * x)
End Function
Function DegArcsin(x As Double) As Double
DegArcsin = degs * Application.Asin(x)
End Function
Function DegArccos(x As Double) As Double
DegArccos = degs * Application.Acos(x)
End Function
Function DegArctan(x As Double) As Double
DegArctan = degs * Atn(x)
End Function
' modified arctan function - returns tangent in right quadrant
Function DegAtan2(y As Double, x As Double) As Double
DegAtan2 = degs * Application.Atan2(x, y)
If DegAtan2 < 0 Then DegAtan2 = DegAtan2 + 360
End Function
' functions below return argument in range 0 to 2 pi or 0 to 360
Private Function range2pi(x)
range2pi = x - tpi * Int(x / tpi)
End Function
Private Function range360(x)
range360 = x - 360 * Int(x / 360)
End Function
' converts degs min sec to decimal deg
Function degdecimal(d, m, s)
degdecimal = d + m / 60 + s / 3600
End Function
' returns Julian day number, greg is switch for Gregorian
' or Julian calendars. Default is gregorian.
Function jday(year As Integer, month As Integer, day As Integer, _
hour As Integer, min As Integer, sec As Double, Optional greg) _
As Double
Dim a As Double
Dim b As Integer
a = 10000# * year + 100# * month + day
If (a < -47120101) Then
MsgBox "Warning: date too early for algorithm"
End If
If (IsMissing(greg)) Then greg = 1
If (month <= 2) Then
month = month + 12
year = year - 1
End If
If (greg = 0) Then
b = -2 + Fix((year + 4716) / 4) - 1179
Else
b = Fix(year / 400) - Fix(year / 100) + Fix(year / 4)
End If
a = 365# * year + 1720996.5
jday = a + b + Fix(30.6001 * (month + 1))
jday = jday + day + (hour + min / 60 + sec / 3600) / 24
End Function
' Returns days after J2000, greg switch for Gregorian or
' Julian calendar - default is Gregorian.
Function day2000(year As Integer, month As Integer, day As Integer, _
hour As Integer, min As Integer, sec As Double, Optional greg) _
As Double
Dim a As Double
Dim b As Integer
If (IsMissing(greg)) Then greg = 1
a = 10000# * year + 100# * month + day
If (month <= 2) Then
month = month + 12
year = year - 1
End If
If (greg = 0) Then
b = -2 + Fix((year + 4716) / 4) - 1179
Else
b = Fix(year / 400) - Fix(year / 100) + Fix(year / 4)
End If
a = 365# * year - 730548.5
day2000 = a + b + Fix(30.6001 * (month + 1))
day2000 = day2000 + day + (hour + min / 60 + sec / 3600) / 24
End Function
' Given polar coords, returns cartesian coords
' r is radius vector, theta is declination like angle, phi is
' RA like angle. Index = 1 returns x coord and so on
Function rectangular(r As Double, theta As Double, phi As Double, _
index As Integer) As Double
Dim r_cos_theta As Double
r_cos_theta = r * DegCos(theta)
Select Case index
Case 1
rectangular = r_cos_theta * DegCos(phi) 'returns x coord
Case 2
rectangular = r_cos_theta * DegSin(phi) 'returns y coord
Case 3
rectangular = r * DegSin(theta) 'returns z coord
End Select
End Function
' given cartesian coords, returns distance
Function rlength(x As Double, y As Double, z As Double) As Double
rlength = Sqr(x * x + y * y + z * z)
End Function
' returns spherical coord given cartesian coords
Function spherical(x As Double, y As Double, z As Double, _
index As Integer) As Double
Dim rho As Double
rho = x * x + y * y
Select Case index
Case 1
spherical = Sqr(rho + z * z) 'returns r
Case 2
rho = Sqr(rho)
spherical = DegArctan(z / rho) 'returns theta
Case 3
rho = Sqr(rho)
spherical = DegAtan2(y, x) 'returns phi
End Select
End Function
' returns obliquity of equator given date in days after J2000.0
Function obliquity(d As Double) As Double
Dim t As Double
t = d / 36525 'julian centuries since J2000.0
obliquity = -(46.815 + (0.00059 - 0.001813 * t) * t) * t / 3600#
obliquity = obliquity + 23.43929111
End Function
' given geocentric ecliptic coordinates, returns geocentric
' equatorial coordinates - all cartesian. NB, d is epoch of
' frame, use d = 0 if working in J2000.0 frame
Function requatorial(x As Double, y As Double, z As Double, _
d As Double, index As Integer) As Double
Dim obl As Double
obl = obliquity(d)
Select Case index
Case 1
requatorial = x
Case 2
requatorial = y * DegCos(obl) - z * DegSin(obl)
Case 3
requatorial = y * DegSin(obl) + z * DegCos(obl)
End Select
End Function
' given geocentric equatorial coords, returns geocentric
' ecliptic coords - all cartesian. NB, d is epoch of
' frame, use d = 0 if working in J2000.0 frame
Function recliptic(x As Double, y As Double, z As Double, _
d As Double, index As Integer) As Double
Dim obl As Double
obl = obliquity(d)
Select Case index
Case 1
recliptic = x
Case 2
recliptic = y * DegCos(obl) + z * DegSin(obl)
Case 3
recliptic = -y * DegSin(obl) + z * DegCos(obl)
End Select
End Function
' returns equatorial coords given geocentric ecliptic coords
Function sequatorial(r As Double, theta As Double, phi As Double, _
d As Double, index As Integer) As Double
Dim x As Double, y As Double, z As Double
x = rectangular(r, theta, phi, 1)
y = rectangular(r, theta, phi, 2)
z = rectangular(r, theta, phi, 3)
sequatorial = spherical(requatorial(x, y, z, d, 1), _
requatorial(x, y, z, d, 2), _
requatorial(x, y, z, d, 3), index)
End Function
' returns ecliptic coords given geocentric equatorial coords
Function secliptic(r As Double, theta As Double, phi As Double, _
d As Double, index As Integer) As Double
Dim x As Double, y As Double, z As Double
x = rectangular(r, theta, phi, 1)
y = rectangular(r, theta, phi, 2)
z = rectangular(r, theta, phi, 3)
secliptic = spherical(recliptic(x, y, z, d, 1), _
recliptic(x, y, z, d, 2), _
recliptic(x, y, z, d, 3), index)
End Function
' approximate presession from one epoch d2 to another d1
Function precess(d1 As Double, d2 As Double, dec As Double, _
ra As Double, index As Integer, _
Optional ddec, Optional dra) As Double
Dim m As Double, n As Double, t As Double
If (IsMissing(dra)) Then dra = 0
If (IsMissing(ddec)) Then ddec = 0
t = d1 / 36525 'years since J2000
m = 0.01281233333333 + 0.00000775 * t
n = 0.005567527777778 - 2.361111111111E-06 * t
t = (d2 - d1) / 365.25
Select Case index
Case 1 'dec
precess = dec + (n * DegCos(ra) + ddec / 3600) * t
Case 2 'ra
precess = ra + (m + n * DegSin(ra) * DegTan(dec) + dra / 3600) * t
End Select
End Function
' returns ecliptic longitude and distance of Sun given
' days after J2000.0. Based on method in AA
Function ssun(d As Double, index As Integer) As Double
Dim G As Double
Dim L As Double
Dim lambda As Double
G = range360(357.528 + 0.9856003 * d)
L = range360(280.461 + 0.9856474 * d)
Select Case index
Case 1
ssun = 1.00014 - 0.01671 * DegCos(G) - 0.00014 * DegCos(2 * G)
Case 2
ssun = 0 'ecliptic latitude of Sun is zero
Case 3
ssun = range360(L + 1.915 * DegSin(G) + 0.02 * DegSin(2 * G))
Case 4 'equation of time
lambda = range360(L + 1.915 * DegSin(G) + 0.02 * DegSin(2 * G))
ssun = -1.915 * DegSin(G) - 0.02 * DegSin(2 * G)
ssun = ssun + 2.466 * DegSin(2 * lambda)
ssun = ssun - 0.053 * DegSin(4 * lambda)
End Select
End Function
' caretsian version of above
Function rsun(d As Double, index As Integer) As Double
Dim x As Double
Dim y As Double
Dim z As Double
rsun = rectangular(ssun(d, 1), ssun(d, 2), ssun(d, 3), index)
End Function
Function sun(d As Double, index As Integer) As Double
sun = sequatorial(ssun(d, 1), ssun(d, 2), ssun(d, 3), d, index)
End Function
' returns ecliptic longitide, latitude and distance of Moon
' (geocentric) given days after J2000.0
Function smoon(ByVal d As Double, index As Integer) As Double
Dim Nm As Double, im As Double, wm As Double, _
am As Double, ecm As Double, _
Mm As Double, em As Double, Ms As Double, _
ws As Double, xv As Double, _
yv As Double, vm As Double, rm As Double, x As Double, _
y As Double, z As Double, lon As Double, lat As Double, _
ls As Double, lm As Double, _
dm As Double, F As Double, dlong As Double, dlat As Double
d = d + 1.5 'epoch of theory is not same as J2000.0
Nm = range360(125.1228 - 0.0529538083 * d) * rads
im = 5.1454 * rads
wm = range360(318.0634 + 0.1643573223 * d) * rads
am = 60.2666 '(Earth radii)
ecm = 0.0549
Mm = range360(115.3654 + 13.0649929509 * d) * rads
em = Mm + ecm * Sin(Mm) * (1# + ecm * Cos(Mm))
xv = am * (Cos(em) - ecm)
yv = am * (Sqr(1# - ecm * ecm) * Sin(em))
vm = Application.Atan2(xv, yv)
rm = Sqr(xv * xv + yv * yv)
x = Cos(Nm) * Cos(vm + wm) - Sin(Nm) * Sin(vm + wm) * Cos(im)
x = rm * x
y = Sin(Nm) * Cos(vm + wm) + Cos(Nm) * Sin(vm + wm) * Cos(im)
y = rm * y
z = rm * (Sin(vm + wm) * Sin(im))
lon = Application.Atan2(x, y)
If lon < 0 Then lon = tpi + lon
lat = Atn(z / Sqr(x * x + y * y))
ws = range360(282.9404 + 0.0000470935 * d) * rads
Ms = range360(356.047 + 0.9856002585 * d) * rads
ls = Ms + ws 'Mean Longitude of the Sun (Ns=0)
lm = Mm + wm + Nm 'Mean longitude of the Moon
dm = lm - ls 'Mean elongation of the Moon
F = lm - Nm 'Argument of latitude for the Moon
Select Case index
Case 1 ' distance terms earth radii
rm = rm - 0.58 * Cos(Mm - 2 * dm)
rm = rm - 0.46 * Cos(2 * dm)
smoon = rm
Case 2 ' latitude terms
dlat = -0.173 * Sin(F - 2 * dm)
dlat = dlat - 0.055 * Sin(Mm - F - 2 * dm)
dlat = dlat - 0.046 * Sin(Mm + F - 2 * dm)
dlat = dlat + 0.033 * Sin(F + 2 * dm)
dlat = dlat + 0.017 * Sin(2 * Mm + F)
smoon = lat * degs + dlat
Case 3 ' longitude terms
dlon = -1.274 * Sin(Mm - 2 * dm) '(the Evection)
dlon = dlon + 0.658 * Sin(2 * dm) '(the Variation)
dlon = dlon - 0.186 * Sin(Ms) '(the Yearly Equation)
dlon = dlon - 0.059 * Sin(2 * Mm - 2 * dm)
dlon = dlon - 0.057 * Sin(Mm - 2 * dm + Ms)
dlon = dlon + 0.053 * Sin(Mm + 2 * dm)
dlon = dlon + 0.046 * Sin(2 * dm - Ms)
dlon = dlon + 0.041 * Sin(Mm - Ms)
dlon = dlon - 0.035 * Sin(dm) '(the Parallactic Equation)
dlon = dlon - 0.031 * Sin(Mm + Ms)
dlon = dlon - 0.015 * Sin(2 * F - 2 * dm)
dlon = dlon + 0.011 * Sin(Mm - 4 * dm)
smoon = lon * degs + dlon
End Select
End Function
' rectangular version of above
Function rmoon(d As Double, index As Integer) As Double
rmoon = rectangular(smoon(d, 1), smoon(d, 2), smoon(d, 3), index)
End Function
' returns spherical equatorial coords of moon given days
' after J2000.0
Function moon(d As Double, index As Integer) As Double
Dim nigel As Double
If index = 4 Then
nigel = smoon(d, 2)
moon = d
Else
moon = sequatorial(smoon(d, 1), smoon(d, 2), smoon(d, 3), d, index)
End If
End Function
' returns true anomaly given eccentricity of orbit, mean anomaly
' and optional precision parameter. Default precision is 10^-8 rads
Function kepler(m As Double, ecc As Double, Optional eps) As Double
Dim delta As Double, E As Double, v As Double
E = m 'first guess
delta = 0.05 'set delta equal to a dummy value
If (IsMissing(eps)) Then eps = 8 'if no eps then assume 10^-8
Do While Abs(delta) >= 10 ^ -eps 'converged?
delta = E - ecc * Sin(E) - m 'new error
E = E - delta / (1 - ecc * Cos(E)) 'corrected guess
Loop
v = 2 * Atn(((1 + ecc) / (1 - ecc)) ^ 0.5 * Tan(0.5 * E))
If v < 0 Then v = v + tpi
kepler = v
End Function
' returns heliocentric rectangular coords of planet
' given days after J2000.o and planet number
Function rplanet(d As Double, pnumber As Integer, _
index As Integer) As Double
Dim x As Double, y As Double, z As Double, v As Double, _
m As Double, i As Double, o As Double, p As Double, _
a As Double, E As Double, _
L As Double, r As Double
element i, o, p, a, E, L, d, pnumber 'calls the sub element
m = range2pi(L - p)
v = kepler(m, E, 8)
r = a * (1 - E * E) / (1 + E * Cos(v))
Select Case index
Case 1 'x coordinate
rplanet = Cos(o) * Cos(v + p - o) - Sin(o) * Sin(v + p - o) * Cos(i)
rplanet = rplanet * r
Case 2 'y coordinate
rplanet = Sin(o) * Cos(v + p - o) + Cos(o) * Sin(v + p - o) * Cos(i)
rplanet = rplanet * r
Case 3 'z coordinate
rplanet = r * (Sin(v + p - o) * Sin(i))
End Select
End Function
' gives equatorial spherical coords of planet given
' days after J2000.0 and planet number. Refeered to
' equinox of date
Function planet(d As Double, pnumber As Integer, _
index As Integer) As Double
Dim x As Double, y As Double, z As Double, v As Double, _
m As Double, i As Double, o As Double, p As Double, _
a As Double, E As Double, L As Double, r As Double, _
xe As Double, ye As Double, ze As Double, s1 As Double, _
si As Double, so As Double, c1 As Double, ci As Double, _
co As Double
element i, o, p, a, E, L, d, pnumber
m = range2pi(L - p)
v = kepler(m, E, 8)
r = a * (1 - E * E) / (1 + E * Cos(v))
s1 = Sin(v + p - o)
si = Sin(i)
so = Sin(o)
c1 = Cos(v + p - o)
ci = Cos(i)
co = Cos(o)
x = r * (co * c1 - so * s1 * ci)
y = r * (so * c1 + co * s1 * ci)
z = r * (s1 * si)
element i, o, p, a, E, L, d, 3
m = range2pi(L - p)
v = kepler(m, E, 8)
r = a * (1 - E * E) / (1 + E * Cos(v))
s1 = Sin(v + p - o)
si = Sin(i)
so = Sin(o)
c1 = Cos(v + p - o)
ci = Cos(i)
co = Cos(o)
xe = r * (co * c1 - so * s1 * ci)
ye = r * (so * c1 + co * s1 * ci)
ze = r * (s1 * si)
x = x - xe
y = y - ye
ecl = 23.429292 * rads 'value for J2000.0 frame
xe = x
ye = y * Cos(ecl) - z * Sin(ecl)
ze = y * Sin(ecl) + z * Cos(ecl)
Select Case index
Case 3
planet = Application.Atan2(xe, ye) * degs
If planet < 0 Then planet = 360 + planet
Case 2
planet = Atn(ze / Sqr(xe * xe + ye * ye)) * degs
Case 1
planet = Sqr(xe * xe + ye * ye + ze * ze)
End Select
End Function
' not accesible as user defined function
' returns elements of orbits of planets given days after
' J2000.0 and planet number
Sub element(i As Double, o As Double, p As Double, _
a As Double, E As Double, L As Double, _
ByVal d As Double, ByVal pnum As Integer)
Select Case pnum
Case 1 'mercury
i = (7.00487 - 0.000000178797 * d) * rads
o = (48.33167 - 0.0000033942 * d) * rads
p = (77.45645 + 0.00000436208 * d) * rads
a = 0.38709893 + 1.80698E-11 * d
E = 0.20563069 + 0.000000000691855 * d
L = range2pi(rads * (252.25084 + 4.092338796 * d))
Case 2 'venus
i = (3.39471 - 0.0000000217507 * d) * rads
o = (76.68069 - 0.0000075815 * d) * rads
p = (131.53298 - 0.000000827439 * d) * rads
a = 0.72333199 + 2.51882E-11 * d
E = 0.00677323 - 0.00000000135195 * d
L = range2pi(rads * (181.97973 + 1.602130474 * d))
Case 3 'earth
i = (0.00005 - 0.000000356985 * d) * rads
o = (-11.26064 - 0.00013863 * d) * rads
p = (102.94719 + 0.00000911309 * d) * rads
a = 1.00000011 - 1.36893E-12 * d
E = 0.01671022 - 0.00000000104148 * d
L = range2pi(rads * (100.46435 + 0.985609101 * d))
Case 4 'mars
i = (1.85061 - 0.000000193703 * d) * rads
o = (49.57854 - 0.0000077587 * d) * rads
p = (336.04084 + 0.00001187 * d) * rads
a = 1.52366231 - 0.000000001977 * d
E = 0.09341233 - 0.00000000325859 * d
L = range2pi(rads * (355.45332 + 0.524033035 * d))
Case 5 'jupiter
i = (1.3053 - 0.0000000315613 * d) * rads
o = (100.55615 + 0.00000925675 * d) * rads
p = (14.75385 + 0.00000638779 * d) * rads
a = 5.20336301 + 0.0000000166289 * d
E = 0.04839266 - 0.00000000352635 * d
L = range2pi(rads * (34.40438 + 0.083086762 * d))
Case 6 'saturn
i = (2.48446 + 0.0000000464674 * d) * rads
o = (113.71504 - 0.0000121 * d) * rads
p = (92.43194 - 0.0000148216 * d) * rads
a = 9.53707032 - 0.0000000825544 * d
E = 0.0541506 - 0.0000000100649 * d
L = range2pi(rads * (49.94432 + 0.033470629 * d))
Case 7 'uranus
i = (0.76986 - 0.0000000158947 * d) * rads
o = (74.22988 + 0.0000127873 * d) * rads
p = (170.96424 + 0.0000099822 * d) * rads
a = 19.19126393 + 0.0000000416222 * d
E = 0.04716771 - 0.00000000524298 * d
L = range2pi(rads * (313.23218 + 0.011731294 * d))
Case 8 'neptune
i = (1.76917 - 0.0000000276827 * d) * rads
o = (131.72169 - 0.0000011503 * d) * rads
p = (44.97135 - 0.00000642201 * d) * rads
a = 30.06896348 - 0.0000000342768 * d
E = 0.00858587 + 0.000000000688296 * d
L = range2pi(rads * (304.88003 + 0.0059810572 * d))
Case 9 'pluto
i = (17.14175 + 0.0000000841889 * d) * rads
o = (110.30347 - 0.0000002839 * d) * rads
p = (224.06676 - 0.00000100578 * d) * rads
a = 39.48168677 - 0.0000000210574 * d
E = 0.24880766 + 0.00000000177002 * d
L = range2pi(rads * (238.92881 + 3.97557152635181E-03 * d))
End Select
End Sub
' returns time of sunrise/set or twilights as angles
' given days after J2000.0, longitude and latitude of observer
' and an optional altitude required for the Sun. default altitude
' corresponds to upper limb on mathematical horizon
Function Sunrise(ByVal day As Double, _
glat As Double, glong As Double, index As Integer, _
Optional altitude) As Double
Dim sinalt As Double, gha As Double, lambda As Double, _
delta As Double, t As Double, c As Double, days As Double, _
utold As Double, utnew As Double, sinphi As Double, _
cosphi As Double, L As Double, G As Double, E As Double, _
obl As Double, signt As Double, act As Double
If (IsMissing(altitude)) Then altitude = -0.833
utold = 180
utnew = 0
sinalt = Sin(altitude * rads)
sinphi = DegSin(glat)
cosphi = DegCos(glat)
If index = 1 Then signt = 1 Else signt = -1
Do While Abs(utold - utnew) > 0.1
utold = utnew 'carry forward the iteration
days = day + utold / 360
t = days / 36525
L = range360(280.46 + 36000.77 * t)
G = 357.528 + 35999.05 * t
lambda = L + 1.915 * DegSin(G) + 0.02 * DegSin(2 * G)
E = -1.915 * DegSin(G) - 0.02 * DegSin(2 * G) _
+ 2.466 * DegSin(2 * lambda) - 0.053 * DegSin(4 * lambda)
obl = 23.4393 - 0.13 * t
gha = utold - 180 + E
delta = DegArcsin(DegSin(obl) * DegSin(lambda))
c = (sinalt - sinphi * DegSin(delta)) / (cosphi * DegCos(delta))
act = DegArccos(c)
If c > 1 Then act = 0
If c < -1 Then act = 180
utnew = range360(utold - (gha + glong + signt * act))
Loop
Sunrise = utnew
End Function
' returns topocentric coords of moon given days after J2000.0
' observer's latitude and longitude.
Function tmoon(days As Double, glat As Double, _
glong As Double, index As Integer) As Double
Dim x As Double, y As Double, z As Double, xt As Double, _
yt As Double, zt As Double, r As Double, xe As Double, _
ye As Double, ze As Double, lst As Double
x = rmoon(days, 1)
y = rmoon(days, 2)
z = rmoon(days, 3)
xe = requatorial(x, y, z, days, 1)
ye = requatorial(x, y, z, days, 2)
ze = requatorial(x, y, z, days, 3)
r = rlength(x, y, z)
lst = gst(days) + glong
xt = xe - DegCos(glat) * DegCos(lst)
yt = ye - DegCos(glat) * DegSin(lst)
zt = ze - DegSin(glat)
tmoon = spherical(xt, yt, zt, index)
End Function
' returns siderial time at longitude zero given days
' after J2000.0
Function gst(days As Double) As Double
Dim t As Double
t = days / 36525
gst = 280.46061837 + 360.98564736629 * days
gst = gst + 0.000387933 * t ^ 2 - t ^ 3 / 38710000
gst = range360(gst)
End Function
' returns time of rise or set of object given RA and DEC
' of object, latitude and longitude of observer and days
' before J2000.0
Function riset(d As Double, dec As Double, ra As Double, _
glat As Double, glong As Double, index As Integer) As Double
Dim lst As Double, tau As Double
lst = gst(d) + glong
tau = DegSin(-34 / 60) - DegSin(glat) * DegSin(dec)
tau = tau / (DegCos(glat) * DegCos(dec))
tau = DegArccos(tau)
Select Case index
Case 1 'rising, degrees on day after midnight chosen
riset = range360(360 - 0.9972696 * (lst - ra + tau))
Case 2 'transit on day after or after
If ra > lst Then
riset = Abs(0.9972696 * (lst - ra))
Else
riset = 360 - 0.9972696 * (lst - ra)
End If
Case 3 'setting, degrees UT after midnight chosen
riset = range360(0.9972696 * (ra + tau - lst))
End Select
End Function
' returns altitude and azimuth of object given RA and DEC of
' object, latitude and longitude of observer and days after
' J2000.0
Function altaz(d As Double, dec As Double, ra As Double, _
glat As Double, glong As Double, index As Integer) As Double
Dim h As Double, lst As Double, a As Double, z As Double, _
sa As Double, cz As Double
h = gst(d) + glong - ra
sa = DegSin(dec) * DegSin(glat)
sa = sa + DegCos(dec) * DegCos(glat) * DegCos(h)
a = DegArcsin(sa)
cz = DegSin(dec) - DegSin(glat) * sa
cz = cz / (DegCos(glat) * DegCos(a))
Select Case index
Case 1 'altitude
altaz = a
Case 2 'azimuth
If DegSin(h) < 0 Then
altaz = DegArccos(cz)
Else
altaz = 360 - DegArccos(cz)
End If
End Select
End Function
' calculates the time of rise and set of the planets given
' days after J2000.0, planet number and longitude and latitude
' of observer
Function prise(d As Double, p As Integer, glat As Double, _
glong As Double, index As Integer) As Double
Dim lst As Double, tau As Double, ra As Double, dec As Double, _
ut0 As Double, ut1 As Double, slt As Double, clt As Double, _
flag As Integer, count As Integer
count = 0
flag = 2 - index
ut0 = 0
ut1 = 180
slt = DegSin(glat)
clt = DegCos(glat)
Do While Abs(ut0 - ut1) > 0.1 And count < 10
count = count + 1
ut0 = ut1
dec = planet(d + ut0 / 360, p, 2)
ra = planet(d + ut0 / 360, p, 3)
tau = (-0.009890037858741 - slt * DegSin(dec)) / (clt * DegCos(dec))
Select Case tau
Case Is >= 1
tau = 180
Case Is <= -1
tau = -180
Case Else
tau = DegArccos(tau)
End Select
lst = gst(d + ut0 / 360) + glong
ut1 = range360(ut0 - 0.9972696 * (lst - ra + flag * tau))
Loop
prise = ut1
End Function