interpolating    mapping    searching

Spheres: Great Circle Intersections.

Great Circle intersects Great Circle: Any two points on a sphere uniquely define a great circle (so long as they are not antipodal). Given two points (lat₁, lon₁), (lat₂, lon₂) the great circle that contains these two points is defined by points (lat_c, lon_c) that satisfy:

lat_c= atan((sin(lat₁)*cos(lat₂)*sin(lon_c - lon₂) - sin(lat₂)*cos(lat₁)*sin(lon_c - lon₁)) / (cos(lat₁)*cos(lat₂)*sin(lon₁ - lon₂)))

Any two Great Circles intersect exactly twice.  So given four points defining two great circles we seek the two points at which those great circles intersect.  Two equations and two unknowns, it should be simple.  But if we do it on the sphere it involves spherical trig, and spherical trig involves a lot of trig.  That is:

Given:

lat_c= atan((sin(lat₁)*cos(lat₂)*sin(lon_c - lon₂) - sin(lat₂)*cos(lat₁)*sin(lon_c - lon₁)) / (cos(lat₁)*cos(lat₂)*sin(lon₁ - lon₂) ) )

and

lat_c= atan((sin(lat₃)*cos(lat₄)*sin(lon_c - lon₄) - sin(lat₄)*cos(lat₃)*sin(lon_c - lon₃)) / (cos(lat₃)*cos(lat₄)*sin(lon₃ - lon₄)))

Solve for (lat_c, lon_c)

Here is one solution for lon_c.  This is one fourth of the complete solution:

lon_c = -ArcCos[ -Sqrt( b² f² Cos[lon₁]² - 2abf² Cos[lon₁] Cos[lon₂] + a²f² Cos[lon₂]² - 2bcef Cos[lon₁] Cos[lon₃] + 2acef Cos[lon₂] Cos[lon₃] + c²e² Cos[lon₃]² + 2bcdf Cos[lon₁] Cos[lon₄] - 2acdf Cos[lon₂] Cos[lon₄] - 2c²de Cos[lon₃] Cos[lon₄] + c²d² Cos[lon₄]²)/ Sqrt( b²f² Cos[lon₁]² - 2abf² Cos[lon₁] Cos[lon₂] + a²f² Cos[lon₂]² - 2bcef Cos[lon₁] Cos[lon₃] + 2acef Cos[lon₂] Cos[lon₃] + c²e² Cos[lon₃]² + 2bcdf Cos[lon₁] Cos[lon₄] - 2acdf Cos[lon₂] Cos[lon₄] - 2c²de Cos[lon₃] Cos[lon₄] + c²d² Cos[lon₄]² + b²f² Sin[lon₁]² - 2abf² Sin[lon₁] Sin[lon₂] + a²f² Sin[lon₂]² - 2bcef Sin[lon₁] Sin[lon₃] + 2acef Sin[lon₂] Sin[lon₃] + c²e² Sin[lon₃]² + 2bcdf Sin[lon₁] Sin[lon₄] - 2acdf Sin[lon₂] Sin[lon₄] - 2c²de Sin[lon₃] Sin[lon₄] + c²d² Sin[lon₄]²)]

It is much easier in Cartesian space.  While any two non antipodal points on a sphere define a great circle any three non colinear points in space define a plane.  The two points that define the great circle, combined with the origin, can be used to uniquely define the plane that contains that great circle.  So given two planes we can find the line defined by the intersection of those planes, and then find the two points where the line intersects the sphere.  It's a lot more steps, but a lot less trig and, overall, a lot less math.

Given two planes (that contain the origin) and a sphere:

(y₀z₁- y₁z₀)x - (x₀z₁- x₁z₀)y + (x₀y₁- x₁y₀)z = 0

x²+ y²+ z²= r²

Solve for (x, y, z)

The algorithm used by the spheres package is as follows:

1. Define the constants in the planar equations. 

Let a = (y₀z₁ - y₁z₀ )

Let b = -(x₀z₁ - x₁z₀ )

Let c = (x₀y₁ - x₁y₀ )

Let d = (y₂z₃ - y₃z₂ )

Let e = -(x₂z₃ - x₃z₂ )

Let f =  (x₂y ₃ - x₃y ₂)



