Find the intersection point of a great circle arc and latitude line

In spherical geometry, I need to know at what longitude λ a great circle arc φ1,λ1-φ2,λ2 has intersected a line of latitude φ.

I have found the equivalent equation for solving latitude φ for an unknown longitude λ on the excellent Williams Aviation Formulary, excerpted here:

It would seem all that is required is to rearrange this equation in terms of lon rather than lat, but performing the rearrangement is beyond me (and, it would seem, Wolfram Alpha).

λ = ?

Edit: Williams Aviation Formulary also explains 'Crossing Parallels' which seems to be what I need, but this finds the crossing points for a whole great circle rather than a great-circle arc. What is the most effective way of finding only the crossing of the arc?

Many thanks for any advise you can lend on this.

It's a somewhat complex problem, unfortunately. I found using vectors easier than spherical geometry.

For the intersection point $\vec{I}$ , you know the latitude and therefore the z-component. You also know it follows the great circle path defined by $\vec{GC} = \frac{\vec{P_1} \times \vec{P_2}}{||\vec{P_1} \times \vec{P_2}||}$.

To find an unknown point along a great circle defined by vector $\vec{GC}$ you can use the vector addition:

$\vec{P1} \cdot \cos{\alpha} + (\vec{GC} \times \vec{P1}) \cdot \sin{\alpha} = \vec{I}$

You can solve for $\alpha$ since you know the z-component of all three vectors.

$(\vec{I})_z = \vec{P}_{1z}\cos{\alpha} + (\vec{P}_{1x}\vec{GC}_y - \vec{P}_{1y}\vec{GC}_x) \sin{\alpha}$

This page nicely shows you how to solve for an equation of form $a \cos{\theta} + b \sin{\theta} = c$

The arc distance then resolves to:

$\alpha = \tan^{-1}\left ( \frac{\vec{P}_{1x}\vec{GC}_y - \vec{P}_{1y}\vec{GC}_x}{\vec{P}_{1z}} \right) - \cos^{-1} \left ( \frac{\vec{P}_{2z}}{\sqrt{(P_{1z}^2 + (\vec{P}_{1x}\vec{GC}_y - \vec{P}_{1y}\vec{GC}_x)^2}} \right )$

And now it's the computer's turn to do the work ;)

EDIT

As this was something I had to do recently for a project, I have some code in R. The last function (LatIntersect) returns the angle in radians of the arc til the desired latitude. Multiply by $R_{earth}$ to get distance, then finding the coordinates should be trivial:

Coords = function(lat, lon)
{
x = cosd(lon) * cosd(lat)
y = sind(lon) * cosd(lat)
z = sind(lat)
return(c(x = x, y = y, z = z))
}

VecCross = function(V1, V2)
{
if(length(V1) != 3 | length(V2) != 3)
stop("Wrong vector size")

x = V1["y"] * V2["z"] - V2["y"] * V1["z"]
y = V1["z"] * V2["x"] - V2["z"] * V1["x"]
z = V1["x"] * V2["y"] - V2["x"] * V1["y"]
x = as.numeric(x)
y = as.numeric(y)
z = as.numeric(z)

return(c(x = x, y = y, z = z))
}

LatIntersect = function(lat1, lon1, lat2, lon2, LonIntercept)
{
Start = Coords(lat = lat1, lon = lon1)
End = Coords(lat = lat2, lon = lon2)

GC = VecCross(Start, End)
N = sum(GC * GC)
GC = GC / sqrt(N)

zProj = sind(LonIntercept)

a = Start["z"]
b = Start["y"] * GC["x"] - Start["x"] * GC["y"]
c = zProj

arc = atan(b/a) - acos(c / sqrt(a^2 + b^2))

return(arc)

}

• Fantastic, thank you for the time and effort to produce this answer. I'll update this thread with an attempt at turning your final formula into executable code, as soon as I'm able. – Chris Hatton Oct 30 '15 at 2:49

I would rewrite the first equation as $\tan \phi = ...$ Then manipulate it to get an expression of the form $k = a \sin(\lambda-c) + b \sin (\lambda-d)$ where $a$, $b$, $c$, $d$ and $k$ do not depend on $\lambda$.

Wolfram Alpha provides solutions to the latter equation.