The mathematical  analysis of lunar eclipses

 

by Byron W Soulsby

Theodore Lunar Observatory, Australia

 

The QuickBasic programs given below have been modified to run under Windows 98 or Windows XP. They have now been provided for downloading to any PC for execution. The CTR Version 2 programs and data files are in compressed form and can be found here. They need to be expanded before use onto any available drive in your root directory into a folder named "CTR2" (without the quotes).  At the present time data files are complete for all lunar eclipses up to and including 2012.

 

A presentation of the mathematical expressions used in my earlier development QuickBasic programs  for Crater Timing Reduction (CTR) with the Macintosh II computer

Introduction

 

The second lunar eclipse in 2005, on October 17 will be my 50th lunar eclipse analysed. While most of these events have not been observed personally and as many crater timings and images were contributed by overseas observers, the occasion seemed appropriate for the publication of the mathematics used in some of my computer programs.

 

Sequence

 

The flow of my analysis is shown in the menu of the Microsoft QuickBasic  (b) program as follows:

 

VALID DRIVES ARE : 1 for ZIP CTR, 2 for HD IIci, or 3 as specified

DRIVE NUMBER: 2

SELECTED PATH IS HD IIci:Dev:MS QB:CTR:

 

17:01:06             MENU for HD IIci:Dev:MS QB:CTR:*.BAS Files     11-17-2004

1-      Lunar Eclipse Circs. Opt. Frame Coords with ILEEPSCPT

2-      Predict data for Future Events

3-      Crater Timing Predictions with Co-ordinates PLOT

4-      Altitude/Azimuth at Primary Contacts

5-      Observers' Timings to data files

6-      Crater Timing Reduction with Co-ordinates (CTRci.BAS)

7-      Output Menu

8-      Oblateness Menu

9-      Reduce Schmidt's Timings with CTRCP

10-  Find angle delta & Fi for video measures with ViaX4

11-  Find Topo Long with TERM

12-  Exit

 

SELECT ONE:

Option 6 is the most critical in the series of reductions as it includes the analysis of timings for primary contacts, mid crater timings, and all four contacts of the umbra with the edges of each crater, known as single crater timings. The latter is found  from an original series of mathematical expressions developed by David Herald.

 

Mathematics (CTR)

 

Check for correct time zone, change to ET by T0' = T0 + DT,

 find TF the time interval

Correct crater coordinates for heliocentric libration in Latitude (y) and Longitude (x) of the crater:

x'c  = cos y0  sin (xo - l'0)

