Given a distance (generally, a large one, say of 850km), a polar coordinate on the earth, and a bearing (with respect to the north pole), I'm using the Haversine formula to calculate a second coordinate. However, as I approach a pole, the shortest distance becomes tangential to that pole (at least I think so–correct me if I'm wrong). I'm not a mathematician, but a programmer, so bear with me while I lay out my formula and the results:
φ2 = arcsin(sin(φ1) * cos(850000/6371000) + cos(φ1) * sin(850000/6371000) * cos(bearing)) * (180 / π)
λ2 = (λ1 + arctan(sin(bearing) * sin(850000/6371000) * cos(φ1), cos(850000/6371000) - sin(φ1) * sin(φ2))) * (λ2 - π) * (180 / π)
I also have a formula that resets the coordinates to 180/-180 degrees when they cross either threshold, but I'm not entirely sure how to represent that mathematically.
Given a starting latitude and longitude (83.8047036171331, -27.8616033804678), I convert to radians and iterate from -180 degrees bearing to 180 degrees bearing, finding coordinates at 850km from my starting coordinate. I'm hoping to have a polygon as an end result. Here are my numbers divided by quadrant with respect to bearing:
-180.0 : 76.1604699668239, -27.8616033804677
-165.0 : 76.2781247730064, -36.2069893102399
-150.0 : 76.6289215559029, -44.5762730745583
-135.0 : 77.2063944341844, -52.9975787393726
-120.0 : 77.9998819868466, -61.5087265207167
-105.0 : 78.994667035189, -70.1659285755627
-90.0 : 80.1720895006312, -79.0598803856636
-75.0 : 81.5094837735314, -88.3495571649968
-60.0 : 82.9795366250911, -98.3431131407785
-45.0 : 84.5477083040733, -109.725327457386
-30.0 : 86.1617477729691, -124.353651321838
-15.0 : 87.6953134517602, -148.975376057487
0.0 : 88.5510627325576, 152.138396619532
15.0 : 87.6953134517602, 93.2521692965516
30.0 : 86.1617477729691, 68.6304445609022
45.0 : 84.5477083040733, 54.0021206964505
60.0 : 82.9795366250911, 42.619906379843
75.0 : 81.5094837735314, 32.6263504040613
90.0 : 80.1720895006312, 23.3366736247281
105.0 : 78.994667035189, 14.4427218146272
120.0 : 77.9998819868466, 5.7855197597812
135.0 : 77.2063944341844, -2.72562802156294
150.0 : 76.6289215559029, -11.1469336863772
165.0 : 76.2781247730064, -19.5162174506955
180.0 : 76.1604699668239, -27.8616033804677
As you can see, once I get to about 85 degrees N, the numbers start to increase precipitously. Now I understand in a very basic sense that the curve of a great circle has to do with its relationship to a pole and that a pole is a kind of coordinate singularity. Now my question is, how do I calculate a point that should, reasonably it seems, transverse a pole? Is the haversine formula inadequate in this regard and if so, is there another methodology for calculating a transpolar location? I apologize in advance if some of my vocab is off! Like I said I'm an amateur with interest in Mathematics but without much of a formal education!
Best Answer
I haven't checked the formula, but you have ${850000 \over 6371000}$ instead of either ${8500000 \over 6371000}$ or ${8500 \over 6371}$.
The point for the -180 bearing in the data above is about 750 km south of the center location. (I believe the correct coordinates are about $(-96.1953, -75.7800)$.)
It is not clear to me how you are computing your new coordinates. The Haversine formula is used to compute the distance given two sets of coordinates and does not have a trivial inverse. Perhaps you can provide a link?
Here is one way to do it, I am sure there must be better ways:
I am performing the computations in the x-y-z plane rather than using latitude/longitude. The mapping between the two is straightforward.
Given a point $p \in \mathbb{R}^3$ on the surface of the earth excluding the poles, a bearing $\beta$ (measured clockwise from north) and a distance $L$, compute a new point $p'$ a distance $L$ away from $p$ at bearing $\beta$ (from $p$).
Note that the poles are excluded so that the bearing makes sense, at either pole a bearing to/from North makes no sense.
I am assuming a perfectly spherical Earth. To reduce notational clutter, define $u: \mathbb{R}^3 \setminus \{0\} \to \mathbb{R}^3$ by $u(x) = {1 \over \|x\|} x$ (that is, normalise $x$).
Let $R= \|p\|$, be the radius (of the idealised Earth), and let $N=(0,0,R)$ be the porth pole.
Note that a distance $L$ corresponds to an angle (at the Earth's center) of $\alpha = {L \over R}$ (radians).
Compute a 'local north direction' $\nu = u(N-{1 \over R^2}\langle p, N \rangle p)$. This is a direction contained in the intersection of the tangent plane at $p$ with the subspace containing $p$ and $N$.
Compute a 'local east direction' $e = u(\nu \times p)$. This is a direction perpendicular to both $\nu$ and $e$ with the appropriate orientation.
(Note that $\nu, p, e$ are mutually orthogonal.)
Compute a 'local direction' $d = (\sin \beta) e + (\cos \beta) \nu$. This is a unit vector in the direction we want to go in on the tangent plane at $p$.
Then we have the new point $p' = (\cos \alpha) p + (\sin \alpha) R d$.