celestial-navigation

Table of contents

1. Introduction
2. Making Sights
   1. Atmospheric Refraction
   2. Dip of Horizon
3. Sight Reduction
   1. Using two sights
   2. Using three or more sights
   3. Running the chicago script
4. Dead Reckoning
   1. Running the sea script
5. A real-life example
6. Terrestrial Navigation
7. Sextant Calibration

1. Introduction

This project contains a toolkit (written in Python) to be used for celestial navigation together with some demos and documentation. Sights (altitude measurements) have to be obtained using a sextant, a nautical almanac and an accurate clock. The toolkit takes care of the sight reduction (conversion to estimated location on Earth), a task that traditionally is performed with careful manual work using special tables (sight reduction tables), pen, dividers and specialized plotting charts. Using a calculator can speed up this manual task, and also reduce possible (human) errors. A computer can do it even quicker, and this toolkit will typically perform a sight reduction in less than one millisecond.
You can use the toolkit on a mobile phone without internet connection. If you also have access to a solar powered battery (powerbank) you have a tool useful while traveling in remote places or on the ocean.

• The script supports stationary observations, i.e. when observations are made from a single position, using multiple sights where the position is determined from the intersection (s) of the small circles of equal altitude. For two sights you will get two possible coordinates. For three or more sights you will get one coordinate (calculated as a mean value).
• There is also support for dead reckoning observations, typically at sea on a moving ship. This also needs a working compass and a chip log or similar. See section 4 below for more information.
• As a bonus there is also support for terrestial navigation. See Section 6 below for more information.

For more information on installation see here.

A digital version of the Nautical Almanac for 2024 is included in this repository.

A more detailed description of celestial navigation can be found here.

If you wonder why I wrote this in the first place, then see this short explanation.

2. Making sights

You create a sight with code like this (for the Sun).

a = Sight (   object_name          = "Sun", \
          time_year            = 2024,\
          time_month           = 5,\
          time_day             = 5,\
          time_hour            = 15, \
          time_minute          = 55, \
          time_second          = 18, \
          gha_time_0_degrees   = 45, \
          gha_time_0_minutes   = 50.4, \
          gha_time_1_degrees   = 60, \
          gha_time_1_minutes   = 50.4, \
          decl_time_0_degrees  = 16, \
          decl_time_0_minutes  = 30.6, \
          decl_time_1_degrees  = 16, \
          decl_time_1_minutes  = 31.3, \
          measured_alt_degrees = 55, \
          measured_alt_minutes = 8, \
          measured_alt_seconds = 0 \
          )

You can also see a complete example in the python script [starfixdata.chicago.py] (starfixdata.chicago.py) and also a corresponding excel file. This sample is built using altitudes taken from a star atlas Stellarium https://en.wikipedia.org/wiki/Stellarium_(software) from a point in central Chicago on May 5th 2024. In other words: No sextant readings were made and the accuracy is very good. (Running this sample will give you an accuracy of just some 100 meters).

The data is picked from your clock, sextant and the Nautical Almanac in the following way

Argument	Description	Remark	Collected From
object	Name of celestial object	Only mnemonic	N/A
time_year	Current year	In UTC	Clock
time_month	Current month	In UTC	Clock
time_day	Current day of month	In UTC	Clock
time_hour	Observation time - Hours (0-23)	In UTC.	Clock
time_minute	Observation time - Minutes (0-59)	In UTC.	Clock
time_second	Observation time - Seconds (0-59)	In UTC.	Clock
gha_time_0_degrees	GHA degrees reading for this hour	For stars use GHA of Aries.	Nautical Almanac
gha_time_0_minutes	GHA minutes reading for this hour	Can be zero (use decimal degrees). For stars use GHA of Aries.	Nautical Almanac
gha_time_1_degrees	GHA degrees reading for next hour	For stars use GHA of Aries.	Nautical Almanac
gha_time_1_minutes	GHA minutes reading for next hour	Can be zero (use decimal degrees). For stars use GHA of Aries.	Nautical Almanac
decl_time_0_degrees	Declination degrees reading for this hour		Nautical Almanac
decl_time_0_minutes	Declination minutes reading for this hour	Can be zero (use decimal degrees)	Nautical Almanac
decl_time_1_degrees	Declination degrees reading for next hour		Nautical Almanac
decl_time_1_minutes	Declination minutes reading for next hour	Can be zero (use decimal degrees)	Nautical Almanac
measured_alt_degrees	Altitude of object in degrees (0-90)		Sextant
measured_alt_minutes	Altitude of object in minutes (0-60)	Can be zero (use decimal degrees)	Sextant
measured_alt_seconds	Altitude of object in seconds (0-60)	Can be zero (use decimal degrees/minutes)	Sextant
sha_diff_degrees	SHA of star vs Aries in degrees	Only use for stars. Otherwise skip	Nautical Almanac
sha_diff_minutes	SHA of star vs Aries in minutes	Only use for stars. Otherwise skip	Nautical Almanac
observer_height	Height of observer above sea level in meters	Only relevant for observations using natural horizon	Height Measurement
sextant	An object defining a specific sextant	See this code sample for details	Calibration

