# Simplest form for locus of latitudes/longitudes equidistant from two given latitudes/longitudes?

Given two latitudes/longitudes (th1,ph1 and th2,ph2), I want to find a simple formula for the locus of th3,ph3 that are equidistant from th1,ph1 and th2,ph2.

Mathematica happily spits out an answer (giving ph3 as a function of th3), but it's unbelivably long (below). Is there a simpler form? I realize I could convert to rectangular coordinates and back, but I'd ultimately just be calculating this result piecemeal.

Mathematica's "Simplify[]" doesn't help much, and "FullSimplify[]" hangs. Restricting the solution to Reals (ie, "Solve[eqn, Reals]") also hangs.

(* define distance using '=' (not ':=') for convenience *)

d2[th1_, ph1_, th2_, ph2_] =
(Sin[ph1]*Cos[th1] - Sin[ph2]*Cos[th2])^2 +
(Sin[ph1]*Sin[th1] - Sin[ph2]*Sin[th2])^2 +
(Sin[ph1] - Sin[ph2])^2;

(* and solve *)

s5 = Solve[{d2[th1, ph1, th3, ph3] == d2[th2, ph2, th3, ph3]}, {th3, ph3}];
InputForm[s5]

InputForm[
{
{th3 -> -ArcCos[(Csc[ph3]^2*(-(Sin[ph3]*(-4*Cos[th1]*Sin[ph1]^3 -
4*Cos[th1]^3*Sin[ph1]^3 + 4*Cos[th2]*Sin[ph1]^2*Sin[ph2] +
4*Cos[th1]^2*Cos[th2]*Sin[ph1]^2*Sin[ph2] + 4*Cos[th1]*Sin[ph1]*
Sin[ph2]^2 + 4*Cos[th1]*Cos[th2]^2*Sin[ph1]*Sin[ph2]^2 -
4*Cos[th2]*Sin[ph2]^3 - 4*Cos[th2]^3*Sin[ph2]^3 +
8*Cos[th1]*Sin[ph1]^2*Sin[ph3] - 8*Cos[th1]*Sin[ph1]*Sin[ph2]*
Sin[ph3] - 8*Cos[th2]*Sin[ph1]*Sin[ph2]*Sin[ph3] +
8*Cos[th2]*Sin[ph2]^2*Sin[ph3] - 4*Cos[th1]*Sin[ph1]^3*
Sin[th1]^2 + 4*Cos[th2]*Sin[ph1]^2*Sin[ph2]*Sin[th1]^2 +
4*Cos[th1]*Sin[ph1]*Sin[ph2]^2*Sin[th2]^2 - 4*Cos[th2]*Sin[ph2]^3*
Sin[th2]^2)) - Sqrt[Sin[ph3]^2*(-4*Cos[th1]*Sin[ph1]^3 -
4*Cos[th1]^3*Sin[ph1]^3 + 4*Cos[th2]*Sin[ph1]^2*Sin[ph2] +
4*Cos[th1]^2*Cos[th2]*Sin[ph1]^2*Sin[ph2] + 4*Cos[th1]*Sin[ph1]*
Sin[ph2]^2 + 4*Cos[th1]*Cos[th2]^2*Sin[ph1]*Sin[ph2]^2 -
4*Cos[th2]*Sin[ph2]^3 - 4*Cos[th2]^3*Sin[ph2]^3 +
8*Cos[th1]*Sin[ph1]^2*Sin[ph3] - 8*Cos[th1]*Sin[ph1]*Sin[ph2]*
Sin[ph3] - 8*Cos[th2]*Sin[ph1]*Sin[ph2]*Sin[ph3] +
8*Cos[th2]*Sin[ph2]^2*Sin[ph3] - 4*Cos[th1]*Sin[ph1]^3*Sin[th1]^
2 + 4*Cos[th2]*Sin[ph1]^2*Sin[ph2]*Sin[th1]^2 +
4*Cos[th1]*Sin[ph1]*Sin[ph2]^2*Sin[th2]^2 - 4*Cos[th2]*Sin[ph2]^
3*Sin[th2]^2)^2 - 4*Sin[ph3]^2*(4*Cos[th1]^2*Sin[ph1]^2 -
8*Cos[th1]*Cos[th2]*Sin[ph1]*Sin[ph2] + 4*Cos[th2]^2*
Sin[ph2]^2 + 4*Sin[ph1]^2*Sin[th1]^2 - 8*Sin[ph1]*Sin[ph2]*
Sin[th1]*Sin[th2] + 4*Sin[ph2]^2*Sin[th2]^2)*
(Sin[ph1]^4 + 2*Cos[th1]^2*Sin[ph1]^4 + Cos[th1]^4*Sin[ph1]^4 -
2*Sin[ph1]^2*Sin[ph2]^2 - 2*Cos[th1]^2*Sin[ph1]^2*Sin[ph2]^2 -
2*Cos[th2]^2*Sin[ph1]^2*Sin[ph2]^2 - 2*Cos[th1]^2*Cos[th2]^2*
Sin[ph1]^2*Sin[ph2]^2 + Sin[ph2]^4 + 2*Cos[th2]^2*Sin[ph2]^4 +
Cos[th2]^4*Sin[ph2]^4 - 4*Sin[ph1]^3*Sin[ph3] -
4*Cos[th1]^2*Sin[ph1]^3*Sin[ph3] + 4*Sin[ph1]^2*Sin[ph2]*
Sin[ph3] + 4*Cos[th1]^2*Sin[ph1]^2*Sin[ph2]*Sin[ph3] +
4*Sin[ph1]*Sin[ph2]^2*Sin[ph3] + 4*Cos[th2]^2*Sin[ph1]*
Sin[ph2]^2*Sin[ph3] - 4*Sin[ph2]^3*Sin[ph3] -
4*Cos[th2]^2*Sin[ph2]^3*Sin[ph3] + 4*Sin[ph1]^2*Sin[ph3]^2 -
8*Sin[ph1]*Sin[ph2]*Sin[ph3]^2 + 4*Sin[ph2]^2*Sin[ph3]^2 +
2*Sin[ph1]^4*Sin[th1]^2 + 2*Cos[th1]^2*Sin[ph1]^4*Sin[th1]^2 -
2*Sin[ph1]^2*Sin[ph2]^2*Sin[th1]^2 - 2*Cos[th2]^2*Sin[ph1]^2*
Sin[ph2]^2*Sin[th1]^2 - 4*Sin[ph1]^3*Sin[ph3]*Sin[th1]^2 +
4*Sin[ph1]^2*Sin[ph2]*Sin[ph3]*Sin[th1]^2 - 4*Sin[ph1]^2*
Sin[ph3]^2*Sin[th1]^2 + Sin[ph1]^4*Sin[th1]^4 +
8*Sin[ph1]*Sin[ph2]*Sin[ph3]^2*Sin[th1]*Sin[th2] -
2*Sin[ph1]^2*Sin[ph2]^2*Sin[th2]^2 - 2*Cos[th1]^2*Sin[ph1]^2*
Sin[ph2]^2*Sin[th2]^2 + 2*Sin[ph2]^4*Sin[th2]^2 +
2*Cos[th2]^2*Sin[ph2]^4*Sin[th2]^2 + 4*Sin[ph1]*Sin[ph2]^2*
Sin[ph3]*Sin[th2]^2 - 4*Sin[ph2]^3*Sin[ph3]*Sin[th2]^2 -
4*Sin[ph2]^2*Sin[ph3]^2*Sin[th2]^2 - 2*Sin[ph1]^2*Sin[ph2]^2*
Sin[th1]^2*Sin[th2]^2 + Sin[ph2]^4*Sin[th2]^4)]))/
(2*(4*Cos[th1]^2*Sin[ph1]^2 - 8*Cos[th1]*Cos[th2]*Sin[ph1]*Sin[ph2] +
4*Cos[th2]^2*Sin[ph2]^2 + 4*Sin[ph1]^2*Sin[th1]^2 -
8*Sin[ph1]*Sin[ph2]*Sin[th1]*Sin[th2] + 4*Sin[ph2]^2*Sin[th2]^2))]},
{th3 -> ArcCos[(Csc[ph3]^2*(-(Sin[ph3]*(-4*Cos[th1]*Sin[ph1]^3 -
4*Cos[th1]^3*Sin[ph1]^3 + 4*Cos[th2]*Sin[ph1]^2*Sin[ph2] +
4*Cos[th1]^2*Cos[th2]*Sin[ph1]^2*Sin[ph2] + 4*Cos[th1]*Sin[ph1]*
Sin[ph2]^2 + 4*Cos[th1]*Cos[th2]^2*Sin[ph1]*Sin[ph2]^2 -
4*Cos[th2]*Sin[ph2]^3 - 4*Cos[th2]^3*Sin[ph2]^3 +
8*Cos[th1]*Sin[ph1]^2*Sin[ph3] - 8*Cos[th1]*Sin[ph1]*Sin[ph2]*
Sin[ph3] - 8*Cos[th2]*Sin[ph1]*Sin[ph2]*Sin[ph3] +
8*Cos[th2]*Sin[ph2]^2*Sin[ph3] - 4*Cos[th1]*Sin[ph1]^3*
Sin[th1]^2 + 4*Cos[th2]*Sin[ph1]^2*Sin[ph2]*Sin[th1]^2 +
4*Cos[th1]*Sin[ph1]*Sin[ph2]^2*Sin[th2]^2 - 4*Cos[th2]*Sin[ph2]^3*
Sin[th2]^2)) - Sqrt[Sin[ph3]^2*(-4*Cos[th1]*Sin[ph1]^3 -
4*Cos[th1]^3*Sin[ph1]^3 + 4*Cos[th2]*Sin[ph1]^2*Sin[ph2] +
4*Cos[th1]^2*Cos[th2]*Sin[ph1]^2*Sin[ph2] + 4*Cos[th1]*Sin[ph1]*
Sin[ph2]^2 + 4*Cos[th1]*Cos[th2]^2*Sin[ph1]*Sin[ph2]^2 -
4*Cos[th2]*Sin[ph2]^3 - 4*Cos[th2]^3*Sin[ph2]^3 +
8*Cos[th1]*Sin[ph1]^2*Sin[ph3] - 8*Cos[th1]*Sin[ph1]*Sin[ph2]*
Sin[ph3] - 8*Cos[th2]*Sin[ph1]*Sin[ph2]*Sin[ph3] +
8*Cos[th2]*Sin[ph2]^2*Sin[ph3] - 4*Cos[th1]*Sin[ph1]^3*
Sin[th1]^2 + 4*Cos[th2]*Sin[ph1]^2*Sin[ph2]*Sin[th1]^2 +
4*Cos[th1]*Sin[ph1]*Sin[ph2]^2*Sin[th2]^2 - 4*Cos[th2]*
Sin[ph2]^3*Sin[th2]^2)^2 - 4*Sin[ph3]^2*
(4*Cos[th1]^2*Sin[ph1]^2 - 8*Cos[th1]*Cos[th2]*Sin[ph1]*Sin[ph2] +
4*Cos[th2]^2*Sin[ph2]^2 + 4*Sin[ph1]^2*Sin[th1]^2 -
8*Sin[ph1]*Sin[ph2]*Sin[th1]*Sin[th2] + 4*Sin[ph2]^2*Sin[th2]^2)*
(Sin[ph1]^4 + 2*Cos[th1]^2*Sin[ph1]^4 + Cos[th1]^4*Sin[ph1]^4 -
2*Sin[ph1]^2*Sin[ph2]^2 - 2*Cos[th1]^2*Sin[ph1]^2*Sin[ph2]^2 -
2*Cos[th2]^2*Sin[ph1]^2*Sin[ph2]^2 - 2*Cos[th1]^2*Cos[th2]^2*
Sin[ph1]^2*Sin[ph2]^2 + Sin[ph2]^4 + 2*Cos[th2]^2*Sin[ph2]^4 +
Cos[th2]^4*Sin[ph2]^4 - 4*Sin[ph1]^3*Sin[ph3] -
4*Cos[th1]^2*Sin[ph1]^3*Sin[ph3] + 4*Sin[ph1]^2*Sin[ph2]*
Sin[ph3] + 4*Cos[th1]^2*Sin[ph1]^2*Sin[ph2]*Sin[ph3] +
4*Sin[ph1]*Sin[ph2]^2*Sin[ph3] + 4*Cos[th2]^2*Sin[ph1]*Sin[ph2]^2*
Sin[ph3] - 4*Sin[ph2]^3*Sin[ph3] - 4*Cos[th2]^2*Sin[ph2]^3*
Sin[ph3] + 4*Sin[ph1]^2*Sin[ph3]^2 - 8*Sin[ph1]*Sin[ph2]*
Sin[ph3]^2 + 4*Sin[ph2]^2*Sin[ph3]^2 + 2*Sin[ph1]^4*Sin[th1]^2 +
2*Cos[th1]^2*Sin[ph1]^4*Sin[th1]^2 - 2*Sin[ph1]^2*Sin[ph2]^2*
Sin[th1]^2 - 2*Cos[th2]^2*Sin[ph1]^2*Sin[ph2]^2*Sin[th1]^2 -
4*Sin[ph1]^3*Sin[ph3]*Sin[th1]^2 + 4*Sin[ph1]^2*Sin[ph2]*Sin[ph3]*
Sin[th1]^2 - 4*Sin[ph1]^2*Sin[ph3]^2*Sin[th1]^2 +
Sin[ph1]^4*Sin[th1]^4 + 8*Sin[ph1]*Sin[ph2]*Sin[ph3]^2*Sin[th1]*
Sin[th2] - 2*Sin[ph1]^2*Sin[ph2]^2*Sin[th2]^2 -
2*Cos[th1]^2*Sin[ph1]^2*Sin[ph2]^2*Sin[th2]^2 +
2*Sin[ph2]^4*Sin[th2]^2 + 2*Cos[th2]^2*Sin[ph2]^4*Sin[th2]^2 +
4*Sin[ph1]*Sin[ph2]^2*Sin[ph3]*Sin[th2]^2 - 4*Sin[ph2]^3*Sin[ph3]*
Sin[th2]^2 - 4*Sin[ph2]^2*Sin[ph3]^2*Sin[th2]^2 -
2*Sin[ph1]^2*Sin[ph2]^2*Sin[th1]^2*Sin[th2]^2 +
Sin[ph2]^4*Sin[th2]^4)]))/(2*(4*Cos[th1]^2*Sin[ph1]^2 -
8*Cos[th1]*Cos[th2]*Sin[ph1]*Sin[ph2] + 4*Cos[th2]^2*Sin[ph2]^2 +
4*Sin[ph1]^2*Sin[th1]^2 - 8*Sin[ph1]*Sin[ph2]*Sin[th1]*Sin[th2] +
4*Sin[ph2]^2*Sin[th2]^2))]},
{th3 -> -ArcCos[(Csc[ph3]^2*(-(Sin[ph3]*(-4*Cos[th1]*Sin[ph1]^3 -
4*Cos[th1]^3*Sin[ph1]^3 + 4*Cos[th2]*Sin[ph1]^2*Sin[ph2] +
4*Cos[th1]^2*Cos[th2]*Sin[ph1]^2*Sin[ph2] + 4*Cos[th1]*Sin[ph1]*
Sin[ph2]^2 + 4*Cos[th1]*Cos[th2]^2*Sin[ph1]*Sin[ph2]^2 -
4*Cos[th2]*Sin[ph2]^3 - 4*Cos[th2]^3*Sin[ph2]^3 +
8*Cos[th1]*Sin[ph1]^2*Sin[ph3] - 8*Cos[th1]*Sin[ph1]*Sin[ph2]*
Sin[ph3] - 8*Cos[th2]*Sin[ph1]*Sin[ph2]*Sin[ph3] +
8*Cos[th2]*Sin[ph2]^2*Sin[ph3] - 4*Cos[th1]*Sin[ph1]^3*
Sin[th1]^2 + 4*Cos[th2]*Sin[ph1]^2*Sin[ph2]*Sin[th1]^2 +
4*Cos[th1]*Sin[ph1]*Sin[ph2]^2*Sin[th2]^2 - 4*Cos[th2]*Sin[ph2]^3*
Sin[th2]^2)) + Sqrt[Sin[ph3]^2*(-4*Cos[th1]*Sin[ph1]^3 -
4*Cos[th1]^3*Sin[ph1]^3 + 4*Cos[th2]*Sin[ph1]^2*Sin[ph2] +
4*Cos[th1]^2*Cos[th2]*Sin[ph1]^2*Sin[ph2] + 4*Cos[th1]*Sin[ph1]*
Sin[ph2]^2 + 4*Cos[th1]*Cos[th2]^2*Sin[ph1]*Sin[ph2]^2 -
4*Cos[th2]*Sin[ph2]^3 - 4*Cos[th2]^3*Sin[ph2]^3 +
8*Cos[th1]*Sin[ph1]^2*Sin[ph3] - 8*Cos[th1]*Sin[ph1]*Sin[ph2]*
Sin[ph3] - 8*Cos[th2]*Sin[ph1]*Sin[ph2]*Sin[ph3] +
8*Cos[th2]*Sin[ph2]^2*Sin[ph3] - 4*Cos[th1]*Sin[ph1]^3*Sin[th1]^
2 + 4*Cos[th2]*Sin[ph1]^2*Sin[ph2]*Sin[th1]^2 +
4*Cos[th1]*Sin[ph1]*Sin[ph2]^2*Sin[th2]^2 - 4*Cos[th2]*Sin[ph2]^
3*Sin[th2]^2)^2 - 4*Sin[ph3]^2*(4*Cos[th1]^2*Sin[ph1]^2 -
8*Cos[th1]*Cos[th2]*Sin[ph1]*Sin[ph2] + 4*Cos[th2]^2*
Sin[ph2]^2 + 4*Sin[ph1]^2*Sin[th1]^2 - 8*Sin[ph1]*Sin[ph2]*
Sin[th1]*Sin[th2] + 4*Sin[ph2]^2*Sin[th2]^2)*
(Sin[ph1]^4 + 2*Cos[th1]^2*Sin[ph1]^4 + Cos[th1]^4*Sin[ph1]^4 -
2*Sin[ph1]^2*Sin[ph2]^2 - 2*Cos[th1]^2*Sin[ph1]^2*Sin[ph2]^2 -
2*Cos[th2]^2*Sin[ph1]^2*Sin[ph2]^2 - 2*Cos[th1]^2*Cos[th2]^2*
Sin[ph1]^2*Sin[ph2]^2 + Sin[ph2]^4 + 2*Cos[th2]^2*Sin[ph2]^4 +
Cos[th2]^4*Sin[ph2]^4 - 4*Sin[ph1]^3*Sin[ph3] -
4*Cos[th1]^2*Sin[ph1]^3*Sin[ph3] + 4*Sin[ph1]^2*Sin[ph2]*
Sin[ph3] + 4*Cos[th1]^2*Sin[ph1]^2*Sin[ph2]*Sin[ph3] +
4*Sin[ph1]*Sin[ph2]^2*Sin[ph3] + 4*Cos[th2]^2*Sin[ph1]*
Sin[ph2]^2*Sin[ph3] - 4*Sin[ph2]^3*Sin[ph3] -
4*Cos[th2]^2*Sin[ph2]^3*Sin[ph3] + 4*Sin[ph1]^2*Sin[ph3]^2 -
8*Sin[ph1]*Sin[ph2]*Sin[ph3]^2 + 4*Sin[ph2]^2*Sin[ph3]^2 +
2*Sin[ph1]^4*Sin[th1]^2 + 2*Cos[th1]^2*Sin[ph1]^4*Sin[th1]^2 -
2*Sin[ph1]^2*Sin[ph2]^2*Sin[th1]^2 - 2*Cos[th2]^2*Sin[ph1]^2*
Sin[ph2]^2*Sin[th1]^2 - 4*Sin[ph1]^3*Sin[ph3]*Sin[th1]^2 +
4*Sin[ph1]^2*Sin[ph2]*Sin[ph3]*Sin[th1]^2 - 4*Sin[ph1]^2*
Sin[ph3]^2*Sin[th1]^2 + Sin[ph1]^4*Sin[th1]^4 +
8*Sin[ph1]*Sin[ph2]*Sin[ph3]^2*Sin[th1]*Sin[th2] -
2*Sin[ph1]^2*Sin[ph2]^2*Sin[th2]^2 - 2*Cos[th1]^2*Sin[ph1]^2*
Sin[ph2]^2*Sin[th2]^2 + 2*Sin[ph2]^4*Sin[th2]^2 +
2*Cos[th2]^2*Sin[ph2]^4*Sin[th2]^2 + 4*Sin[ph1]*Sin[ph2]^2*
Sin[ph3]*Sin[th2]^2 - 4*Sin[ph2]^3*Sin[ph3]*Sin[th2]^2 -
4*Sin[ph2]^2*Sin[ph3]^2*Sin[th2]^2 - 2*Sin[ph1]^2*Sin[ph2]^2*
Sin[th1]^2*Sin[th2]^2 + Sin[ph2]^4*Sin[th2]^4)]))/
(2*(4*Cos[th1]^2*Sin[ph1]^2 - 8*Cos[th1]*Cos[th2]*Sin[ph1]*Sin[ph2] +
4*Cos[th2]^2*Sin[ph2]^2 + 4*Sin[ph1]^2*Sin[th1]^2 -
8*Sin[ph1]*Sin[ph2]*Sin[th1]*Sin[th2] + 4*Sin[ph2]^2*Sin[th2]^2))]},
{th3 -> ArcCos[(Csc[ph3]^2*(-(Sin[ph3]*(-4*Cos[th1]*Sin[ph1]^3 -
4*Cos[th1]^3*Sin[ph1]^3 + 4*Cos[th2]*Sin[ph1]^2*Sin[ph2] +
4*Cos[th1]^2*Cos[th2]*Sin[ph1]^2*Sin[ph2] + 4*Cos[th1]*Sin[ph1]*
Sin[ph2]^2 + 4*Cos[th1]*Cos[th2]^2*Sin[ph1]*Sin[ph2]^2 -
4*Cos[th2]*Sin[ph2]^3 - 4*Cos[th2]^3*Sin[ph2]^3 +
8*Cos[th1]*Sin[ph1]^2*Sin[ph3] - 8*Cos[th1]*Sin[ph1]*Sin[ph2]*
Sin[ph3] - 8*Cos[th2]*Sin[ph1]*Sin[ph2]*Sin[ph3] +
8*Cos[th2]*Sin[ph2]^2*Sin[ph3] - 4*Cos[th1]*Sin[ph1]^3*
Sin[th1]^2 + 4*Cos[th2]*Sin[ph1]^2*Sin[ph2]*Sin[th1]^2 +
4*Cos[th1]*Sin[ph1]*Sin[ph2]^2*Sin[th2]^2 - 4*Cos[th2]*Sin[ph2]^3*
Sin[th2]^2)) + Sqrt[Sin[ph3]^2*(-4*Cos[th1]*Sin[ph1]^3 -
4*Cos[th1]^3*Sin[ph1]^3 + 4*Cos[th2]*Sin[ph1]^2*Sin[ph2] +
4*Cos[th1]^2*Cos[th2]*Sin[ph1]^2*Sin[ph2] + 4*Cos[th1]*Sin[ph1]*
Sin[ph2]^2 + 4*Cos[th1]*Cos[th2]^2*Sin[ph1]*Sin[ph2]^2 -
4*Cos[th2]*Sin[ph2]^3 - 4*Cos[th2]^3*Sin[ph2]^3 +
8*Cos[th1]*Sin[ph1]^2*Sin[ph3] - 8*Cos[th1]*Sin[ph1]*Sin[ph2]*
Sin[ph3] - 8*Cos[th2]*Sin[ph1]*Sin[ph2]*Sin[ph3] +
8*Cos[th2]*Sin[ph2]^2*Sin[ph3] - 4*Cos[th1]*Sin[ph1]^3*
Sin[th1]^2 + 4*Cos[th2]*Sin[ph1]^2*Sin[ph2]*Sin[th1]^2 +
4*Cos[th1]*Sin[ph1]*Sin[ph2]^2*Sin[th2]^2 - 4*Cos[th2]*
Sin[ph2]^3*Sin[th2]^2)^2 - 4*Sin[ph3]^2*
(4*Cos[th1]^2*Sin[ph1]^2 - 8*Cos[th1]*Cos[th2]*Sin[ph1]*Sin[ph2] +
4*Cos[th2]^2*Sin[ph2]^2 + 4*Sin[ph1]^2*Sin[th1]^2 -
8*Sin[ph1]*Sin[ph2]*Sin[th1]*Sin[th2] + 4*Sin[ph2]^2*Sin[th2]^2)*
(Sin[ph1]^4 + 2*Cos[th1]^2*Sin[ph1]^4 + Cos[th1]^4*Sin[ph1]^4 -
2*Sin[ph1]^2*Sin[ph2]^2 - 2*Cos[th1]^2*Sin[ph1]^2*Sin[ph2]^2 -
2*Cos[th2]^2*Sin[ph1]^2*Sin[ph2]^2 - 2*Cos[th1]^2*Cos[th2]^2*
Sin[ph1]^2*Sin[ph2]^2 + Sin[ph2]^4 + 2*Cos[th2]^2*Sin[ph2]^4 +
Cos[th2]^4*Sin[ph2]^4 - 4*Sin[ph1]^3*Sin[ph3] -
4*Cos[th1]^2*Sin[ph1]^3*Sin[ph3] + 4*Sin[ph1]^2*Sin[ph2]*
Sin[ph3] + 4*Cos[th1]^2*Sin[ph1]^2*Sin[ph2]*Sin[ph3] +
4*Sin[ph1]*Sin[ph2]^2*Sin[ph3] + 4*Cos[th2]^2*Sin[ph1]*Sin[ph2]^2*
Sin[ph3] - 4*Sin[ph2]^3*Sin[ph3] - 4*Cos[th2]^2*Sin[ph2]^3*
Sin[ph3] + 4*Sin[ph1]^2*Sin[ph3]^2 - 8*Sin[ph1]*Sin[ph2]*
Sin[ph3]^2 + 4*Sin[ph2]^2*Sin[ph3]^2 + 2*Sin[ph1]^4*Sin[th1]^2 +
2*Cos[th1]^2*Sin[ph1]^4*Sin[th1]^2 - 2*Sin[ph1]^2*Sin[ph2]^2*
Sin[th1]^2 - 2*Cos[th2]^2*Sin[ph1]^2*Sin[ph2]^2*Sin[th1]^2 -
4*Sin[ph1]^3*Sin[ph3]*Sin[th1]^2 + 4*Sin[ph1]^2*Sin[ph2]*Sin[ph3]*
Sin[th1]^2 - 4*Sin[ph1]^2*Sin[ph3]^2*Sin[th1]^2 +
Sin[ph1]^4*Sin[th1]^4 + 8*Sin[ph1]*Sin[ph2]*Sin[ph3]^2*Sin[th1]*
Sin[th2] - 2*Sin[ph1]^2*Sin[ph2]^2*Sin[th2]^2 -
2*Cos[th1]^2*Sin[ph1]^2*Sin[ph2]^2*Sin[th2]^2 +
2*Sin[ph2]^4*Sin[th2]^2 + 2*Cos[th2]^2*Sin[ph2]^4*Sin[th2]^2 +
4*Sin[ph1]*Sin[ph2]^2*Sin[ph3]*Sin[th2]^2 - 4*Sin[ph2]^3*Sin[ph3]*
Sin[th2]^2 - 4*Sin[ph2]^2*Sin[ph3]^2*Sin[th2]^2 -
2*Sin[ph1]^2*Sin[ph2]^2*Sin[th1]^2*Sin[th2]^2 +
Sin[ph2]^4*Sin[th2]^4)]))/(2*(4*Cos[th1]^2*Sin[ph1]^2 -
8*Cos[th1]*Cos[th2]*Sin[ph1]*Sin[ph2] + 4*Cos[th2]^2*Sin[ph2]^2 +
4*Sin[ph1]^2*Sin[th1]^2 - 8*Sin[ph1]*Sin[ph2]*Sin[th1]*Sin[th2] +
4*Sin[ph2]^2*Sin[th2]^2))]}}]
-

