Lunar4 is a free celestial navigation sight reduction program and almanac for the Windows desktop. Its almanac function produces barycentric, geocentric, and topocentric coordinates. Its sight reduction function produces azimuth, altitude, and intercept for the Marcq St. Hilaire position line method. It also solves for time, given
The available solar system bodies are the Sun, Moon, and planets, including Pluto. Their positions and velocities are obtained from a Jet Propulsion Laboratory planetary ephemeris. An internal star catalog contains all stars down to third magnitude. There's also a manual entry form for stars not in the catalog.
Lunar4 is copyright © 2019 by Paul S. Hirose. Nonprofit redistribution of the program is permissible if you give me credit.
Lunar4 was created mainly for my use on a Windows 8.1 system. I do not guarantee successful installation in other Windows versions.
The current version is 4.4. To install, download the .zip file (270 kB) and extract its files. By default, when you do "Extract all" a new directory "Lunar4_4" is created in the directory that contains the .zip. To uninstall the program delete its directory. Shortcuts on the desktop and start menu must added by hand. In exchange for that inconvenience you enjoy better security since system administrator privileges are not necessary to install the software. (They may be necessary to install the program in Program Files. I'm not sure since I always install my own programs in a subdirectory of Documents. I recommend you do the same.)
A JPL solar system ephemeris is necessary to operate the program but is not provided with the installation .zip. You can download this DE441 ephemeris (9.0 MB), which is valid for the whole 21st century. It's also possible to build an ephemeris yourself with ASCII files dowloaded from JPL as described later in the Ephemeris Tools section.
Click the Body #1 button. Click Open Ephemeris and select the JPL ephemeris. The default body is the Moon. That's the correct choice for this example. Click OK to exit. Now you're back in the main window. Click the Body #2 button and accept the default choice of the Sun. Next click Topocenter and enter 40N 100W, zero height. The boxes for minutes and seconds of latitude and longitude can be left blank. That's the same as zero. OK out of the dialog.
Click the Time button. Enter 2019-04-27 15:00:00. Click the "leap second table" button, navigate to the directory that contains your Lunar4 installation, and load leapSecTable.txt. After the table loads, the UTC button, which was disabled, is enabled. Select it. In precise work it's necessary to enter UT1-UTC, but for this demonstration the default of zero is close enough. OK out of the time dialog.
Finally, click OK in the main window. Since you have not entered any data in this window, the program runs in "ephemeris mode," where it simply computes and displays the positions of body #1 and #2 at the chosen time. It does not attempt to solve for time. Note that the limb selection drop-downs for the bodies always affect the output, even in ephemeris mode. For example, if you select Moon upper limb, altitude is computed for that limb.
The Tests section has several time solution examples you can try.
The main window appears first. At startup, it forces you to visit the four buttons at the top in sequence before the OK button is enabled. This prevents you from forgetting to set the bodies, position, and time.
Observed angles and limbs may be entered in this window. (They are not required.) Decimal values are OK in the degrees, minutes, and seconds boxes. An empty box is the same as zero. If any component (degrees, minutes, or seconds) is negative, the whole angle is assumed to be negative.
Both bodies have a "center of light" check box. If this is checked, all apparent positions (topocentric and geocentric) of that body are calculated for its center of light instead of center of mass. This box is disabled unless you select "center" for the altitude limb and lunar distance limb, and the body is the Moon or a planet.
Angle formats D (decimal degrees), DM (degrees and decimal minutes), and DMS (degrees, minutes, decimal seconds) may be selected. That affects only output. You can input in any format regardless of the setting.
The precision setting controls the rounding of angle and time to an appropriate display resolution. Sometimes this causes anomalies in the displayed times, when Lunar4 seems to not use the value you entered. Actually it does. But the appropriate rounding unit for time is based on a 1 hour to 15 degrees ratio. So, if 1 arcminute precision is selected, times are rounded to the nearest multiple of 1/15 minute or 4 seconds, which may be a little different from the value you entered. That anomaly affects only the display. Computations are performed with the actual time you enter and a resolution of about one nanosecond.
Near the OK button are four check boxes. If the "solve for time" box is checked, the time solution mode depends on the data:
In all cases the time must be approximately correct to obtain a solution. For a classic "lunar" the position you entered on the topocenter page must be approximately correct too. Time within one hour and position within 10 degrees of latitude and longitude generally assures success.
Note: in a classic lunar solution Body #1 must be the Moon. If another body is selected it is possible to check the "solve for time" box, but the solution will fail. To help avoid mistakes, the Moon is not available in the Body #2 dialog box.
If "solve for time" in not checked, the program simply reports the positions of the bodies. If the observed altitude of a body has been entered, the (observed - computed) difference is displayed, likewise if observed lunar distance was entered.
The "dip correction" check box corrects observed altitudes for dip of the sea horizon with a formula based on the optical principle of the refractive invariant. It's affected by air pressure and temperature at the observer. A temperature lapse rate of -0.009 °C per meter is assumed. This duplicates the Nautical Almanac dip table at 10 C and 1010 mb.
The "solar eclipse corrections" box applies the standard corrections of +.5″ and -.25″ to longitude and latitude of the Moon "to help correct for the difference between the center of figure and the center of mass." (see Astronomical Almanac eclipse section.) Also, the radius of the Moon is changed from the default IAU radius of 1737.4 km to a value such that k = .2725076. This affects the geocentric and topocentric coordinates. With this box checked, contact times for solar eclipses will closely match the values in the Almanac.
The "JPL Horizons compatible" box applies the 1976 IAU precession model and 1980 nutation model. Corrections to the celestial pole computed thereby are published by the IERS, but Lunar4 simulates them by assuming the IAU 2006/00B precession nutation model yields the true pole. (In the years 1950 - 2050 this is true to within a millisecond of arc.) When this button is checked, azimuth, elevation, right ascension, and declination will very closely match of the Horizons values within several centuries of the present. However, it won't be very accurate in extreme epochs.
It's possible to run two or more instances of the program simultaneously. That can be convenient for comparing two different ephemerides, for example.
This appears if you click the Body #1 or Body #2 button at the top of the main window. All controls in this window are disabled until you open a JPL ephemeris with the button at top left. If the program has been used previously, it remembers the last ephemeris and offers that as the default.
The Ephemeris Tools button is explained elsewhere.
Selection of a solar system body is self explanatory. The Moon is only available as Body #1. Otherwise, the windows to select Body #1 and #2 are identical. There are two ways to select a star. You can search for a star designation in the built-in star catalog, which is a subset of the Hipparcos catalog, complete to magnitude 3. This is a substring search, so if you enter "Deneb", the list of results contains Deneb Algedi, Debebola, and Deneb. Select the desired star and OK out of the window.
At the bottom of the window you may enter any real or fictitious star by its catalog data: right ascension, declination, parallax, proper motion, and radial velocity at a specified epoch. That window is explained elsewhere.
All stars of third magnitude or brighter are in the program's internal catalog. (Data are from the Hipparcos catalog new reduction.) Enter a designation in the box and click Search. Generally, if a star has a well known proper name, try that first, then the Beyer designation (with the Greek letter coding and constellation abbreviation in SIMBAD format), then the Flamsteed designation. Many stars have more than one entry. For example, either polaris or alf umi finds Polaris. (The search is not case sensitive.) A substring may be entered. For example, zuben will find Zubenelgenubi. In a few cases there is a suffix in the designation, such as alf02 CVn, so try adding 01 or 02 after the Greek letter if a search fails. Finally, the SIMBAD query by identifiers page may help you guess what identifier I used. Or, examine the lookup table with a text editor. It's in the same directory as the executable.
All search hits are listed in the "search results" box. Magnitude of each star is displayed to assist identification. Select one star and click OK.
Lunar4's star catalog is two separate XML files in the same directory as the executable. One converts a designation, such as "Polaris" or "alf UMi", to Hipparcos catalog number. The other converts Hipparcos catalog number to the catalog data (RA, declination, etc.). These files are XML and can be viewed and modified with any text editor. Instructions are in the files.
Latitude and longitude are accepted in the D, DM, and DMS formats regardless of the format setting in the main window. For example, you can enter integer degrees, fractional minutes, and leave the seconds box blank. The coordinates are assumed to be in the ITRS (practically the same as the WGS84) and height is assumed to be with respect to the ellipsoid.
If the standard temperature and pressure boxes are checked, the values at the observer are assumed to be the ICAO standard atmosphere temperature and pressure at the observer's height (15 C and 1013.25 millibars at sea level). Note that if pressure is entered manually, you must use the altimeter setting, not the actual air pressure at the observer ("station pressure"). I use altimeter setting since it's much easier to obtain. (The "barometric pressure" quoted on TV news is actually the altimeter setting.)
Refraction is affected by temperature more than pressure. It's unusual to experience a 2% change in pressure from nominal, say from 30.0 inches Hg to 29.4. The small temperature rise from 50 to 60 F has the same effect on refraction.
This dialog box is accessed via either body select window. It converts planetary ephemerides from the ASCII format at the JPL download site to the binary format the program needs. See the JPL Planetary and Lunar Ephemerides page for information on the available ephemides. For example, if you want the DE406 ephemeris for years +2000 through +2099, download ascp2000.406. (The "p2000" means +2000 vs. "m2000" for -2000.) Download the DE406 ASCII header file (header.406) too. The files may be put in any convenient location. Select the ASCII to binary option, then choose the input and output files with the buttons at the bottom of the window. For a binary file name I recommend the same name as the ASCII file, except "bin" instead of "asc". If the conversion is successful it runs silently and the dialog box disappears. In recent ephemerides such as DE431, the ASCII file is 300 MB, so the conversion may take several seconds and a "Program Not Responding" message may appear briefly at the top of the window. This is normal.
Some header files have two versions. One is for backward compatibility with old software. That's not an issue for this program, so either header is good. However, there is an important point:
When JPL provides two headers, they have different file extensions, e.g., header.430_229 and header.430_572. But the program requires the header file and ASCII ephemeris file to have the same extension, so you must rename the header.
If the read-only file property is true, two or more instances of Lunar4 can share the same binary ephemeris. The property must be controlled via Windows; you have no access to it from Lunar4.
This dialog box also has functions to merge two binary ephemerides into a third ephemeris, or extract a time subset of an ephemeris. In the merge function, the input files must abut or overlap in time.
Star data may be entered manually in this dialog box. The coordinate system must be the ICRS. Default epoch, 1991.25, is that of the Hipparcos catalog. Other catalogs usually have epoch 2000.0. Right ascension may be entered in hours or degrees. All flavors of sexagesimal are legal, e.g., RA may be entered in decimal degrees with minutes and seconds boxes blank. Proper motion in RA may be milliseconds of arc or milliseconds of time per year, and expressed in coordinate units or great circle units. The minimum usable parallax is 0.1 nanosecond of arc; if the value is less (including zero or negative) the program automatically puts the star at 10 gigaparsecs to avoid numeric difficulty.
Radial velocity has practically no effect unless the star has very high proper motion and parallax. In the preparation of the Hipparcos catalog, the only stars 6th magnitude or brighter whose radial velocity had to be taken into account were mu Cas, e Eri, omicron Eri, alpha Cen, 61 Cyg, and epsilon Ind.
Calendar selection (Gregorian or Julian) is totally manual. Either calendar can be used for any date within the range of the calendar algorithm, which is good for the years -5 million to +5 million. Either the "BC" or algebraic year convention is legal. For example, check the BC box and enter year 1, or do not check the box and enter year 0. Either way, it's the same year.
The time input window also provides for Julian date. For utmost resolution, to read this input the program splits it into two strings: the integer part and the fraction. Both are converted into separate double precision variables.
Time may also be entered as Besselian epoch (e.g., B1950.0) or Julian epoch such as J2000.0.
Next in the time input window are the buttons to select time scale: UTC, UT1, TT (Terrestrial Time), GAT (Greenwich apparent time), LAT (local apparent time), or LMT (local mean time). The UTC button is disabled until you click the Leap Second Table button and load the file.
There is a subtle distinction between GAT and LAT. For GAT, the program computes the geocentric apparent place of the Sun, whereas LAT is based on its topocentric apparent place. That includes the parallax due to the offset from the geocenter, and even the diurnal aberration, but not refraction.
LMT is identical to UT1 plus an offset due to the observer's longitude.
You can manually enter ∆T, but the default is the long term model of Stephenson et al (DOI: 10.1098/rspa.2016.0404). Outside the -719 to +2015 span of the spline, an integration of the length of day expression in the paper is used. This includes the 1500 year oscillation. Biases are applied so there is less than 2 millisecond discontinuity where the integral connects to either end of the spline.
The polar motion x and y angles can be entered. These are never more than a few tenths of a second of arc and for most purposes are negligible.
The angle format seting (D, DM, DMS) in the main window controls output only. Input can use any format because there are separate boxes for degrees (or hours), minutes, and seconds. For example, to enter an angle in decimal degrees, use the degrees box and enter zeros in the minutes and seconds boxes. (If a box is blank, that's identical to zero.) The boxes are wide enough for any normal precision. However, if necessary the number can overflow the box. In that case it won't be possible to see the entire number without scrolling, but the program will read it. Internally, your input is converted to radians and stored in a double precision variable with about 15 significant digits.
Time is stored as a Julian date in two double precision components for high resolution. One part is the JD at midnight. That is, it's an integer plus .5. The other part is the fractional day after midnight. Resolution is about a tenth of a nanosecond.
All four contact times for a total or annular eclipse can be computed. To match the Astronomical Almanac predictions be sure to check the "solar eclipse corrections" box on the main window. That applies standard corrections for the Moon semidiameter and the offset between its center of mass and center of figure. In my tests the contacts were accurate within a tenth second of arc. In many cases you can enter zero lunar distance, set the appropriate limbs (see below), check the "solve for time" box, and let Lunar4 get the time of contact. However, it often happens that the rate of distance is sufficiently nonlinear to defeat the program. In that case you must use trial and error.
To compute the first and fourth contacts, select the near limbs of both bodies. However, to obtain the second and third contact times (automatically or by trial and error) you must understand the Lunar4 definition of near limb: it is the semicircle that is nearest to the center of the other body. In a total eclipse the Moon has greater semidiameter than the Sun, so at the instant of second or third contact its near limb touches the Sun far limb. (To form a mental picture of this, it helps if you imagine the Moon is several times the diameter of the Sun.) On the other hand, in an annular eclipse the Moon far limb touches the Sun near limb at second and third contact.
The classic lunar distance problem is to determine Greenwich time from sextant observations of the separation angle and altitudes of two bodies from an unknown position. To test the Lunar program, create simulated observations with the JPL Horizons online calculator. Time 2019 Jan 1 19:00 UTC. Position 40°N 100°W. Zero height above the ellipsoid.
Apply refraction from the IAU SOFA Library Refco routine. Air pressure 30.3 inches Hg (1026 mb), temperature 0 C, humidity 50%, wavelength .55 microns (about the center of the visual spectrum).
To apply refraction to the semidiameter of the Moon in the direction of the Sun, create a vector to (position angle, -90° + semidiameter), or in this case (67.3, -89.7437). Y-rotate the coordinate system by 90° + the unrefracted altitude of the Moon center to obtain (-.2451, 15.2304). The "longitude" angle is the azimuth (relative to the center of the Moon) of the point on the limb. The latitude angle is its unrefracted altitude.
Add refraction (0.0612) to the latitude angle. Result (-.2451, 15.2916) is the refracted point on the limb that is on the great circle between the centers of the Moon and the Sun. Create a second vector to (0, refracted altitude of center), or in this case, (0, 15.1932), to represent the Moon center. The separation angle between the vectors is the refracted Moon semidiameter in the direction of the Sun: .2561°. Compare that to unrefracted semidiameter .2563.
(The solution is not rigorous. Imagine a clock dial outside the atmosphere with the 12 mark oriented to your zenith. A refracted semidiameter computation for position angle 30°, as described above, would give the semidiameter at the 11 mark. But the dial edge is distorted to a non-circular shape due to refraction, and so the position angle of that mark is not exactlly 30°. But for our purposes the error is negligible, especially since the effect of refraction on semidiameter is so small anyway.)
With the same procedure, compute the Sun's refracted semidiameter in the direction of the Moon: .2710. That's only .0001 less than its unrefracted semidiameter.
From the values calculated above, the near limb to near limb refracted lunar distance = 48.9770 - .2561- .2710 = 48.4499.
Launch the Lunar program and input the simulated observation just created.
In the pre-solution results, the Sun altitude residual (observed altitude minus the predicted value from Lunar) is -8.8866° and the Moon altitude residual is -10.1062°. The lunar distance residual is -0.5599°.
The post-solution residuals in altitude and lunar distance are no greater than .0001, which was the requested accuracy. The solution is 18:59:59.93 UTC, compared to the correct 19:00:00.
Create a simulated observation with the USNO MICA 2.2.2 program at 2019 Jan 10 00:00 UT1, 40°N 100°W, 0 feet above the ellipsoid.
Apply refraction for temperature 0 C, air pressure 1026 mb (30.30 inches Hg), 50% humidity.
Calculate the refracted semidiameter of the Moon. The position angle from the Moon to Aldebaran is 45.8. However, the direction of sunlight is such that the far limb is used for the lunar distance observation, so the position angle of the observed point on the limb = 45.8 + 180 = 225.8. Create a vector to (225.8, -90 + unrefracted SD). Y-rotate the coordinate system 90 + Moon center unrefracted altitude. Apply refraction to that point. Now it represents the refracted vector to the limb point. Create a second vector to (0, refracted altitude of center) to represent the Moon center. The separation angle between the vectors is the refracted semidiameter, .24688, which is .00021 less than the unrefracted SD.
The refracted lunar distance from the far limb to Aldebaran is 98.37954 (the refracted distance from the center, plus refracted semidiameter).
With all data of a simulated observation in hand, get a solution from the Lunar program.
The solution is 23:59:59.880 UT1 vs. the correct 00:00:00.
The simulated Moon - Sun lunar distance observation already computed can also test the determination of time at a known location from a lunar distance only. The result is 19:00:00.10 UTC vs. the correct 19:00:00 if the initial time setting is one hour in error. The program iterates three times.
We can solve a classic time sight with the Aldebaran altitude calculated previously. With the initial time one hour in error the result is perfect at 00:00:00.000 UT1. To achieve this accuracy Lunar iterates three times.
To solve for time, given a lunar distance and two altitudes, the program iteratively seeks a time and position where the observed altitdes and lunar distance are duplicated. When that condition is attained, the computed effects of parallax and refraction on lunar distance virtually duplicate the observation conditions. But the algorithm has difficulty if the bodies have identical or reciprocal azimuths since the observer position becomes indeterminate. For example,
When you solve for time, a message box says the solution does not have the requested accuracy. The problem is the Sun altitude residual (observed minus computed) of .4′. However the altitude residuals are small enough that the solution is satisfactory: 00:00:03.2 UTC vs. the correct 00:00:00.
Here is a traditional time sight with the Sun, sadistically contrived to fall in the middle of the leap second at the end of 2016. Let's see what the Lunar program can do. The first challenge is to generate the problem with JPL Horizons, which does not accept UTC in a leap second. For instance, 2016 December 31 23:59:60 is rejected. Thus we must input time to Horizons in the TT scale. I find it most comfortable to begin half a second after the end of the leap second, at 2017 Jan 1 00:00:00.5 UTC. At that time, TAI - UTC = 37 seconds. Also, there's the constant TT - TAI = 32.184 s. Thus, TT is 1m 9.184s ahead of UTC, so 00:00:00.5 UTC = 00:01:09.684 TT. The center of the leap second is one second earlier at 00:01:08.684 TT.
In Lunar make these settings,
Result is almost exactly in the middle of the leap second as it should be: 23:59:60.492 UTC, an error of 8 milliseconds. If the polar motion angles are not entered, the error reduces to 3 milliseconds, but that's just good luck. The error due to disregarding polar motion partly cancels other errors. In general, for best possible azimuth and altitude accuracy you should apply polar motion.
The accuracy of Lunar has been spot checked against some authoritative references.
Patrick Wallace, Example application of the IAU 2000 resolutions concerning Earth orientation and rotation. In this reduction to apparent place of a fictitious star, the barycentric, geocentric astrometric, GCRS, geocentric apparent intermediate, and topocentric horizontal unrefracted coordinates are less than 1 mas in error.
All tests are highly satisfactory, but it does not follow that the program can duplicate physical reality to that accuracy. What the tests really show is that Lunar4 has a high degree of mathematical consistency with several authoritative references. Imperfect calculation of refraction will be responsible for most of the discrepancies between the real world and the program. To refract the simulated observations I used the SOFA Library Refco routine, which is completely different from the Cassini model in the program. The two models are nevertheless very consistent.
One thing that is not modeled is the departure of the lunar limb from a perfect circle. A high quality occultation program does that. However, a lunar distance observation is very different from an occultation, which "measures" the separation angle with extreme accuracy, the extinction of a star being virtually instantaneous. Also, the point on the limb can be calculated with exactness. Thus, a model of the limb profile makes a noticeable improvement in occultation predictions.
Such refinement cannot be expected from a sextant lunar distance observation. In particular, it's not made at a specific point on the limb. The observer oscillates the sextant about the line of sight to the Moon and seeks a setting such that the star's movement is exactly tangent to the limb. Although the tangent point can be calculated as in an occultation, a graphic of the limb profile reveals rapid variation of height with respect to position angle. This is not noticeable in the hand held low power telescope of a sextant. Thus, semidiameter at the tangent point is unlikely to match the observer's perception. For that reason I believe application of the limb profile would not be beneficial.
The core algorithm was invented independently by me, though I don't claim to be first. It iteratively determines time and position by solving a set of linear equations. There are three observables in the problem: one separation angle and two altitudes. There are three unknowns: time, latitude, and longitude. You supply the observables and estimates of the unknowns. From the estimates Lunar4 predicts the angles that ought be observed at that time and place. In addition, it evaluates how they respond to small changes in each unknown, i.e., computes partial derivatives. This yields nine coefficients in a set of three linear equations in three unknowns. Solving the equations yields adjustments to the unknowns. After applying them, the predicted angles equal to the observed angles. In reality, the assumption that each angle responds linearly to an adjustment in time, latitude, and longitude is not strictly true. Therefore, the process must be repeated until there's no significant change in the solution, to the accuracy requested by the user
In other words, Lunar4 iteratively seeks the time and place that cause its predicted altitudes and lunar distance to equal the observed angles. When that occurs, the parallax, refraction, and all other factors that affect the observations have automatically been applied, as rigorously as they can be computed. The disadvantage of this method is the quantity of computation, which is impractical without electronic assistance.
Between 15° and 10° altitude I use a weighted average of Cassini and the low altitude formula in section B of The Astronomical Almanac. The latter "fades in" linearly, and takes over completely below 10°.
(page last modified 2023-07-04)