This is intended to give an understanding into the math that allows a navigator with two measurements of the angle of the sun (using, say, a sextant) to calculate his/her position on the earth. I’ve written small python utility that does just this.

## Math

Given two measurements of the angle between the horizon and the sun and two different times of day. Let’s call the angle `theta_1`

and `theta_2`

and the times `tau_1`

and `tau_2`

.

The zenith angle, the compliement of `theta`

, is the angle between the sky’s zenith and the sun. These angles are:

```
phi_1 = pi/2 - theta_1
phi_2 = pi/2 - theta_2
```

If we assume the sun and position vectors are normalized then these angles are the arccosine of the dot product of vector pointing to the sun and the vector pointing to your location.

Using the assumption that your movement on the earth’s surface is insignificant between `tau_1`

and `tau_2`

this gives us

```
cos(pi/2-theta_1) = a1 = sun_1 . pos
cos(pi/2-theta_2) = a2 = sun_2 . pos
```

Fortunately `sun_1`

and `sun_2`

are predictable functions of time, which we know (`tau_1`

and `tau_2`

).

Now we just have to solve for `pos`

as a function of `sun_1`

, `sun_2`

, `a1`

, and `a2`

. To do this we’ll use the vector triple product which states

```
a X (b X c) = b(a dot c) - c(a dot b)
```

Let’s use the above with the vectors of interest for this problem and introduce some new helper variables to make the problem more tractable

```
pos X (sun_1 X sun_2) = sun_1(pos dot sun_2) - sun_2(pos dot sun_1)
pos X (sun_1 X sun_2) = pos X b = -c = sun_1*a2 - sun_2*a1
b X pos = c
```

Figuring out pos is the last step here. With some understanding of cross products and assuming each vector involved here is normalized we can arrive at

```
pos = (c X b)/(b.b) + t*b
pos = d + t*b;
```

Now we can make use of the fact that `||pos||=1`

which means `pos . pos = 1`

.

```
(d+tb).(d+tb) = 1
t^2(b.b) + 2t(b.d) + d.d - 1 = 0
```

And we can solve for t

```
t = [-2(b.d) +- sqrt(4(b.d)^2 - 4(b.b)*(d.d-1))]/(2b.b)
```

And to recap

```
b = sun_1 X sun_2
a1 = cos(pi/2 - theta_1)
a2 = cos(pi/2 - theta_2)
d = (c X b)/(b.b) = (sun_2*a1 - sun_1*a2) X (sun_1 X sun_2) / (b.b)
pos = d + t*b
```

You’ll have two solutions for `t`

and thus two solutions for `pos`

. This will usually mean you need to know which hemisphere you’re in.