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).

Spherical Trigonometry


Image from Wikipedia

In a nutshell, with planar trigonometry, a triangle is described by the lengths of its sides and/or its internal angles. The internal angles always add up to 180°.

With spherical trigonometry, the lengths of the sides are described by their subtended angles instead of distances. The interior angles always add up to more than 180°.

Most navigational problems are really problems in spherical trigonometry, typically describing one triangle side in degrees of latitude and another in degrees of longitude. This makes a right triangle, and we're generally trying to solve for the hypotenuse.

The sight reduction tables are really just a collection of solutions to problems in spherical trigonometry.

A great circle on a sphere is the largest circle you can make that touches the sphere. Its plane passes through the center of the sphere, bisecting it. It has the same radius as the sphere.

A lesser circle is any circle on the surface of a sphere which is not a great circle.

Let a,b,c be the sides of the triangle, as subtended angles.

Let A,B,C be the internal angles of the triangle, opposite the corresponding sides. A,B,C also designate the vertices.

O is the center of the sphere. R is the radius of the sphere.

Cosine Rule:

Cosine rule allows the length of an arc to be determined if the other two arcs and their enclosed angle are known.

  cos(a) = cos(b) * cos(c)  +  sin(b) * sin(c) * cos(A)
  cos(b) = cos(a) * cos(c)  +  sin(a) * sin(c) * cos(B)
  cos(c) = cos(a) * cos(b)  +  sin(a) * sin(b) * cos(C)

Sine rule:

Sine rule allows you to find an angle if two sides and an angle are known, OR to find a side if two angles and a side are known.

  sin(a)   sin(b)   sin(c)
  ------ = ------ = ------
  sin(A)   sin(B)   sin(C)

Spherical trig is often used in navigation. In such cases, point 'A' is typically the north pole and 'B' and 'C' are two navigation points.

Example: distance and direction from London and Baghdad

  North Pole:   90.00°N,  0.00°    A
  London:       51.30°N,  0.10°W   B
  Baghdad:      33.20°N, 44.26°E   C

Draw edges c,b from the North Pole to these two cities. The lengths of the edges are derived from the latitudes and are 38.7° and 56.8°. The angle between these edges, A, is derived from the longitudes, and is 44.36°.

The distance is edge a, which we can get from:

    cos(a) = cos(b) * cos(c)  +  sin(b) * sin(c) * cos(A)

Giving 36.74° between them. Converting to radians and multiplying by R gives 2564.73 miles. (Or compute circumference * 36.74/360 to get the same result.)

Note that sin(x) = cos(90-x). We can use this to apply the cosine rule to latitudes directly:

    cos(a) = sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) * cos(dlon)

The direction to Baghdad from London is angle B. Since we have all the sides and angle A, the sine rule can give us B:

             sin(A)
    sin(B) = ------ * sin(b)
             sin(a)

    sin(B) = 0.98

    B = 77.96 or 102.04.

Correct answer is 102.04. I'm not sure how you choose between the different answers other than by looking at the diagram.

References: