frosch03.de/posts/2015-10-11-little_vor_navigation_helper

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:

{% highlight python lineno %} def distance((lat1, lon1), (lat2, lon2)): r\_m = 6.371e6 r\_nm = 3443.92 phi1 = np.radians(lat1) phi2 = np.radians(lat2) dphi = np.radians(lat2-lat1) dlbd = np.radians(lon2-lon1) a = np.sin(dphi/2) * np.sin(dphi/2) +
np.cos(phi1) * np.cos(phi2) *
np.sin(dlbd/2) * np.sin(dlbd/2) c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1-a)) d = r\_nm * c return d {% endhighlight %}

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

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.

{% highlight python lineno %} def bearing((lat1, lon1), (lat2, lon2)): phi1 = np.radians(lat1) phi2 = np.radians(lat2) dlbd = np.radians(lon2-lon1) x1 = np.sin(dlbd)/np.cos(phi2) x2 = (np.cos(phi1)/np.sin(phi2)) -
(np.sin(phi1)/np.cos(phi2)/np.cos(dlbd)) brg = np.degrees(np.arctan2(x1, x2)) return np.fmod((360 + brg), 360) {% endhighlight %}

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

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

{% highlight xml lineno %} DE BAYREUTH BAY 49.9866676331 11.6383333206 487 110.06 2.47592 FALSE {% endhighlight %}

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

{% highlight python lineno %} navaids\_tree = ET.parse('openaip\_navaid\_germany\_de.aip') navaids\_root = navaids\_tree.getroot() {% endhighlight %}

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:

{% highlight python lineno %} vor\_navaids = navaids\_root[0].findall(".//[@TYPE='VOR']")

• navaids\_root[0].findall(".//[@TYPE='DVOR']")
  
• navaids\_root[0].findall(".//[@TYPE='DVOR-DME']")
  
• navaids\_root[0].findall(".//[@TYPE='VORTAC']") {% endhighlight %}

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.

{% highlight python lineno %} def aipGetLatLon(navaid): lat = float(navaid.find('GEOLOCATION/LAT').text) lon = float(navaid.find('GEOLOCATION/LON').text) return (lat, lon)

def aipGetId(navaid): return (navaid.find('ID').text)

def aipGetFrequency(navaid): return (float(navaid.find('RADIO/FREQUENCY').text)) {% endhighlight %}

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:

{% highlight python lineno %} vors = map(lambda n: ((aipGetId(n), aipGetFrequency(n)), aipGetLatLon(n)), vor\_navaids) {% endhighlight %}

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:

{% highlight xml lineno %} 400.2 A\_EDSP POLTRINGEN {% endhighlight %}

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

{% highlight python lineno %} route\_tree = ET.parse('route.gpx') route\_root = route\_tree.getroot() {% endhighlight %}

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

{% highlight python lineno %} route = route\_root[0].findall("{http://www.topografix.com/GPX/1/1%7Drtept") {% endhighlight %}

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

{% highlight python lineno %} def gpxGetLatLon(pt): lat = float(pt.attrib['lat']) lon = float(pt.attrib['lon']) return (lat, lon)

def gpxGetName(pt): return (pt.findall("{http://www.topografix.com/GPX/1/1%7Dname")[0].text) {% endhighlight %}

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:

{% highlight python lineno %} route = map(lambda pt: (gpxGetName(pt), gpxGetLatLon(pt)), route\_) {% endhighlight %}

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.

{% highlight python lineno %} def getNavaidsEnroute(navaids, route): return map(lambda (sr, r): map(lambda ((sn, fn), n): ((sr + "->" + sn, fn), calcBrgDst(r, n)), navaids), route) {% endhighlight %}

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

{% highlight python lineno %} def calcBrgDst((lat1, lon1), (lat2, lon2)): brng = bearing((lat1, lon1), (lat2, lon2)) dist = distance((lat1, lon1), (lat2, lon2)) return (brng, dist) {% endhighlight %}

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.