y'c  = cos b0 sin y0 - sin b'0 cos y0 cos (x0 -  l'0)

 

 where heliocentric librations are l'0 = l0 + TFDl0 and b'0 = b0 + TFDb0,

                                                                                                (where Db0 not = 0)

Find position of crater in polar coordinates rc, Qc

            where Q"c = Qc + c'

            where c' = c + TF (Dc/13) , values of c and Dc are at opposition

from corrected polar coordinates (Q"c, rc) and find (xc, yc)

find Moon's centre  coordinates        x = x1 + TFDx

                                                y = y1 + TFDy

translate coordinates to umbra centre x = x  -  xc S'c , (xc is negative, see note)

                                                y = y + ycS'c, where S'c = S'c1 + TFDS'c

from (x, y) find polar coordinates of crater (rc, Y)

Find umbral oblateness            a = 0 .9983247 +  0.001676 cos 2Y, or:

use Jean Meeus expression a = 1  -  (cos2d0 sin2Y) / 298.3294

find, pc = pc1 + TF Dpc and L' = S'c(1 -  rc)0.5 tan(S'0 -  p0)

find theoretical umbral radius, rt = F2 + L' where F2 = a pc  -  S'0 + p0

find  % umbral enlargement,    E = ( rt  -  ro) / rt * 100

 

Terms are:

 L' is the umbral enlargement due to projection from the crater's position on the Moon's surface to the fundamental plane, F2 is the umbral semi-diameter at the plane, rc is the apparent distance of the crater from the Moon's centre, in units of the Moon's semi-diameter Sc.

 

Notes:

 xc is positive for Besselian elements to the East, as the Moon's longitude is negative to the East, hence sign of x is changed.

s'0 = s0  -  allowance for irradiation

S'c = 0.000022o  -  0.272453 pco

For change in C (Moon's axis PA) use DC / 13 as this is an approximation for Sun/Earth - Earth/Moon velocity ratio, using rate of hourly motion in RA of Moon/Sun to achieve heliocentric rate of change from geocentric value tabulated in the Astronomical Ephemerides.

Crater Timing Prediction

 

In option 3 of the menu predictions are computed for crater timings for each eclipse. A complete data base of lunar features has been prepared to give latitude, longitude and diameter (in miles). The mathematics used in the program are as follows:

 

Correct crater coordinates for libration x'c = cos y sin (x  -  l')

                                                            y'c = cos b sin y  -  sin b cos y cos(x  -  l')

where l' = l + TFDl .*

find position of crater (rc, Qc) and correct for position angle convention difference

correct for rotational libration, c' = c + TF(DC / 13) : Q'c = 90 - Qc

Q"c = Q'c + C' and find (xc, yc) for crater

find Moon's centre coordinates and translate origin to umbra centre,

x' = x + TFDx,  y' = y + TF Dy,  x = x'  -  xc Sc,  y = y' + yc Sc

where Sc = Scm + TFDSc *

Find approximate contact times for crater,

n2 = (Dx)2 + (Dy)2

t = -1/n2(x Dx + y Dy),   D = + 1/n (x Dy  -  y Dy)

tB = Topp + t  -  1/n (F22  -  D2)0.5

tE = Topp + t + 1/n (F22  -  D2)0.5

Find position angle of contact to check for umbral oblateness later:

                        b = tan-1 (Dy / Dx) where yo = y / cos b

                        YB = b + cos-1 (yo / F2)

                        YE = b - cos-1 (yo / F2)

Notes:

The correction for crater coordinates is approximate

The PA for contacts (for oblateness check) are N to E + ve, N to W - ve.

* indicates these expressions are no longer required as program uses time of opposition in UT

 

Reference

Astronomical Ephemeris Explanatory Supplement, 1961, Section E pages 257 to 262.

 

In option 1 of the Menu the Besselian elements of the Moon and Umbra are found as follows:

Find the Besselian elements of the centre of the Moon with the centre of the umbra as origin:

All data is for times T1 and T2

For the Sun:           a = a0  -  180o  if a0 > 12 hours

                                    a = a0 + 180o  if a0 < 12 hours

                                    d = - d0

For the Moon:        x = (ac  -  a) cos dc

                                                      y = dc  -  d + 1/4 (ac  -  a) sin 2d sin (ac  -  a)

Find the radius of the umbra based on an enlargement of 2%:

                        F2 = 1.02 (p1  -  S'0 + p0),

              with S'0 = S' -  0.000432o , irradiation allowance

                        p1 = 0.998323 pc

and find the Moon's true semi-diameter, Sc = 0.272453 pc +  0.000022o

find rate of change of Besselian elements:

                        Dx = x2  -  x1 and Dy = y2  -  y1

The times T1 and T2 are required to bracket the eclipse time (this was originally computed for each hour, but the program now interpolates Besselian elements at any given time over the eclipse).

Option 6 of the menu incorporates the reduction of single crater timings, where the umbra edge touches the edges of a timed crater. The mathematics were prepared by David Herald and incorporated into my reduction program CTRci.BAS:

 

            Find Q, from cos Q = UM2  -  r2m  -  r2u / 2ru rm

where UM = (x2m + y2m)0.5, rm = (x2 + y2)0.5 and Ru = (x2c + y2c)0.5 ,

            Find F, from cos F = cos (longitude) cos (latitude) of crater

Find a, from tan a = tan Q / cos2F, for 1C and 3C a is + ve,  #

                                    for 2C and 4C a is - ve and is in decimal degrees. 

Find crater radius, rc = crater radius (miles) / 1080 Sc in degrees

Find crater radius, r1 = rc (tan2Q + cos4F)0.5 / (tan2Q + cos2F)0.5  #

Find radius increase, Dru = r1 cos(a - q)  #

Find umbral enlargement, %E = ((ru  +  Dru) -  rt) / rt * 100

Notes:

At primary contacts F = 90o

The scan below shows the variables used in the above mathematics.

Indicate unknown crater radius as - 1 for Prime Contacts, mountains and where crater radius figure is not available.

# Formulae derived by David Herald, Canberra Astronomical Society in February 1982.

 

This scan illustrates the terminology used in the reduction of single crater timings. The mathematics have been incorporated into my QuickBASIC programs.

Oblateness Estimates

Use data produced from the Crater Timing Reduction programs to find statistics for probability of linear regression from:

 

y = a + bx = a + b cos 2y, or r = a + b cos 2y

where              y = equatorial umbral radii (ro, rc)

                        b = linear regression slope  (-  ve)

                        a = linear regression intercept

            y = position angle, 90o at equator

            ro = observed radius

                        rc = computed theoretical radius

                        x = cos 2y

Sxx = n Sx2  -  (Sx)2 , Syy = n Sy2  -  (Sy)2 , Sxy = n Sxy  -  Sx Sy

a = y  -  bx where y = S y / n and x = S x / n , b = Sxy / Sxx ,

Se2 = (Sxx Syy  -  (Sxy)2) / (n (n  -  2) Sxx), and fo = re / (re  -  rp) = 2 b / a

y = (a + bxo) + t a/2 Sc (1/n +  (n (xo  -  x)2 ) / Sxx + 1)0.5

 where probability = 1  -  a

tbo = ( b  -  bo) / Sc (Sxx / n)0.5 with bo = 0,

reject null hypothesis IF tbo < t a/2 or | tbo | > | t a/2 |

r = Sxy / (Sxx Syy)  and Z = (n  -  3)0.5 / 2 LN (1 + r) / (1  -  r) ;

|Z| > | t a/2| and Z > tbo

Use the above expressions to also compute the theoretical umbral oblateness from fc .

 

The mean or expected oblateness by Meeus is

                                     fm = 1/298.3 (p0 + pc) / (p0+ pc  -  S'0)

Notes:

All output for f is reciprocal

Tabulated values for t a/2 can be found in Table IV of Probability and Statistics for Engineers 1965 by Miller and Freud, Chapter 12,

The value of fm from Sky and Telescope, Jean Meeus , April 1979 page 333.

From a private communication from the late Wim Nijenhuis who passed away on August 1, 2000, a best fit ellipse approach has also been applied to crater timing values of re and rp.  This has also been used in my programs in option 8 to demonstrate acceptable observed oblateness.

The program menu also includes option 7, an output menu to allow computations to be output in a form suitable for tabulation. This is done for each observer and a summary table is prepared to list all results and the mean %E for immersion crater timings and for emersion crater timings.

Option 8 includes an oblateness output menu to produce data for tabulation and graphical output for linear regression best fit values and for best fit ellipse tabulation. The data is limited to a larger number of crater timings within the acceptable range of 0 < %E < 4.

Option 9 allows for computation of  19th century crater timings made by Julius Schmidt where data files for the Sun and Moon are in a different form (generated by the USNO Interactive Computer Ephemeris which does not provide polynomial coefficients for the Moon).

Option 10 allows computation of variables for the analysis of images using a topocentric approach. This program generates slant angles and best fit umbral semi-diameter from measures taken of lunar eclipse images during partial phases of the eclipse. It has been found that the umbral semi-diameter is reduced from the expected F2 value due to curvature of the Moon's surface and possible change in size during the eclipse.

Option 11 allows for computation of the Moon's terminator to assist in determining surface features best placed for imaging  - some further development of this program is required.

 

Produced by the Theodore Lunar Observatory  - 2005 July 1, updated 2007 March  14.

 

---