'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