# Idea

So this time i had the problem of getting the correct vor radial for navigation. To get from one airport to another, i usually set up a route i want to fly by. My thought was, to create a little program that tells you for each point of the route the two nearest vor stations with the radial. With this information it is then possible to triangulate the point of the route by radio navigation.

A quick search revealed openaip.net, an free list of navigation aids that you could download. You get a download for navigation aids per country. The date is presented as .aip xml file.

I then used .gpx as input format for the flight route, since this format is also xml based.

As it is pretty straigt forward within python to work with data in xml format i choose python as the language for this task. I also use the numpy library as well a few more librarys, so my import section looks like this:

The task itself is basically not more then calculating for each point of the route the distance to every vor station, sort the result by distance and drop all but the two nearest ones. Calculate the bearing of thouse two also.

A finished version with a little more functionality to use it from a command line is available at github

The navigation aid files could be downloaded from openaip.net. Note that you need to create an account on that page to actually get to the download link.

# Geometry

## Distance

To calculate the distince of two points on an euclidian surface, one has just to follow the theorem of pythagoras. Calculate square root of the sum of the square of the absolute value of the x coordinate differences and the y coordinate differences. Or for short: \[d = \sqrt{|x_1 - x_2|^2 + |y_1 - y_2|^2}\]

Unfortunatly this does not apply on a sphere like the earth. Since i didn't want to derive that formula on my own, i just look it up. To calculeate the distance between two points on a sphere the haversine formula could be used: \[a = \sin^{2}(\frac{\Delta\phi}{2}) + \cos(\phi_{1}) * \cos(\phi_{2}) * \sin^{2}(\frac{\Delta\lambda}{2}) \]

\[c = 2 * \arctan2(\sqrt(a), \sqrt(1-a)) \]

\[d = c * R \]

- \(R\) is just the mean radius of the earth: ~6.371km
- \(\phi\) is the value of latitude
- \(\lambda\) is the value of longitude

I belive that the intermediate value \(c\) is used, since it is unit free. Only the last step, by multiplying it with the earth's radius introduces the unit. If you multiply it with the radius measured in metres you end up with the distance in meter. The same is true if you multiply \(c\) with the radius in miles.

All right, so the corresponding python code, that calculates the distance looks pretty similar:

Lines 2 and 3 are just the radius in meters and nautical miles. Within Lines 4-7 i calculate the radians for (\phi_{1,2}\), \(\Delta\phi\) and \(\Delta\lambda\). Line 8-10 calculates then \(a\), \(c\) and \(d\) according to the haversine formulas.

## Bearing

For the bearing, there is the same situation. The euclidian formulas do not apply, so i needed to look up the formula for a sphere.

\[ x_1 = \sin(\Delta\lambda) * \cos(\phi2)\] \[ x_2 = \cos(\phi_1) * \sin(\phi_2) - \sin(\phi_1) * \cos(\phi_2) * \cos(\Delta\lambda)\] \[ brg = \arctan2(x_1, x_2) \]

Within the code i also need to first calculate the radian values. Then i calculate the two opperands of the \(\arctan2()\) function and finaly get the result back in degrees. To avoid negative bearings the result is added by 360 and is then divided by flotingpoint modulo 360.

So with these two functions, it is possible to calculate the distance of two points and their relative bearing to each other. From a mathematical standpoint the problem is solved at this point ;-)

# XML Handling

So the basic xml parsing is done by the `xml.etree.ElementTree`

library. One parses a xml file with `ET.parse('filename.xml')`

and
gets back the root tree from which the root is returned by
`getroot()`

. So the navaids of an .aip file are available via:

## aip files

An .aip file consists of a list of navaids. A navaid has a type and various other elements. It looks like the following:

Now there are different types of vor's together with ndb typed navaids
within a .aip file. Since i would like to navigate with vors, i need
to fetch all vor's from the .aip file
(e.g. `openaip_navaid_germany_de.aip`

).

Within an aip file there are vor's of the type: - VOR - DVOR - DVOR-DME - VORTAC

Therefore fetching all vor's is fetching all four vor types and unioning them together:

The main data type for navigation aids used later on consists of bearing, distance, frequency and id. So i wrote getter functions, to read these information out of aip files.

With these getters and the list of navaids from the .aip file, the
list of vor's could be calculated in a format, that is much more
suitable for the given task. A tuple of a *description* and the
latitude/longitude will do the job perfectly. The *description* itself
is a tuple of an id string and a frequency float.

```
((id_string, frequency) , (latitude, longitude))
```

The following maps a lambda expression over the navaid list and greps the data in exactly this format:

## gpx files

The flight path or route is stored within an .gpx file. The xml is pretty basic and the main information is stored within a list of route points:

The file is processed just as the .aip file. First the tree parsed and the root returned:

To get the list of route points only the *rtept*'s of the gpx must be
grepped:

Also getters for the information of the gpx file are needed:

The route is then stored internaly in a tuple of an id string and the latitude/longitude information:

```
(id_string , (latitude, longitude))
```

This list of route points is created by mapping the corresponding lambda expression over the list of route points:

# Calculations along the Route

## Navigation Aids along the Route

The previous sections showed, how the information of navaids and route
points are accessed. This part shows, how to combine them together to
get for every point on the road, the two nearest vor's with their
distance and their bearing. Since we do not know the nearest vor of
each waypoint, we consider every known vor as *along the route*.

The function `calcBrgDst`

calculates just the bearing and the distance
of two points and returns both as tuple:

The `getNavaidsEnroute`

function takes a list of navaids and a
route. It then maps, for each route point, the `calcBrgDst`

function
with the route point over each navaid position.

The result is a list of waypoints, each contianing a list of navaid discription together with bearing and distance.

## Uniqify, Sort and Shorten

### Unique

A sucessful triangulation could only be done with two distinct navigation aids. Therefore it is importend that no navaid is multiple times within the list of navaids along the route. A fast way to unique a list of elements, i found at stackoverflow:

### Sort

A list of everything could be sorted, if one could present a key
function to the library function `sorted`

. All possible key functions
for the navaid data are listed here:

The sort function is then pretty straigt forward:

### Shorten

All could be combined to the function `getNearestNavaidsEnroute`

. This
function takes a list of navaids and a list of waypoints and
optionally a count of how many navaids shold be left in the nearest
list.

It maps the combined `uniqFast(sortDistance())`

function over the
navaids along the route and drops all but `count`

navaids. This
results in the desired list of waypoints, together with a list of a
few nearest navigation aids.

# Showing the Result

Last but not least, a show function helps to visualize the calculation result:

The result of a short trip around stuttgard would result in this output:

```
A_EDSP->STG(116.850): 53°, 15.32NM
A_EDSP->SUL(116.100): 230°, 15.54NM
--
B_N 48->SUL(116.100): 254°, 12.78NM
B_N 48->STG(116.850): 37°, 19.66NM
--
C_N 48->SUL(116.100): 228°, 4.42NM
C_N 48->STG(116.850): 52°, 26.43NM
--
D_N 48->KRH(115.950): 331°, 10.16NM
D_N 48->STG(116.850): 112°, 23.48NM
--
E_EDTQ->LBU(109.200): 56°, 5.45NM
E_EDTQ->STG(116.850): 172°, 10.1NM
--
```