The Moon is about half a degree in diameter. The routine shown here will give you the apparent position of the moon good to within an eighth of the diameter, and will be within 1/15 of a diameter in 4 cases out of five. This accuracy is good enough for finding the Moonrise and set, and for calculations of the phase, paralactic angle, pole angle and approximate libration points.
The method shown here is taken from Paul Schlyter's page "How to compute planetary positions". The algorithm is based on an article by van Flandern and Pulkkinen, hence I refer to this method as 'vFPS' for short.
I compared the accuracy of this vFPS method, and the method found on page D76 of the Astronomical Almanac with the output from the Interactive Computer Ephemeris (ICE) for an 18.6 year period either side of J2000.0. The Moon approximately repeats the positions in its complex orbit over this so-called 'Saros' cycle, so a simple statistical approach to the errors may be feasible.
The QBASIC program will ask you for the year, month, day and hour
and minute of the instant you want the position for, and will then
calculate the geocentric RA and DEC of the Moon for the equinox of
and obliquity of date. Time should be in UT (strictly TDT, but the
difference of about 1 min is not significant at the accuracy level
used here).
The easiest way to extract the program is as follows;
**** below
DEF FNday function definition,
which should all be on one line including the last /24 statement
'moon2.bas' or some similar title (remember to click the
'all files *.*' option in the 'Save as type' box)
QBASIC
'************************************************************************************ ' Moon positions to a few minutes of arc ' From Paul Schlyter's page at ' http://hotel04.ausys.se/pausch/comp/ppcomp.html ' ' The QBASIC code here is a literal translation of Paul's ' formulas, you can probably devise a more efficient ' version! ' ' definitions and functions DEFDBL A-Z pr$ = "####.###" pr2$ = "#########.######" pi = 4 * ATN(1) tpi = 2 * pi twopi = tpi degs = 180 / pi rads = pi / 180 ' ' Get the days to Dec 31st 0h 2000 - note, this is NOT same ' as J2000 ' h is UT in decimal hours ' FNday only works between 1901 to 2099 - see Meeus chapter 7 ' DEF FNday (y, m, d, h) = 367 * y - 7 * (y + (m + 9) \ 12) \ 4 + 275 * m \ 9 + d - 730530 + h / 24 ' ' define atn2 ' ' the atn2 function below returns an angle in the range 0 to two pi ' depending on the signs of x and y. ' DEF FNatn2 (y, x) a = ATN(y / x) IF x < 0 THEN a = a + pi IF (y < 0) AND (x > 0) THEN a = a + tpi FNatn2 = a END DEF ' ' the function below returns the true integer part, ' even for negative numbers ' DEF FNipart (x) = SGN(x) * INT(ABS(x)) ' ' the function below returns an angle in the range ' 0 to two pi ' DEF FNrange (x) b = x / tpi a = tpi * (b - FNipart(b)) IF a < 0 THEN a = tpi + a FNrange = a END DEF ' CLS ' ' get the time (strictly TDT, but UT will do) and date ' INPUT " year : ", y INPUT " month : ", m INPUT " day : ", day INPUT "hour UT : ", h INPUT " minute : ", mins h = h + mins / 60 d = FNday(y, m, day, h) PRINT " days : "; d ' moon elements Nm = FNrange((125.1228 - .0529538083# * d) * rads) im = 5.1454 * rads wm = FNrange((318.0634 + .1643573223# * d) * rads) am = 60.2666 '(Earth radii) ecm = .0549 Mm = FNrange((115.3654 + 13.0649929509# * d) * rads) ' sun elements Ns = 0! isun = 0! ws = FNrange((282.9404 + 4.70935E-05 * d) * rads) asun = 1! '(AU) ecs = .016709 - 1.151E-09 * d Ms = FNrange((356.047 + .9856002585# * d) * rads) 'position of Sun 'Es = Ms + Es * SIN(Ms) * (1! + ecs * COS(Ms)) 'xv = COS(Es) - ecs 'yv = SQR(1! - ecs * ecs) * SIN(Es) 'vs = FNatn2(yv, xv) 'rs = SQR(xv * xv + yv * yv) 'lonsun = vs + ws 'xs = rs * COS(lonsun) 'ys = rs * SIN(lonsun) 'xe = xs 'ye = ys * COS(ecl) ''ze = ys * SIN(ecl) 'ras = FNatn2(ye, xe) 'decs = FNatn2(ze, sqr(xe * xe + ye * ye)) ' position of Moon Em = Mm + ecm * SIN(Mm) * (1! + ecm * COS(Mm)) xv = am * (COS(Em) - ecm) yv = am * (SQR(1! - ecm * ecm) * SIN(Em)) vm = FNatn2(yv, xv) rm = SQR(xv * xv + yv * yv) xh = rm * (COS(Nm) * COS(vm + wm) - SIN(Nm) * SIN(vm + wm) * COS(im)) yh = rm * (SIN(Nm) * COS(vm + wm) + COS(Nm) * SIN(vm + wm) * COS(im)) zh = rm * (SIN(vm + wm) * SIN(im)) ' moons geocentric long and lat lon = FNatn2(yh, xh) lat = FNatn2(zh, SQR(xh * xh + yh * yh)) ' perturbations ' first calculate arguments below, which should be in radians 'Ms, Mm Mean Anomaly of the Sun and the Moon 'Nm Longitude of the Moon's node 'ws, wm Argument of perihelion for the Sun and the Moon 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 ' then add the following terms to the longitude ' note amplitudes are in degrees, convert at end dlon = -1.274 * SIN(Mm - 2 * dm) '(the Evection) dlon = dlon + .658 * SIN(2 * dm) '(the Variation) dlon = dlon - .186 * SIN(Ms) '(the Yearly Equation) dlon = dlon - .059 * SIN(2 * Mm - 2 * dm) dlon = dlon - .057 * SIN(Mm - 2 * dm + Ms) dlon = dlon + .053 * SIN(Mm + 2 * dm) dlon = dlon + .046 * SIN(2 * dm - Ms) dlon = dlon + .041 * SIN(Mm - Ms) dlon = dlon - .035 * SIN(dm) '(the Parallactic Equation) dlon = dlon - .031 * SIN(Mm + Ms) dlon = dlon - .015 * SIN(2 * F - 2 * dm) dlon = dlon + .011 * SIN(Mm - 4 * dm) lon = dlon * rads + lon ' latitude terms dlat = -.173 * SIN(F - 2 * dm) dlat = dlat - .055 * SIN(Mm - F - 2 * dm) dlat = dlat - .046 * SIN(Mm + F - 2 * dm) dlat = dlat + .033 * SIN(F + 2 * dm) dlat = dlat + .017 * SIN(2 * Mm + F) lat = dlat * rads + lat ' distance terms earth radii rm = rm - .58 * COS(Mm - 2 * dm) rm = rm - .46 * COS(2 * dm) ' next find the cartesian coordinates ' of the geocentric lunar position xg = rm * COS(lon) * COS(lat) yg = rm * SIN(lon) * COS(lat) zg = rm * SIN(lat) ' rotate to equatorial coords ' obliquity of ecliptic of date ecl = (23.4393 - 3.563E-07 * d) * rads xe = xg ye = yg * COS(ecl) - zg * SIN(ecl) ze = yg * SIN(ecl) + zg * COS(ecl) ' geocentric RA and Dec ra = FNatn2(ye, xe) dec = ATN(ze / SQR(xe * xe + ye * ye)) PRINT USING pr$; ra * degs / 15 PRINT USING pr$; dec * degs END '************************************************************************************
For 1998, August 10th at 0h UT, the routine gives the following output;
lat : -1.91667
long : 52.5
year : 1998
month : 8
day : 10
hour UT : 0
minute : 0
days : -508
22.947 (RA geocentric)
-7.798 (Dec geocentric)
These results compare with Dr Ahmed Monsur's Mooncalc 4.0 working in geocentric mode as follows;
Moon Dec: -7.816 Moon RA: 22.948
As you can see, the results are within 1 arcmin in this case. I have investigated the accuracy more thoroughly below.
I wrote very simple QBASIC programs implementing two different 'low precision' algorithms for finding the Moon's position. Each of these programs was modified to generate daily positions for about 18.6 years (roughly one 'Saros') either side of the year 2000.
I then compared each of the two sets positions with the positions generated by the Interactive Computer Ephemeris (ICE), and I have tried to analyse the main features of the resulting error 'signal' using a simple statistical approach.
The ICE produces apparent geocentric positions in Right Ascension and Declination for the moon, and so code was added to the QBASIC routines to transform the co-ordinates to RA and DEC. A reference which produces the geocentric ecliptic coordinates would be better, as the resulting error signals could be related directly to the longitude and latitude series used.
The table below shows the errors in seconds of time (RA) and arcseconds (Dec) between the Interactive Computer Ephemeris moon position, and the two approximate methods discussed here.
vFPS D76
RA Dec RA Dec
max 27 265 97 811
rms 7 66 22 224
RA errors are in seconds of time, and Dec errors are in
arcseconds.
'Max' refers to the largest (absolute) error seen in any of the 13871 daily positions calculated. The term 'rms' refers to the standard deviation of the errors - which is similar to the 'root mean square' amplitude of a signal. The large ratio between the peak of the error 'signal' and the 'rms' level indicates that peaks are extreme 'spikes' on an otherwise smoother curve. This impresion is borne out by a graph of part of the RA error.
The table below shows the fraction of the approximate positions which fell within a certain interval of the ICE positions;
interval vFPS D76 (arcmin) RA Dec RA Dec -1 < e <1 0.44 0.60 0.14 0.21 -2 < e <2 0.78 0.94 0.28 0.41 -4 < e <4 0.99 1.00 0.53 0.72 -8 < e <8 ---- ---- 0.86 0.97
We could characterise an algorithm's performance by specifying a desired accuracy (i.e. 4 arcmin for rise and set calculations) and seeing what percentage of the positions fell within the interval. If the fraction is higher than 0.85 or so, that algorithm might be judged as 'good enough' for the purpose. We must take positions over at least one major period to prevent any systematic bias - here I have taken positions over two 18.6 year periods.
Using the vFPS method, you can be 99% sure that your Moon position is within 4 arcminutes of the apparent geocentric position as shown by the ICE, and the position will be within 2 arcminutes 'most of the time'. You have a roughly fifty-fifty chance of your position being within 1 arcmin.
The D76 method will give positions which are within 8 arcminutes of the apparent geocentric position 'most of the time'. This error is eqivalent to roughly one quarter of the width of the moon. You have a fifty-fifty chance of the position being within 4 arcmin. I would characterise this algorithm as roughly 'half to a quarter' as accurate as the vFPS version.
As the error time series has a periodic structure (the low precision methods are truncated Fourier-like series), this statistical viewpoint may not be entirely valid, but it does give a repeatable framework for assessing the accuracy of position algorithms.
[Root]
Last Modified 10th August 1998
Keith Burnett