Approximate astronomical positions

Moon position to 4 arcmin

[Root]

Contents

Overview

[Top]

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.

QBASIC program

[Top]

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;

'************************************************************************************
'   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.

Accuracy and comparison with the D76 method

[Top]

Method

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.

Results

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.

Conclusions

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

keith@xylem.demon.co.uk