2.i. Atmospheric refraction

The measured altitude values (attributes measured_alt_degrees, measured_alt_minutes and measure_alt_seconds) are corrected for atmospheric refraction using Bennett's empirical formula

$R = \cot \left( h_a + \frac{7.31}{h_a + 4.4} \right)$

Where $R$ is the refraction in arc minutes and $h_a$ is the measured angle from zenith in degrees.

2.ii. Dip of horizon

If you specify the observer_height parameter you will correct for the dip of the horizon. This is useful for observations from a ship deck at sea, or from a hill/mountain with flat surroundings. The dip is calculated using this formula

$a_{\text{diff}}= \arctan{ \frac{d}{r}}$

where

• $a_{\text{diff}}$ is the calculated dip (in radians)
• $d = \sqrt{h \left( 2r + h \right)}$.
• $h$ is the height of the observer (in meters).
• $r$ is the radius of the Earth (in meters).

3. Sight reduction

3.i. Using two sights

Using two star fixes a sight reduction can be done in the following way

from starfix import Sight, SightCollection, get_representation

a = Sight (....Parameters....)
b = Sight (....Parameters....)

collection = SightCollection ([a, b])
intersections = collection.get_intersections ()
print (get_representation(intersections,1))

The result will be a tuple of two coordinates (intersections of two circles of equal altitude). These intersections can be located far away from each other. You will have to decide which one of them is the correct observation point, based on previous knowledge of your location.

The intersections are calculated using an algorithm based on this article
This is a short outline of the algorithm. Two circles $A$ and $B$ define the circles of equal altitude defined from the sighting data as described above. The circles relate to a sight pair $S_{p_{1,2}} = {s_1, s_2}$ which we will come back to later.

$A = \lbrace p \in \mathbb{R}^3 \mid p \cdot a = \cos \alpha \land |p| = 1 \rbrace$
$B = \lbrace p \in \mathbb{R}^3 \mid p \cdot b = \cos \beta \land |p| = 1 \rbrace$

(From now on we assume all coordinates/vectors are located on the unity sphere, i.e. $\lbrace p \in \mathbb{R}^3 \mid |p| = 1 \rbrace$, i.e. the Earth is three-dimensional and its surface has "radius = 1")

We aim for finding the intersections $p_1$ and $p_2$ for the circles $A$ and $B$ and the point $q$ being the midpoint between $p_1$ and $p_2$.

Using the Pythagorean Theorem for a Sphere it is easy to see this:

$\cos aq \cos pq = \cos \alpha$
$\cos bq \cos pq = \cos \beta$

From which we derive this

$q \cdot (a \cos \beta - b \cos \alpha) = 0$

Applying two cross-products we can get the value for $q$

$q = \mathrm{normalize}((a \times b) \times (a \cos \beta - b \cos \alpha))$

Now we can find the intersection points by rotating $q$ for an angle of $\rho$ along a rotation axis $r$.
$\rho$ and $r$ are calculated this way:

$r = (a \times b) \times q$
$\rho = \arccos(p \cdot q) = \arccos\left(\frac{\cos \alpha}{a \cdot q}\right)$

The final rotation is accomplished using Rodrigues/Gauss rotation formula.

$p_{\mathrm{rot}} = q \cos \rho + \left( r \times q \right) \sin \rho + r \left(r \cdot q \right)\left(1 - \cos \rho \right)$

