- Introduction
- Making Sights
- Sight Reduction
- Dead Reckoning
- A real-life example
- Terrestrial Navigation
- Sextant Calibration
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 toolkit 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 and usage of the Python scripts see here.
A short explanation of the logical steps and algorithms used in this toolkit can be found here.
Digital versions of the Nautical Almanac for 2024 and 2025 are included in this repository.
A more detailed description of celestial navigation can be found here.
An even more detailed description, with a lot of references to various algorithms etc. can be found here
If you wonder why I wrote this in the first place, then see this short explanation.
If you want to contribute to the project then see this page.
You create a sight with code like this (for the Sun). You specify data from your sextant and clock (chronometer). You also add tabular data from the Nautical Almanac. This data is given from the current hour of the observation, and the next hour. (You don't have to enter linear factors etc. from the almanac).
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 a python script 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 | Can be skipped for stars. | Nautical Almanac |
decl_time_1_minutes | Declination minutes reading for next hour | Can be skipped for stars. Can be zero (use decimal degrees) | Nautical Almanac |
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. 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 |
observer_height | Height of observer above sea level or ground 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 |
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
Where
This is a very simple correction.
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
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).
Currently the tool does not aim for very high accuracy and more elaborate corrections are currently elminated.
- Refraction correction is simple (see above), and does not take into account temperature, elevation etc.
- Earth oblateness is ignored.
- Handling of parallax is ignored. It may be needed for Moon sights though, but causes very small errors for all other celestial objects.
Future implementations may contain more precise corrections though. See more about future plans.
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.
For both measurements take note of the measured altitude
(from your sextant),
Using your clock register the corresponding times
Define angles
From the Nautical Almanac, using the timestamps
Now we can define two circles of equal altitude,
The circles relate to a sight pair
(From now on we assume all coordinates/vectors are located on the unity sphere,
i.e.
We aim for finding the intersections
Using the Pythagorean Theorem for a Sphere it is easy to see this:
From which we derive this
Applying two cross-products and a normalization we can get the value for
where
Now we can find the intersection points by rotating
The final rotation is accomplished using Rodrigues/Gauss rotation formula.
Apply the formula above for
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.
Note: I have chosen to use an algorithm based on 3D cartesian vectors. Standard literature on sight reduction typically uses 2D spherical coordinates, such as Longhand haversine sight reduction. Such calculations in 2D are easier to carry out by hand but more complex for a computer. The 3D/cartesian approach is more structurally simple and easier to convert to well-functioning software.
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)
Now set up a set of sight pairs
It is easy to see that the number of sight pairs (the cardinality) can be calculated like this
For each sight pair we now collect the two corresponding intersection points
(
This will result in a set of intersection points
The cardinality can easily be shown as
For each pair of intersection points we calculate the distance. The distance is easily calculated using this formula.
This will give us the following set
The cardinality of
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
The final result will be a single mean value of the extracted intersection points.
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′)
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
The two small circles define a sight pair
We now define a function on a rotation angle
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 geographic position of sight$s_2$ -
$\alpha_{s2}$ is the altitude measured for$s_2$ (see above for calculation)
The solution of the equation
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
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.
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
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.
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 points towards the calculated correct position (intersection).
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.