Previous | Home | Next

Math for Celestial Navigation

With an electronic calculator and a few simple formulae, you can replace several volumes of tables.

For many more on-line calculators, see the National Geopatial-Intelligence Agency which has a full suite of javascript calculators which you can download and use locally on your computer.

Dead Reckoning

Dead Reckoning is simply an exercise in vector addition. You represent your current location as an X,Y point (longitude, latitude) and your velocity as an X,Y vector. Multiply the vector by time and add to the initial point. If a current is involved, that's merely another vector to add to your calculations.

The problem is made more complex by the fact that you are traveling across the surface of a sphere. This means your velocity vector is actually a curve (a segment of a great circle) rather than a straight line. True dead reckoning calculations are actually exercises in spherical trigonometry.

However, for short distances not near the poles, you can treat the Earth's surface as flat and your velocity vector as a straight line.

In the simplest case, you make a rectangular projection of the Earth's surface, where one minute of latitude equals one nautical mile, and one minute of longitude equals cos(latitude) nautical miles.

Dead Reckoning, rectangular projection

Here is a simple Javascript calculator for you to play with.

Start Latitude: degrees
Start Longitude: degrees
Speed: kts
True Heading: degrees
Time: hours
Latitude:0
Longitude:0

Treat north latitude as positive numbers and south latitude as negative numbers. Treat east longitude as positive numbers and west longitude as negative numbers. Use degrees for all calculations instead of degrees and minutes.

Given:

Bs = Latitude at start, degrees
Ls = Longitude at start, degrees
V = speed, kts
hdg = heading
T = time, hours

Compute:

Vx = V ˇ sin(hdg)
Vy = V ˇ cos(hdg)
B = Bs + T ˇ Vy
L = Ls + T ˇ Vx / cos(B)

If the difference between Bs and B is significant, you might want to average the two in that fourth equation.

To compute current, just repeat the above steps, using the current speed and direction.

Sight Reduction

Given Assumed Position, AP (lat, lon) and Geographic Position of a body, GP (declination, gha), the computed altitude (Hc) is found by:

Here is a simple Javascript calculator for you to play with.

Latitude: Hc: 0
Declination: Z: 0
LHA:
lha = gha - lon sin(Hc) = sin(lat) * sin(dec) + cos(lat) * cos(dec) * cos(lha)

and the azimuth (Z) is found by:

sin(dec) - sin(Hc) * sin(lat) cos(Z) = ----------------------------- cos(Hc) * cos(lat) if( lha <= 180 ) Z = 360 - Z

Publication Ho229 gives these formulae:

sin(Hc) = sin(lat) * sin(dec) + cos(lat) * cos(dec) * cos(lha) cos(dec) * sin(LHA) tan(Z) = ---------------------------------------------------- cos(lat) * sin(dec) - sin(lat) * cos(dec) * cos(lha)

Dip and distance to horizon

If you search for references on computing distance to the horizon or sextant dip caused by the Earth's curvature, you'll probably find some excellent references on the mathematics involved, including derivations. (Image from Wikipedia article).

Example:

distance to horizon = d = √[2 r h]
dip (minutes) = 0.531 √h (feet)

However, these calculations don't account for the refraction of the Earth's atmosphere, and therefore cannot be used. Bowditch gives this formula which approximately compensates for refraction:

ro = radius of earth, 3440.1 nm
hf = height of eye, feet
Βo = refraction factor = 0.8279
D = distance to horizon in nm

D = √[(2rohf)/(6076.1Βo)]

Simplified:

D = 1.169 √hf    (nautical miles)
D = 1.345 √hf    (statute miles)

Refraction

Atmospheric refraction is zero at the zenith, 1' at 45° altitude, 5' at 10° altitude, and 34' at the horizon. Since the diameter of the Sun is about 30', this means that when you see the Sun touching the horizon at sunset, in reality it has already set.

Refraction is affected by temperature and pressure. Add 1% to the refraction for every 3°C colder than 10°C or subtact if hotter. Add 1% for every 0.9kPa higher pressure; subtract if lower.

Calculating refraction in the Earth's atmosphere is too complicated to cover here, since it involves having a mathematical model of the atmosphere, and the use of differential equations. See this web page at SDSU for more information.

For altitudes above 30°, a reasonable approximation is

r ≈ k tan z

Where z is the zenith angle (i.e. zero at the zenith = 90°-ht), and k is about 60.3". k actually varies with wavelength, being 60.4" at blue, and 57.2" at red.

For altitudes less than 30° above the horizon, you really need to use a published table.

Apparent altitude, happ is true altitude, ht plus refraction. This equation allows you to build a table that converts apparent altitude to true altitude.

Calculated Position

See Celestial Fix Nautical Almanac Sight Reduction algorithm for n LoPs by Andrés Ruiz San Sebastián-Donostia for more information on this section.

On paper, you determine your position based on two or more fixes. For each fix, you compute a line of position. Where the lines of position intersect, that's your location.

But what if you have more than two lines of position? It would be almost a miracle if all three intersected at a single point. In practice, you'll get an intersection for every possible pair of lines of position. With a little luck, the intersections will all cluster around a small region, which you can take as your fix.

One possible approach to solving this dilemma is to compute all the possible intersections and take an average position. An improved algorithm is the least squares algorithrm.

Definition: p is the difference between Ho and Hc, (p = Ho-Hc). Z is the azimuth of the sight. If difference is positive, the line of position is drawn along the azimuth, at a distance of p. If p is negative, the line of position is drawn in the opposite direction.

Given a series of p and Z values, pi and Zi, we compute A, B, C, D, E, and F as follows:

A = Σ cos² Zi

B = Σ cos Zi - sin Zi

C = Σ sin² Zi

D = Σ pi cos Zi

E = Σ pi sin Zi

F = Σ pi²

And:

G = A C - B²

Check: A+C = n

(Note: The symbols for latitude and longitude are B and L. Try not to confuse the B term above with the B symbol for latitude.)

Given the estimated position Be, Le, the improved fix is given by:

dB = (CD - BE)/G

dL = (AE - BD)/(G cos Be)

Bv = Be + dB

Lv = Le + dL

The estimated error (standard deviation) is:

σ = 60 √[(F - D dB - E dL cos Be)/(n-2)]

σB = σ √[C/G]

σL = σ √[A/G]

If you want, you can substitute the new fix position for the assumed position and start over from the top for an improved estimate.

Confidence Ellipse for axes a,b is:

Prob = e.g. 0.95 for 95% confidence.

k = √[-2Ln(1-Prob)]

tan 2θ = 2B/(A-C)

a = σ k / √[n/2 + B/sin 2θ]

b = σ k / √[n/2 - B/sin 2θ]

See the cited paper for more details and a sample implementation in C.

Moon's semidiameter

This can be found in the almanac, but it can also be easily calculated from the Horizontal Parallax.

SD = 0.2724*HP

This follows simply from the fact that the diameter of the moon is 0.2724 that of the Earth.

Augmentation of Moon's semidiameter

When the Moon is directly overhead, it's about 4000 miles closer than when it's on the horizon. This makes it appear larger, and the effect is known as the "augmentation" of the semidiamter. You can calculate this as 0.3'*sin(H_moon).