Spherical Trigonometry Related to Distance Computation 

3. SPHERICAL TRIGONOMETRY FORMULAS 3.1 Notations 3.2 Length computation 4. INTERSECTION OF TWO ARCS OF GREAT CIRCLES 4.1 Principle 4.1.1 Assumptions 4.1.2 Principle 4.2 Calculus 4.2.1 Check null arcs 4.2.2 Converting geographical location from degrees to radians 4.2.3 Converting geographical location into Cartesian system 4.2.4 Equation of a plane 4.2.5 Coordinates of great circles crossing points on the earth 4.2.6 Converting Cartesian system into geographical location 4.2.7 Check location of crossing points with arcs of great circles 4.3 Summary
This paper describes
the content of the current procedures developed by This document is given by EUROCONTROL for information purposes only. Its content is NOT the final version that will be used for operational purposes. EUROCONTROL cannot accept any responsibility due to a possible misuse of this text. The content of this document can be copied freely, provided this disclaimer notice is accepted, and that the information contained is not sold or included into a product developed for commercial purposes. EUROCONTROL is the European Organisation for the Safety of Air Navigation DED is The Directorate of EATCHIP Implementation CRCO is the Central Route Charges Office EATCHIP is the European Air Traffic Control Harmonization and Integration Program. The basic classical formulas are based on a spherical triangle, which is a figure on the surface of the sphere enclosed by three arcs of great circles (see Fig. 1.) Notation: in the triangle ABC, the angle ABC=B, angle BAC=A, angle BCA=C and side AC=b, side CB=a and AB=c. We will also call the points at the three vertices of the triangle A, B and C (‘A’ will represent the measure of the angle BAC as well as the name of the point). The same convention will be used for sides: for instance ‘a’ will be the name of the side BC and will also represent the measure of the orthodromic length between point B and C. The length of each side a,b,c is expressed either in nautical miles (nm) or as an angle (because a, b, c are great circles 1 nm = 1 minute of angle).
There are two main formulas expressing relations between angles and sides in a spherical triangle:
For computing purposes, we can set the top point of the spherical triangle at the North pole and call it P, then the sides a and b can be expressed in radians coming directly from the colatitude of points B and A. (see Fig.2.) Angle P (at the North Pole) can be expressed as a difference of longitude between point A and B. Side m is the length between points located at angles A and B. ZA is the azimuth at A of the side m with the meridian AP. Note 1: a colatitude is the angle between North Pole and the latitude. It is:
Note 2: we can use both latitude or colatitude in formulas, assuming that:
The orthodromic length between two points can be found directly by using a cosine formula. For instance, searching for m by using Fig. 2 we can express:
The following example is a function programmed in "C" data processing language. It processes a distance expressed in nautical miles from two coordinates given in parameter and expressed in minutes of angle. This distance is always positive and will never exceed half of the circumference of the earth which is 10800 nm. Notes:
#define PI 3.14159265358979323846 double distance(double l1, double g1, double l2, double g2) /* expressed in minutes */ { double lati1 = PI*l1/10800.0; /* conversion in radian ... */ double lngi1 = PI*g1/10800.0; double lati2 = PI*l2/10800.0; double lngi2 = PI*g2/10800.0; double v = sin(lati1)*sin(lati2) + cos(lati1)*cos(lati2)*cos(lngi1lngi2); if (v >= 1.0) return 0.0; if (v <= 1.0) return 10800.0; return (10800.0 * acos(v) / PI); /* distance is expressed in nautical mile */ } This part of the document describes the content of the current function defined by DED4 Section 5 and applied by the CRCO to calculate the intersection of two arcs of great circles. One arc is a segment on the route of the aircraft, the other is a FIR boundary segment. This part is subdivided in three sections:
The shortest path between two points on the surface of the spheroid is along the arc of a geodesic curve. On the surface of a sphere the geodesic curves are the great circles. We suppose that the earth is a sphere. We suppose that aircraft are always looking for the shortest path, and then fly along an arc of great circle between two route points. It also means that the longest part of a great circle is always less or equal to half of the perimeter of the earth. The purpose is to find the intersection point between two arcs of a great circle, that can represent for instance a segment on the route of the aircraft, and a segment of the FIR boundary. The calculations are made in Cartesian coordinates and use classical geometry in third dimension (vectors, line, planes). The arcs are along great circles which, by definition, are a plane cutting of the earth by its center. Two arcs of great circles generate two great circles, which generate two planes which both cut the earth by its center. We already note that these two planes can only:
If the two planes intersect, it is along a line in space which crosses the surface of the earth at two points (exactly opposite to each other). We must check if one of these two points belongs to both the arcs of great circles. if yes, we have a solution, otherwise we have not. If the two planes are identical, it means that the two arcs of great circle are aligned. In this case, even if the arcs have a common part there is no solution (in fact this case is handled outside this function in the context of the CRCO/RSO project) Fig. 1
The different steps of the calculus are identified within different paragraphs. The logical succession of these steps are respected to the case where we found a solution (the two arcs cross). See the summary at the end of the document to get an overview. The function, generally called "iagc" for "Intersection d’Arc de Grands Cercles" (in French) will return a code and may be a result. The different cases are detailed below: Table 1
Before starting any calculation, we check that the two arcs have a nonnull length. In case of null value, there is no solution. By making the assumption that geographical coordinates are expressed in degrees with decimals, the transformation into radian is:
The two arcs of a great circle are represented by four points P1, P2, P3 and P4 located on the earth (two points for each arc). These points are expressed in latitude and longitude (that must be transformed in radian before any calculation). We must convert them into x,y,z coordinates. Knowing the equation of a sphere having a radius of one, we have for the first point P1(lat1,lon1):
Note: The center of the earth is at the origin and OZ is pointing "upward" to the North pole. The generic equation of a plane containing the origin of axes (which is a plane crossing the center of the earth) is:
This type of plane is created from only two points which are the extremities of an arc of great circle (default third point is the center of the earth). So, from P1 and P2 (extremities of the first arc) we create the vector V1 (vector director of the first plane, perpendicular to it) and from P3 and P4 (extremities of the second arc) we calculate V2 (vector director of the second plane, perpendicular to it). V1 and V2 are in fact results of a vectorial product. The expression of V1 is:
In order to facilitate the test to identify if the two planes are identical, we calculate the unit vector U1 and U2 from V1 and V2. For instance the calculus for U1 is:
The test checks if the two vectors U1 and U2 are identical. Note the use of epsilon which is positive and very small.
If the two planes are different, they cross along a line in space which is represented by the following equation:
This line has a vector director D which is the result of the vectorial product of the two vector director perpendicular to the planes:
Because the sphere representing the earth has a radius of one, again we calculate the unit vector of D, so that it touches the surface of the sphere. This vector S1 and its opposite S2 directly give the coordinates of the crossing point of the two nonoverlapping great circles on the sphere.
We know two ways of doing it. The first method, presented here, is the exact inverse of the procedure used in the paragraph "Converting geographical location from degrees to radians". The second one, which may be longer, uses calculation of angle from projected length with a rotation around the origin. The first method processes the latitude from the z coordinate (remember OZ axe is pointing "upward"). We use arc sine which returns an angle in the range from PI/2 to PI/2 which corresponds respectively to South Pole and North Pole:
The longitude is calculated from, for instance, x coordinate using arc cosine. Because arc cosine returns an angle only from 0 to PI, we must check the sign of the angle return from arc sine calculated from y coordinate. If the sign is negative the angle for longitude must be in the range from 0 to PI (also negative).
Once S1 and S2 are known, we must make sure, in order to have a solution, that the location of one of them is present on both arcs of a great circle. The check of the location of S1 and S2 is simply based on orthodromic length, calculated with geographical coordinates, between extremities of arcs and S1 or S2.
L1t is the total length of the first arc, while l1a is the length between one extremity of the arc to S1, and l1b the length from the other extremity to S1. "Test1" equals 0 means that S1 is somewhere on the first arc. To check if S1 is also on the second arc let’s calculate "test2":
If S1 is located on the first arc and also on the second arc (test1=test2=0), S1 is the solution. If not, the same check must done with the other point S2 (the opposite of S1). If none of the two checks are positive, there is no solution: arcs do not intersect. Improvements (in respect of accuracy and processing time) are currently being studied, and will be implemented in future versions.
I got a good headache after reading it too, good stuff though.