(In the following the equator has latitude $\theta=0$.)

You are given two points $$x_i\ :=\ (\cos\theta_i\cos\phi_i,\cos\theta_i\sin\phi_i,\sin\theta_i)\ \in S^2\qquad(i=1,2)\ .$$ If $\theta_1=\theta_2$ then the two poles of $S^2$ both belong to the locus $\gamma$ in question. In this case $\gamma$ is the meridian circle $\phi={\rm const.}\$ where $$\phi\in {\phi_1+\phi_2+2\pi{\mathbb Z} \over2}\ .$$ Assume now that $\theta_1\ne\theta_2$. In this case $\gamma$ is a great circle going around the $x_3$-axis; therefore it has a parametric representation of the form $$\gamma:\quad \phi\ \mapsto\ x(\phi)=\bigl(\cos\theta(\phi)\cos\phi,\ \cos\theta(\phi)\sin\phi,\ \sin\theta(\phi)\bigr)\qquad (0\leq\phi\leq2\pi)$$ where the function $\phi\mapsto\theta(\phi)$ remains to be determined. To this end we finally use the geometric characterization of $\gamma$: All points $x(\phi)$ must satisfy the equation $$x(\phi)\cdot x_1\ =\ x(\phi)\cdot x_2\ ,$$ where $\cdot$ denotes the scalar product. Doing the calculation one finds that $\tan\theta(\phi)$ is given by $$\tan\theta(\phi)={\cos(\phi-\phi_2)\cos\theta_2-\cos(\phi-\phi_1)\cos\theta_1 \over \sin\theta_1 -\sin\theta_2}\ ,$$ and taking the $\arctan$ on both sides one gets an explicit expression for $\theta(\phi)\in \bigl]-{\pi\over2}, {\pi\over2}\bigr[\$.

-

The problem is a lot easier in $\vec{x} = (x,y,z)$ form, (i.e. vector notation), so convert your two points to $\vec{x}_1$ and $\vec{x}_2$. Now the point midway between these is simply half of $\vec{x}_1+\vec{x}_2$. Normalizing this vector gives a point on the sphere which is equidistant: $\vec{m}=(\vec{x}_1+\vec{x}_2)/|\vec{x}_1+\vec{x}_2|$. Note that this fails when the two points are on opposite sides of the sphere, as is expected.

To get the remainder of the equidistant points you need to rotate this point around the sphere on an axis that is parallel to the two points. That is, the axis you want to rotate around is $\vec{x}_2-\vec{x}_1$. So normalize this vector (which fails when the two points are coincident, as expected), to get $\vec{u} = (\vec{x}_2-\vec{x}_1) /| \vec{x}_2-\vec{x}_1|$.

The rotation needed is defined by the $SO(3)$ rotation $$\vec{u}'(\theta) = \exp (\theta\begin{pmatrix}0&u_3&-u_2\\-u_3&0&u_1\\u_2&-u_1&0\end{pmatrix})\;\vec{m}$$ for $\theta$ from $0$ to $2\pi$.

To get this to work you need to have a way of converting back and forth to Cartesian form but that should be easy to figure out.

-
OK, but I'd like a formula that goes from polar to polar. Ultimately, following your procedure will do this, but is the result any simpler than the one I gave above? – barrycarter Apr 6 '11 at 15:53
I doubt it is simpler, but I like a formula that is understandable, and I wouldn't be surprised if it uses fewer multiplications. Also, I believe that the exponential function above can be solved exactly by looking at the powers of the matrix. – Carl Brannen Apr 6 '11 at 19:40
I might've mis-stated my question. I agree that a step-by-step process is better. In fact, that's how I got my hideous answer. Maybe I should tag this question "math-golf". I'm just curious how short an explicit polar-to-polar answer can be. I'm trying to convert this to Perl (without necessarily using a matrix multiplying library). – barrycarter Apr 7 '11 at 1:35