Now compute the orbital elements of the planet of interest. If you
want the position of the Sun or the Moon, you only need to compute
the orbital elements for the Sun or the Moon. If you want the
position of any other planet, you must compute the orbital elements
for that planet and for the Sun (of course the orbital
elements for the Sun are really the orbital elements for the Earth;
however it's customary to here pretend that the Sun orbits the
Earth). This is necessary to be able to compute the geocentric
position of the planet.

Please note that a, the semi-major axis, is given in Earth radii for
the Moon, but in Astronomical Units for the Sun and all the planets.

When computing M (and, for the Moon, when computing N and w as well),
one will quite often get a result that is larger than 360 degrees, or
negative (all angles are here computed in degrees). If negative,
add 360 degrees until positive. If larger than 360 degrees, subtract
360 degrees until the value is less than 360 degrees. Note that, in
most programming languages, one must then multiply these angles with
pi/180 to convert them to radians, before taking the sine or cosine
of them.

atan2() is a function that converts an x,y coordinate pair to the
correct angle in all four quadrants. It is available as a library
function in Fortran, C and C++. In other languages, one has to write
one's own atan2() function. It's not that difficult:

(since the Sun always is in the ecliptic plane, zs is of course zero).
xs,ys is the Sun's position in a coordinate system in the plane of
the ecliptic. To convert this to equatorial, rectangular, geocentric
coordinates, compute:

xe = xs
ye = ys * cos(ecl)
ze = ys * sin(ecl)

Finally, compute the Sun's Right Ascension (RA) and Declination (Dec):

HA is usually given in the interval -12 to +12 hours, or -180 to +180
degrees. If HA is zero, the object can be seen directly to the south.
If HA is negative, the object is to the east of south, and if HA is
positive, the object is to the west of south. IF your computed HA
should fall outside this interval, add or subtract 24 hours (or 360
degrees) until HA falls within this interval.

Now it's time to convert our objects HA and Decl to local azimuth and
altitude. To do that, we also must know lat, our local latitude. Then
we proceed as follows:

This completes our calculation of the local azimuth and altitude.
Note that azimuth is 0 at North, 90 deg at East, 180 deg at South
and 270 deg at West. Altitude is of course 0 at the (mathematical)
horizon, 90 deg at zenith, and negative below the horizon.

(Note that if decl is exactly 90 deg, cos(Decl) becomes zero and we
get a division by zero when computing topRA, but that formula breaks
down only very close to the celestial poles anyway and we never see
the Moon there. Also if gclat is precisely zero, g becomes zero too,
and we get a division by zero when computing topDecl. In that case,
replace the formula for topDecl with

topDecl = Decl - mpar * rho * sin(-Decl) * cos(HA)

which is valid for gclat equal to zero; it can also be used for gclat
extremely close to zero).

This correction to topocentric position can also be applied to the
Sun and the planets. But since they're much farther away, the
correction becomes much smaller. It's largest for Venus at inferior
conjunction, when Venus' parallax is somewhat larger than 32 arc
seconds. Within our aim of obtaining a final accuracy of 1-2 arc
minutes, it might barely be justified to correct to topocentric
position when Venus is close to inferior conjunction, and perhaps
also when Mars is at a favourable opposition. But in all other cases
this correction can safely be ignored within our accuracy aim. We
only need to worry about the Moon in this case.

If you want to compute topocentric coordinates for the planets too,
you do it the same way as for the Moon, with one exception: the
Moon's parallax is replaced by the parallax of the planet (ppar),
as computed from this formula:

ppar = (8.794/3600)_deg / r

where r is the distance of the planet from the Earth, in astronomical
units.

cbrt() is the cube root function: cbrt(x) = x**(1.0/3.0). The
formulae has been devised so that both g+h and g-h always are
positive. Therefore one can here safely compute cbrt(x) as
exp(log(x)/3.0) . In general, cbrt(-x) = -cbrt(x) and of course
cbrt(0) = 0.

Instead of trying to compute some eccentric anomaly, we compute the
true anomaly and the heliocentric distance directly:

v = 2.0 * atan(s)
r = q * ( 1.0 + s*s )

When we know the true anomaly and the heliocentric distance, we can
proceed by computing the position in space (section 7).