Julian Day (valid from 1900/3/1 to 2100/2/28)

Julian day: 86400 s, Julian year: 365.25 d, Julian Century: 36525 d

double JulianDay (int date, int month, int year, double UT)
{

if (month<=2) {month=month+12; year=year-1;}
return (int)(365.25*year) + (int)(30.6001*(month+1)) - 15 + 1720996.5 + date + UT/24.0;

}

Calculation of Lunar ecliptic coordinates can be found in:

Montenbruck, Oliver / Pfleger, Thomas: Astronomie mit dem Personal Computer; mit CD-ROM,
Springer Berlin, 4. Aufl. 2004, ISBN 3-540-21204-3

Montenbruck, Oliver / Pfleger, Thomas: Astronomy on the Personal Computer. with CD-ROM;
Springer Berlin, 4th ed. 2004, ISBN 3-540-67221-4

Meeus, Jean: Astronomical Algorithms
Willmann-Bell; Hardcover (1st ed. 1991), ISBN: 0943396352
Willmann-Bell; Hardcover (2nd ed. 1998), ISBN: 0943396611

Meeus, Jean: Astronomical Formulae for Calculators;
Willmann-Bell; Softcover, ISBN: 0943396220

Convert ecliptic latitude B and longitude L to right ascension RA and declination delta

obliquity of ecliptic:
T = (JD-2451545.0)/36525.0;
eps = 23.0 + 26.0/60.0 + 21.448/3600.0 - (46.8150*T+ 0.00059*T*T- 0.001813*T*T*T)/3600;

X = cos(B)*cos(L)
Y = cos(eps)*cos(B)*sin(L) - sin(eps)*sin(B)
Z = sin(eps)*cos(B)*sin(L) + cos(eps)*sin(B)
R = sqrt(1-Z*Z);

delta = (180/PI)*arctan(Z/R); // in degrees
or
delta = (180/PI)*arcsin[sin(eps)*cos(B)*sin(L) + cos(eps)*sin(B)]; // in degrees

RA = (24/PI)*arctan(Y/(X+R)); // in hours
or
RA = (180/PI)*arctan[(sin(L)*cos(eps) - tan(B)*sin(eps))/cos(L)]; // in degrees
        if (arctan<0) add PI to arctan

Compute sidereal time at Greenwich (according to: Jean Meeus: Astronomical Algorithms)

T = (JD - 2451545.0 ) / 36525;

theta0 = 280.46061837 + 360.98564736629*(JD-2451545.0) + 0.000387933*T*T - T*T*T/38710000.0; // degrees

Local sidereal time:

theta = theta0 + longitude (eastern longitudes positive, western negative)

Hour angle: tau = theta - RA

Convert (tau,delta) to horizon coordinates of the observer (altitude h, azimuth az)

sin(h) = sin(beta )*sin(delta) + cos(beta)*cos(delta)*cos(tau)

tan(az) = - sin(tau) / (cos(beta)*tan(delta) - sin(beta)*cos(tau))

Compute the parallax p in altitude:

horParal = 8.794 / (moonDistance/149.59787E6); // horizontal parallax (arcseconds), Meeus S. 263

sin(p) = cos(h)*sin(horParal/3600); // parallax in altitude (degrees)