Apply the formula above for $\rho$ and $-\rho$ and you will get the two intersection points $p_1$ and $p_2$.

Note: The algorithm will only work if at least one of the circles is a small circle. It cannot be used for calculating intersections of two great circles.

3.ii. Using three or more sights

Using three (or more) sights a sight reduction can be done in the following way

from starfix import Sight, SightCollection, getRepresentation

a = Sight (....Parameters....)
b = Sight (....Parameters....)
c = Sight (....Parameters....)

collection = SightCollection ([a, b, c]) # Add more sights if needed
intersections = collection.getIntersections ()
print (getRepresentation(intersections,1))

A sight is defined as a collection of data as described in the section 1 above, i.e. recorded data for a celestial object for a specific time.

A collection contains of a set of sights (star fixes) $S$

$S_{\mathrm{sights}} = \lbrace s_1, s_2, \dots s_n \rbrace $

Now set up a set of sight pairs

$S_p = \lbrace {S_p}_{i,j} | i<=n \land j<=n \land j>i+1 \rbrace$

It is easy to see that the number of sight pairs (the cardinality) can be calculated like this

$|S_p| = \frac{n^2 + n}{2} $

For each sight pair we now collect the two corresponding intersection points ($L$ = left, $R$ = right) using the algorithm described in 2.1 above.

$S_{p,i,j} \to \lbrace I_{p,i,j,L},I_{p,i,j,R} \rbrace$

This will result in a set of intersection points

$I_p = \lbrace I_{p,i,j,O} | O \in {L,R } \land i<=n \land j<=n \land j>i+1 \rbrace$

The cardinality can easily be shown as

$|I_p| = n^2 + n $

For each pair of intersection points we calculate the distance. The distance is easily calculated using this formula.

$d\left(x,y\right) = \arccos (x \cdot y) $

This will give us the following set

$D = \lbrace d\left(I_k,I_l\right) | k < |I| \land l < |I| \land k <> l \rbrace$

The cardinality of $D$ can be calculated

$|D| = {|I|}^2 - {|I|} = n^4 + 2n^3 - n$

We now need to eliminate all "false" intersections and only choose those close to the probable location of the observer.

The final part of the algorithm sorts the set $D$ to extract a maximum of $\frac{n^2-n}{2}$ intersection points, and also applying a maximal allowed distance limit (which defaults to 100 km).

The final result will be a single mean value of the extracted intersection points.

3.iii. Running the starfixdata.chicago script

This picture shows the small circles defined in the starfixdata.chicago.py sample

When we move in closer we can clearly see a precise intersection.

Note: The mapping software used in the images above is not precise, and the actual intersections are even "tighter" than shown above.

The output of the script will be like this: First we show the two intersection points from two small circles. The second one is within Chicago.

((N 7°,40.8′;W 94°,14.0′);(N 41°,51.2′;W 87°,38.6′))

Then we add another small circle and show the best three intersections points

BEST COORDINATES
(N 41°,51.2′;W 87°,38.6′)
(N 41°,51.3′;W 87°,38.5′)
(N 41°,51.5′;W 87°,38.6′)

And finally we show the mean value calculation. It is easy to see our values being just 0.1-0.2 nautical miles apart. This accuracy cannot be expected if you use real sextant readings.

MEAN VALUE COORDINATE from multi-point sight data.
(N 41°,51.3′;W 87°,38.6′)

4. Dead Reckoning

When sailing (or moving on the ground) you can use this technique to support dead reckoning where repeated sights (typically of the Sun) will give extra accuracy. You do this by defining a trip segment like this.

from starfix import SightTrip, Sight, LatLon

# We are sailing from point s1 to point s2, in the Baltic Sea.  
# We have a rough estimate of an initial position of 59N;18E to start with
# This estimate is used for selecting the correct intersection point on Earth.
s1LatLon = LatLon (59, 18)

# We define two star fixes  
s1 = Sight (......) # This is your sight at the start of this trip segment. 
s2 = Sight (......) # This is your sight at the end of this trip segment.

See above for how to create a sight object