So the planar equations simplify to:
ax + by + cz = 0
dx + ey + fz = 0


2. Solve for x:

ax + by + cz = 0 so
x = (-by - cz) / a


3. Solve for y:
   

dx + ey + fz = 0
d(((-by - cz) / a)) + ey + fz = 0
(-dby/a) - (dcz/a) + ey + fz = 0
ey - (dby/a) = (dcz/a) - fz
y(e-(db/a)) = z((dc/a) - f)
y((ea - db)/a) = z((dc-fa)/a)
y = z * ((dc-fa)/a) * (a/(ea-db))
y = z * ((dc-fa)/(ea-db))

Let h = (dc-fa)/(ea-db) so y = hz



4. Get x in terms of z:
   
x = (-by - cz)/a
x = (-bhz - cz)/a
x = z * ((-bh - c)/a)

Let g = ((-bh - c)/a) so x = gz



5. Solve for z:
   
x²+ y²+ z²= r²
(gz)²+ (hz)²+ z²= r²
z²(g²+ h²+ 1) = r²
z²= r²/(g²+ h²+ 1)
z = +- sqrt(r²/(g²+ h²+ 1))



6. Find the intersects:
   

Let k = sqrt(r²/(g²+ h²+ 1))

The two great circles intersect at the points: (gk, hk, k) and (-gk, -hk, -k)

Convert those points back to spherical coordinates

lat_C= ArcSin(k/r)
lon_C= ArcTan(hk/gk)
or
lat_C= ArcSin(-k/r)
lon_C= ArcTan(-hk/-gk)

Great Circle intersects Latitude:   The equation for the plane that intersects the sphere at a given latitude is: z = z₀where z₀ = r*sin(lat₀).  Finding the points where a great circle intersects a parallel is just a variation of the method above.

Given two planes and a sphere:

(y₀z₁ - y₁z₀)x - (x₀z₁ - x₁z₀)y + (x₀y₁ - x₁y₀)z = 0

x²+ y²+ z²= r²

Solve for (x, y, z)

The algorithm used by the spheres package is as follows:

1. Define the constants in the planar equations. 

Let a = (y₀z₁ - y₁z₀)

Let b = -(x₀z₁ - x₁z₀)

Let c = (x₀y₁ - x₁y₀ )

So ax + by + cz = 0



2. Solve for x:

z = z₂ and ax + by + cz = 0 so
x = (-by - cz₂) / a

3. Solve for y:

x²+ y²+ z²= r²
((-by - cz₂) / a)² + y²+ z₂² = r²
((b²y²+2bycz₂ +c²z₂²) / a²) + y²+z₂² = r²
y²((b²+ a²) / a²) + y(2bcz₂/ a²) = r²- z₂² - ((c²z₂²) / a² )
y² + 2y(bcz₂/ a²) * (a²/(b²+ a²)) = (r²- z₂² - ((c²z₂²) / a² )) * (a²/(b²+ a²))

It's getting messy so
Let d = (bcz₂/ a²) * (a²/(b²+ a²)) = (bcz₂)/(b²+ a²)
Let rhs = (r²- z₂² - ((c²z₂²) / a² )) * (a²/(b²+ a²))
Let e = sqrt(rhs + d²)
so e² = rhs + d²
so e² - d² = rhs


then y²+ 2dy = e² - d²

Complete the square
y²+ 2dy + d² - d² = e² - d²
(y + d)²- d²= e²- d²
(y + d)²= e²- d²+ d²
(y + d) = +- e
y = -d +- e



4. Find the intersects:
   

The great circle and the parallel intersect at the points:

((bd - be -cz₂) / a, -d + e, z₂)  and   ((bd + be - cz₂) / a, -d - e, z₂)

Convert those back to spherical coordinates.

lat_C= ArcSin(z₂ / r)
lon_C= ArcTan((-d + e)/((bd - be -cz₂) / a))
or
lat_C= ArcSin(z₂ / r)
lon_C= ArcTan((-d - e)/((bd + be - cz₂) / a))

Great Circle intersects Small Circle:   A small circle is defined as any circle on a sphere that is not a great circle.  

This algorithm is not yet implemented by the spheres package.