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

Question

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))]}}]

Answer 1

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.

Answer 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[\ $.