Getting around:
In this installation of Calculating Distances by Hand, spherical distances are examined. This will be another function of calculating an arc distance, however, this installation deals with arcs at unknown angles. The basics of the final calculation is similar to the arc distance function of central angle * radius.
This set of formulas will work for things such as calculating the distance a reflection travels across a glass marble, calculating the distance between land masses on earth or even the distance a storm travels across the surface of a faraway planet.
The first installation dealt with 2d and 3d distances along a line, the second installation continued with 2d arcs and circles, this installation deals with 3d arcs and nearly enables us to accurately calculate distances along a planet’s surface. To finalize distance calculations two further shapes will need to be examined; 1) arcs along an ellipse 2) an elliptical spheroid. For this installment, the discussion will be dealing with a Great Circle and associated formulas and ideals within the context of a sphere.
A great circle is a circle that bisects a sphere whose center point is coincidental with the sphere center point.
Getting things rolling;
Building on the last installation of calculating distances along an arc / circle, in this installation I will take a look at another arc distance. This article will begin discussing calculating GPS distances however it should be noted that a Great Circle is the least accurate approximation of a planetary body. Due to the fight between gravity pulling the planet into a sphere and the planet’s spin squashing it into a disk we find balance in an elliptical spheroid.
This article will focus on two different methods of calculating distances along a sphere’s surface; 1) a Great Circle distance 2) Haversine formula.
We will need to define the radius of the planetoid we are interested in. For this example, we will use Earth as our planetoid. Any planetoid would work as well giving us the ability to plot a course on any planet or spheroid’s surface. [ Use Jupiter instead! ]
Due to the spin of the Earth, the planet bulges at the equator which makes our planetary shape an oblate spheroid.
As discussed, the planetoid is not a sphere and in order to use the Great Circle calculations, an approximation needs to be calculated. For this example; mean radius will be used. This will provide us with a radius derived from the WGS-84 ellipsoid. The approximation is defined as: R = 2a+b/3.
a = Equatorial radius (red) (6,378.1370km) = rA
b = Polar radius (cyan & green) (6,356.7523km) = rB
r1 = 2a+b/3 = (yellow) 6,371.0088km (3,958.7613 mi)
A couple of minor functions we will need is unit conversions.
Conv1) Convert Degrees Minutes Seconds to Decimal Degrees [Optional] Conv2) Decimal Degrees to Radians
Conv1: DD° MM’ SS” to Decimal Degrees;
Positive Values == DD = (Seconds/3600) + (Minutes/60) + Degrees
Negative Values == DD = - (Seconds/3600) - (Minutes/60) + Degrees
Conv2: Degrees To Radians;
radians = degrees * π/180
Understanding Latitude and Longitude:
Latitude and Longitude are simply another format in which you can express angles.
Longitude; λ = The vertical lines on a globe; the angular distance of a place east or west of the meridian. In the diagram, the great circles in green and cyan are examples of longitudinal lines.
Latitude; φ = the angular distance of a place north or south of the earth’s equator. In the diagram, the great circle in red is an example of latitudinal lines.
For this example, we will calculate the distance between the coastal areas of:
pA, [Castillo de Salgar of Colombia] and pB, [Hoddevika of Norway].
Technically, considering the great circles in use for this approximation, it is possible to calculate this distance using the information in the previous article discussing arc lengths. As a sphere is in use we can draw a circle whose center point lies at the center of the sphere and has a circumference that passes through pA and pB. Bisecting and aligning a plane to these 3 points will provide for a simple 2d coordinate system. Given the angle between the points pA and pB [δ], along with the radius of the sphere, a distance can be calculated easily along the arc.
Returning to the example on a 3D distance along a sphere’s surface [A 3d Arc]; We will be using the following 7 variables; [in radians]
Note; Abs() = the absolute value of a number, which returns a positive value always.
φA, Latitude degrees of pA = 11.0182000 = 0.1923038676
λA, Longitude degrees of pA = -74.9416610 = -1.30797873136
φB, Latitude degrees of pB = 62.1234567 = 1.084258862132
λB, Longitude degrees of pB = 05.1606150 = 0.0900697231776
φΔ, Equals the Abs( φB – φA) = 51.1052567 = 0.8919549944897
λΔ, Equals the Abs( λB – λA) = 80.102276 = 1.39804845454
δΔ, Equals the angle between pA and pB = [Central Angle] =
ArcCos (sin(φA) * sin(φB) + cos(φA) * cos(φB) * cos(λΔ))
This is the first spherical law of cosines.
cos(φA) = 0.98156652347
sin(φA) = 0.191120799545
cos(φB) = 0.467567964374
sin(φB) = 0.883957124918
cos(λΔ) = 0.171889967964
sin(φA) * sin(φB) = 0.16894259247782760256231
cos(φA) * cos(φB) * cos(λΔ) = 0.07888873943993095754112481693615992
x = 0.24783133191756743930388981693615992
ArcCos(x) = 1.3203552170134638 = δΔ [in radians] = 75.65° +/-
The final step in obtaining an approximated distance along an approximated elliptical spheroid is to simply multiply the central angle between pA and pB [δΔ] times the radius of the spheroid [r1]. Please keep in mind that this approximation is accurate to within 0.5%.
δΔ = ArcCos(x)
d = δΔr1 => [d= δΔ * r1] δΔ * r1 (6371.0088km)
= 8411.99470671868758828144km δΔ * r1 (6371.0088km) * 0.621371
= 5226.96956290849762541802665424mi
In previous articles, the calculations for sin(x) and cos(x) were calculated using infinite number series, power series, and factorials. To calculate for ArcCos another infinite number series will be used called Taylor series and specifically, the subset referred to as Maclaurin series.
A thought about accuracy.
It seems that 0.5% is a small degree of error however in this example 0.5% translates to an error possibility of 26 1/8th miles. Now imagine that you have no method of obtaining a new GPS location and you need to plot the next portion of your total journey. Every consecutive calculation will continue to add errors to both pA and pB of your calculations.
There is a more accurate solution available through Vincenty’s formulae that work on an elliptical spheroid, specifically of the type; oblate. Vincenty’s formulae will return values accurate to within 0.5mm [0.020”].
The great circle and distances less than 1 kilometer+/-:
There is another formula that returns greater accuracy for shorter distances and works slightly better than the great circle distance. This method utilizes the Haversine formula to return the central angle between pA and pB [δ].
The haversine formula is calculated as :
hav(φB – φA) + cos(φA)*cos(φB)*hav(λB – λA)
Where the function hav(θ) = sin2(θ/2) = 1 – cos(θ)/2.
This is referred to as the law of haversines.
Utilizing the above data in the haversine formula;
First, we will solve both hav() functions;
φB = 1.084258862132
λB = 0.0900697231776
φA = 0.1923038676
λA = -1.30797873136
φB - φA = φΔ = 0.891954994532
λB - λA = λΔ = 1.3980484545376
hA = (1-cos(φΔ))/2
hB = (1-cos(λΔ))/2
hA = 0.1860541731415
hB = 0.414055016017
X = hA + cos(φA) * cos(φB) * hB
X = 0.37608433405934156959705105953706226
δΔ = Arcsin(sqrt(x)) = 0.66017761
r1 = 6371.0088km
d = 2δΔr1 => [d = 2 * δΔ * r1]
2 * δΔ = 75.65° +/-
d = 8411.994725745936km
d * 0.621371 = 5226.969574731477998256mi
This value is 0.0000002261% different from calculating the central angle from the law of cosines.
A truer answer is closer to 8414.846978km / 5228.741881566838mi utilizing Vincenty formulae.The Vincenty value shows that the great circle distance is off approximately 0.0338% as predicted [<0.05% error].
Where do we go now?
Let’s now look at a problem in terms of selecting a bearing, or angle and a distance to find a new latitude, longitude pair. In the previous article rotation along the arc was a rudimentary calculation, however, in the case of a sphere and a 3d arc path, a more complex system is needed. For this example, vectors have been chosen to calculate a rotation in 3d space. This will give us a destination point by solving the ‘direct geodetic problem’.
Conv3) Latitude, Longitude => n-Vector
Conv4) n-Vector => Latitude, Longitude
Conv5) Distance => Angular Distance
Conv3) Latitude, Longitude => n-Vector;
|X| | cos(φ) * cos(λ) |
v{x,y,z} = |Y| = | cos(φ) * sin(λ) |
|Z| | Sin(φ) |
Conv4) n-Vector => Latitude, Longitude;
Latitude = φ = atan2(vz, sqrt(vx² + vy²))
Longitude = λ = atan2(vy, vx)
Conv5) Distance => Angular Distance;
Angular Distance = δ = Distance / Radius of great circle
Avg Commercial Flight Speed = 840km/h {14km/min} {0.2333km/sec}
@ 840km/h along a radius of 6371km :
t δ/001 = 0.131847251568699763…radians/hour
t δ/060 = 0.002197454192181166…radians/minute
t δ/360 = 0.000036624236546861…radians/second
Estimated travel time = δ/tδ => distance traveled / rate of travel
Preparing to calculate a vector rotation;
Setting some variables:
vN = {0,0,1} A north pointing vector for a north axis.
vA = {See Above Conv3} Vector representation of start point.
δ = 12345 The angular distance to travel.
φ = 0.2001 Bearing from North [radians].
ê = {computed} East unit vector.
ĵ = {computed} North unit vector.
d = {computed} Directional Vector in the direction of φ.
pB = {computed} Vector representing the destination point.
The formula;
ê = vN * vA
ĵ = vA * ê
d = ĵ * cos(φ) + ê * sin(φ)
pB = vA * cos(δ) + d * sin(δ)
Let’s look at the first step of this formula; calculating the east unit vector:
ê =
[ i j k ] [ 0 * vA.Z + 1 * vA.Y ]
ê = vN * vA = [ 0 0 1 ] = [ 0 * vA.Z + 1 * vA.X ]
[ vA.X vA.Y vA.Z ] [ 0 * vA.Y + 0 * vA.X ]
Let’s look at the next step of this formula; calculating the north unit vector:
ĵ = vA * ê
[ i j k ] [ ê.Y * vA.Z + ê.Z * vA.Y]
ĵ = ê * vA = [ ê.X ê.Y ê.Z ] = [ ê.X * vA.Z + ê.Z * vA.X]
[ vA.X vA.Y vA.Z ] [ ê.X * vA.Y + ê.Y * vA.X]
Let’s look at the next step of this formula;
d = ĵ * cos(φ) + ê * sin(φ)
[ i j k ]
ĵ * cos(φ) = [ ĵ.X*cos(φ) ĵ.Y*cos(φ) ĵ.Z*cos(φ) ]
ê * sin(φ) = [ ê.X*sin(φ) ê.Y*sin(φ) ê.Z*sin(φ) ]
+
d = [ Col. i sum Col. j sum Col. k sum ]
Let’s look at the final step of this formula;
pB = vA * cos(δ) + d * sin(δ)
[ i j k ]
vA * cos(δ) = [ vA.X * cos(δ) vA.Y * cos(δ) vA.Z * cos(δ) ]
d * sin(δ) = [ d.X * sin(δ) d.Y * sin(δ) d.Z * sin(δ) ]
+
pB = [ Col. i sum Col. j sum Col. k sum ]
Utilizing conv4 from above the latitude and longitude can be found by converting from the n-vector (pB)
Why is the Earth an Oblate Spheroid?
The reality is that the planetary body has an odd shape due to the forces acted upon it. There are (3) Three major forces that influence the general shape of the planet.
We live on a spinning planet that orbits our central star which orbits the galactic center of the Milky Way. Our planet travels in a corkscrew,
Firstly, Gravity.
Gravity pulls the planet into a sphere as all points are pulled towards the center of mass. This creates the
Secondly, Centrifugal Forces.
As the planet spins the tendency is to flatten into a disc much like spinning pizza dough. The forces of gravity prevent us from having a flat, disc-shaped planetary body. This creates a bulge at the planet’s equator.
Thirdly, Velocity.
As our solar system is spun around the center of the galaxy the planetary bodies are dragged along by the gravity of our star, the Sun. Imagine if you had a water balloon and you pulled it quickly to the side from the balloon knot. You will notice that the side opposite of the pull becomes the area with the most water (mass). In the same way, the southern hemisphere of our planet also is slightly larger (more mass) due to the same effects motion had on your water balloon.
In the same way that you pulled on the water balloon, our planet is also pulled around by the Sun on its path around the galactic center.
Calculating ArcSin / sin-1 / Inverse Sine:
In order to calculate the inverse value of sine(x), we will utilize another infinite number series referred to as Taylor series. The purpose of this series is to give us the value of functions derivatives at a certain point. The point we will use for our calculations is 0,0. All Taylor series centered at 0,0 are referred to Maclaurin series.
The below Taylor series is for the inverse trigonometric function ArcSine. Let’s begin to understand the layout.
Σ The upper-case Sigma is used as a summation operator.
∞ The infinity symbol above tells us that this will be an
infinite series of polynomial terms.
n=0 The variable below is set to the initial value and will
increase by one for every iteration.
x = sqrt((1-cos(φB - φA))/2 + cos(φA) * cos(φB) * (1-cos(λB - λA))/2) = 0.613257…
The formula to the right of the summation operator Sigma, will be solved using the variables x and n.
Looking at the third iteration of this summation, the polynomial term produced is:
[n=3] ((2 * 3)! / 2^2 * 3 * (3!)^2 * (2 * 3 + 1)) * x^2 * 3 + 1
When simplified; 5x7/112 From above (5(sqrt(x))7)/112
value = 0.00145629826628458304316076327052635345305
[n=5]
value = 0.00008+/-
[n=6]value = 0.00002+/-
As you can see, with every iteration, the value of the polynomial term quickly falls as the math involved to calculate rises exponentially. Even as soon as the fifth iteration this calculation deals with numbers in the millions.
A simplified series can be formulated, the first 8 can be found below.
ArcSin =
x1 + x3/6 + 3x5/40 + 5x7/112 + 35x9/1152 + 63x11/2816 + 231x13/13312 + 143x15/10240 + …
Calculating ArcCos / cos-1 / Inverse Cosine:
ArcCos(x) = π/2 - ArcSin(x)
Calculations for ATan2(y,x):
Considerations before calculating Atan2 need to be accounted for. Looking at the input values of y and x the following formulas apply when using the standard ArcTan function;
If x > 0
then: ATan2(y,x) = ArcTan(y/x)
If x < 0 AND y >= 0
then: ATan2(y,x) = ArcTan(y/x) + π
If x < 0 AND y < 0
then: ATan2(y,x) = ArcTan(y/x) – π
If x = 0 AND y > 0
then: ATan2(y,x) = + π/2
If x = 0 AND y < 0
then: ATan2(y,x) = - π/2
If x = 0 AND y = 0
then: ATan2(y,x) = an undefined value.
The ATan2 function calculates the ArcTangent of the variables y and x with the additional functionality of placing the result in the correct quadrant by taking into consideration the positive/negative values of the variables in use. Arctan calculations have been covered in earlier articles.