# We reach s2 by applying about 175 degrees for 1 hour
# (time between taking of sights) with a speed of 20 knots. 
C_COURSE = 175
SPEED = 20    
st = SightTrip (sight_start = s1, sight_end = s2,\
                estimated_starting_point        = s1LatLon,\
                course_degrees                  = C_COURSE,\
                speed_knots                     = SPEED)

Now you can calculate the coordinates for this trip.

intersections = st.get_intersections ()
print ("Starting point = " + str(get_representation(intersections[0],1)))
print ("End point = " + str(get_representation(intersections[1],1)))

The algorithm is a calculation based on distance calculations on segments of the small circles related to $s_1$ and $s_2$.

The two small circles define a sight pair $S_p$ but the two sights are taken at different times. There are two intersection points $\lbrace p_1,p_2 \rbrace$. We select the intersection $p_i$ point which is closest to the estimated starting point $p_e$ by finding the minimum value of $\arccos (p_n \cdot p_e)$

We now define a function on a rotation angle $\rho$ which we will apply on the circle $s_1$

$K(\rho) = \arccos ( (r(s_1,\rho) + t(\phi,\tau))\cdot p_{\text{GPs2}} ) - (\frac{\pi}{2} - \alpha_{s2})$

where

• $r$ is a rotation function (based on Rodrigues formula, see above)
• $t$ is a straight movement function based on approximate course $\phi$ with distance $\tau$.
• $p_{\text{GPs2}}$ is the vector of the ground point of sight $s_2$
• $\alpha_{s2}$ is the altitude measured for $s_2$ (see above for calculation)

The solution of the equation $K(\rho) = 0$ is computed using Newton's method.

The supplied script sample starfixdata.sea.py contains a demo for a short trip at sea (in the Baltic Sea).

This is a picture of the small circles generated by the sample. The larger circle corresponds to the first observation with a lower Sun altitude. The smaller circle is the final observation with a higher Sun altitude.

When we move in closer we can clearly see the intersection $p_i$.

And the course (with our positions) can easily be found using a paralellogram adjustment where we "squeeze in" a route of 20 NM, course 175 degrees, starting at the first small circle and ending at the final circle. The classical method of doing this is of course using a chart and proper plotting equipment and assume linearity of the circle segments.

4.i. Running the starfixdata.sea script

The script outputs the estimated starting and ending points for our trip segment (see the red arrow in the map above)

Starting point = (N 58°,46.1′;E 18°,0.1′)
End point = (N 58°,26.2′;E 18°,3.5′)

In addition we get some diagnostic information. First the radius and GP coordinate of the small circle of the first observation.

S1 radius = 6581.3
S1 GP     = 23.4367,86.7554

And for the final observation.

S2 radius = 5722.4
S2 GP     = 23.4367,71.7571

If you want to plot the trip segment in Google Maps (GM) you have the coordinates here.

Starting point GM = 58.7684,18.0023
Ending   point GM = 58.4364,18.0583

5. A real-life example

You can also see a real-life measurement I recently made using a simple plastic sextant (Davis Mark III), a standard watch and an artifical horizon.

The sample can be found here The resulting position is just 1.45 nautical miles away from my real position, which I consider being an excellent result given the simple equipment and my modest level of training. But there are some question marks about this accuracy, and I will have to make more sights since I need to get more training.

6. Terrestrial Navigation

A sextant can be used for terrestrial navigation too, if you orient it horizontally. Typically you take sights of lighthouses when performing a landfall towards a coast. This sample shows an example of this. The underlying maths are quite similar to sight reduction of star fixes. You need to find the intersection of two circles representing equal angle to two terrestrial points.

The following picture shows how the sample results in two circles of equal angle. The three small circles are centered on three lighthouses and you have measured the observed angle between them from your observation point with a sextant. The red arrow is the calculated correct position.

7. Sextant Calibration

There are many technical aspects of handling and calibrating a sextant and we will not mention all these things here, with one exception. A sextant may show a gradation error which can cause errors for larger measured angles. The little plastic sextant I have used (a Davis Mark III) was suspected by me, and I decided to measure it by taking a terrestrial angle fix and compare it to the data from my map. And yes, there was an error of about 2 minutes / 10 degrees. For the details see this sample where a measurement of a local view is used as input to a calibration parameter of the used sextant.