Simpler closed form for spherical coordinates of points on great circle?

Question

Short version: how to simplify this expression given these conditions? Note that expr[[3]] is equal to 1, although Mathematica won't simplify it as such.

conds = {-Pi < tha, tha < Pi, -Pi < thb, thb < Pi, -Pi/2 < pha, pha <
Pi/2, -Pi/2 < phb, phb < Pi/2};

expr = {ArcTan[Cos[pha]*Cos[tha]*Cos[r*ArcCos[Cos[pha]*Cos[phb]*Cos[tha - thb] + 
         Sin[pha]*Sin[phb]]] + (((-Cos[pha])*Cos[tha]*Sin[pha]*Sin[phb] + 
       Cos[phb]*(Cos[thb]*Sin[pha]^2 + Cos[pha]^2*Sin[tha]*Sin[tha - thb]))*
      Sin[r*ArcCos[Cos[pha]*Cos[phb]*Cos[tha - thb] + Sin[pha]*Sin[phb]]])/
     Sqrt[Abs[Cos[pha]*((-Cos[phb])*Cos[tha - thb]*Sin[pha] + 
           Cos[pha]*Sin[phb])]^2 + (Cos[pha]*Cos[tha]*Sin[pha]*Sin[phb] - 
         Cos[phb]*(Cos[thb]*Sin[pha]^2 + Cos[pha]^2*Sin[tha]*Sin[tha - thb]))^
        2 + (Cos[pha]*Sin[pha]*Sin[phb]*Sin[tha] + Cos[pha]^2*Cos[phb]*
          Cos[tha]*Sin[tha - thb] - Cos[phb]*Sin[pha]^2*Sin[thb])^2], 
   Cos[pha]*Cos[r*ArcCos[Cos[pha]*Cos[phb]*Cos[tha - thb] + 
         Sin[pha]*Sin[phb]]]*Sin[tha] - 
    ((Cos[pha]*Sin[pha]*Sin[phb]*Sin[tha] + Cos[pha]^2*Cos[phb]*Cos[tha]*
        Sin[tha - thb] - Cos[phb]*Sin[pha]^2*Sin[thb])*
      Sin[r*ArcCos[Cos[pha]*Cos[phb]*Cos[tha - thb] + Sin[pha]*Sin[phb]]])/
     Sqrt[Abs[Cos[pha]*((-Cos[phb])*Cos[tha - thb]*Sin[pha] + 
           Cos[pha]*Sin[phb])]^2 + (Cos[pha]*Cos[tha]*Sin[pha]*Sin[phb] - 
         Cos[phb]*(Cos[thb]*Sin[pha]^2 + Cos[pha]^2*Sin[tha]*Sin[tha - thb]))^
        2 + (Cos[pha]*Sin[pha]*Sin[phb]*Sin[tha] + Cos[pha]^2*Cos[phb]*
          Cos[tha]*Sin[tha - thb] - Cos[phb]*Sin[pha]^2*Sin[thb])^2]], 
  ArcTan[
   Sqrt[(Cos[pha]*Cos[tha]*Cos[r*ArcCos[Cos[pha]*Cos[phb]*Cos[tha - thb] + 
            Sin[pha]*Sin[phb]]] + (((-Cos[pha])*Cos[tha]*Sin[pha]*Sin[phb] + 
          Cos[phb]*(Cos[thb]*Sin[pha]^2 + Cos[pha]^2*Sin[tha]*
             Sin[tha - thb]))*Sin[r*ArcCos[Cos[pha]*Cos[phb]*Cos[tha - thb] + 
             Sin[pha]*Sin[phb]]])/
        Sqrt[Abs[Cos[pha]*((-Cos[phb])*Cos[tha - thb]*Sin[pha] + 
              Cos[pha]*Sin[phb])]^2 + (Cos[pha]*Cos[tha]*Sin[pha]*Sin[phb] - 
            Cos[phb]*(Cos[thb]*Sin[pha]^2 + Cos[pha]^2*Sin[tha]*Sin[
                tha - thb]))^2 + (Cos[pha]*Sin[pha]*Sin[phb]*Sin[tha] + 
            Cos[pha]^2*Cos[phb]*Cos[tha]*Sin[tha - thb] - Cos[phb]*Sin[pha]^2*
             Sin[thb])^2])^2 + 
     (Cos[pha]*Cos[r*ArcCos[Cos[pha]*Cos[phb]*Cos[tha - thb] + 
            Sin[pha]*Sin[phb]]]*Sin[tha] - 
       ((Cos[pha]*Sin[pha]*Sin[phb]*Sin[tha] + Cos[pha]^2*Cos[phb]*Cos[tha]*
           Sin[tha - thb] - Cos[phb]*Sin[pha]^2*Sin[thb])*
         Sin[r*ArcCos[Cos[pha]*Cos[phb]*Cos[tha - thb] + Sin[pha]*Sin[phb]]])/
        Sqrt[Abs[Cos[pha]*((-Cos[phb])*Cos[tha - thb]*Sin[pha] + 
              Cos[pha]*Sin[phb])]^2 + (Cos[pha]*Cos[tha]*Sin[pha]*Sin[phb] - 
            Cos[phb]*(Cos[thb]*Sin[pha]^2 + Cos[pha]^2*Sin[tha]*Sin[
                tha - thb]))^2 + (Cos[pha]*Sin[pha]*Sin[phb]*Sin[tha] + 
            Cos[pha]^2*Cos[phb]*Cos[tha]*Sin[tha - thb] - Cos[phb]*Sin[pha]^2*
             Sin[thb])^2])^2], 
   Cos[r*ArcCos[Cos[pha]*Cos[phb]*Cos[tha - thb] + Sin[pha]*Sin[phb]]]*
     Sin[pha] + (Cos[pha]*((-Cos[phb])*Cos[tha - thb]*Sin[pha] + 
       Cos[pha]*Sin[phb])*Sin[r*ArcCos[Cos[pha]*Cos[phb]*Cos[tha - thb] + 
          Sin[pha]*Sin[phb]]])/
     Sqrt[Abs[Cos[pha]*((-Cos[phb])*Cos[tha - thb]*Sin[pha] + 
           Cos[pha]*Sin[phb])]^2 + (Cos[pha]*Cos[tha]*Sin[pha]*Sin[phb] - 
         Cos[phb]*(Cos[thb]*Sin[pha]^2 + Cos[pha]^2*Sin[tha]*Sin[tha - thb]))^
        2 + (Cos[pha]*Sin[pha]*Sin[phb]*Sin[tha] + Cos[pha]^2*Cos[phb]*
          Cos[tha]*Sin[tha - thb] - Cos[phb]*Sin[pha]^2*Sin[thb])^2]], 
  Sqrt[
   Abs[Cos[pha]*Cos[tha]*Cos[r*ArcCos[Cos[pha]*Cos[phb]*Cos[tha - thb] + 
            Sin[pha]*Sin[phb]]] + ((Cos[phb]*Cos[thb]*Sin[pha]^2 - 
          Cos[pha]*Cos[tha]*Sin[pha]*Sin[phb] + Cos[pha]^2*Cos[phb]*Sin[tha]*
           Sin[tha - thb])*Sin[r*ArcCos[Cos[pha]*Cos[phb]*Cos[tha - thb] + 
             Sin[pha]*Sin[phb]]])/
        Sqrt[Abs[Cos[pha]*((-Cos[phb])*Cos[tha - thb]*Sin[pha] + 
              Cos[pha]*Sin[phb])]^2 + (Cos[phb]*Cos[thb]*Sin[pha]^2 - 
            Cos[pha]*Cos[tha]*Sin[pha]*Sin[phb] + Cos[pha]^2*Cos[phb]*
             Sin[tha]*Sin[tha - thb])^2 + (Cos[pha]*Sin[pha]*Sin[phb]*
             Sin[tha] + Cos[pha]^2*Cos[phb]*Cos[tha]*Sin[tha - thb] - 
            Cos[phb]*Sin[pha]^2*Sin[thb])^2]]^2 + 
    Abs[Cos[pha]*Cos[r*ArcCos[Cos[pha]*Cos[phb]*Cos[tha - thb] + 
            Sin[pha]*Sin[phb]]]*Sin[tha] - 
       ((Cos[pha]*Sin[pha]*Sin[phb]*Sin[tha] + Cos[pha]^2*Cos[phb]*Cos[tha]*
           Sin[tha - thb] - Cos[phb]*Sin[pha]^2*Sin[thb])*
         Sin[r*ArcCos[Cos[pha]*Cos[phb]*Cos[tha - thb] + Sin[pha]*Sin[phb]]])/
        Sqrt[Abs[Cos[pha]*((-Cos[phb])*Cos[tha - thb]*Sin[pha] + 
              Cos[pha]*Sin[phb])]^2 + (Cos[phb]*Cos[thb]*Sin[pha]^2 - 
            Cos[pha]*Cos[tha]*Sin[pha]*Sin[phb] + Cos[pha]^2*Cos[phb]*
             Sin[tha]*Sin[tha - thb])^2 + (Cos[pha]*Sin[pha]*Sin[phb]*
             Sin[tha] + Cos[pha]^2*Cos[phb]*Cos[tha]*Sin[tha - thb] - 
            Cos[phb]*Sin[pha]^2*Sin[thb])^2]]^2 + 
    Abs[Cos[r*ArcCos[Cos[pha]*Cos[phb]*Cos[tha - thb] + Sin[pha]*Sin[phb]]]*
        Sin[pha] + (Cos[pha]*((-Cos[phb])*Cos[tha - thb]*Sin[pha] + 
          Cos[pha]*Sin[phb])*Sin[r*ArcCos[Cos[pha]*Cos[phb]*Cos[tha - thb] + 
             Sin[pha]*Sin[phb]]])/
        Sqrt[Abs[Cos[pha]*((-Cos[phb])*Cos[tha - thb]*Sin[pha] + 
              Cos[pha]*Sin[phb])]^2 + (Cos[pha]*Cos[tha]*Sin[pha]*Sin[phb] - 
            Cos[phb]*(Cos[thb]*Sin[pha]^2 + Cos[pha]^2*Sin[tha]*Sin[
                tha - thb]))^2 + (Cos[pha]*Sin[pha]*Sin[phb]*Sin[tha] + 
            Cos[pha]^2*Cos[phb]*Cos[tha]*Sin[tha - thb] - Cos[phb]*Sin[pha]^2*
             Sin[thb])^2]]^2]}

Longer version:

Given two points A and B on the unit sphere identified by longitude and latitude {tha, pha} and {thb, phb} (where th is the longitude), we can do the following:

(* using Wolfram Cloud *)

$Version

12.0.0 for Linux x86 (64-bit) (April 7, 2019)

(* using my own version of xyz2sph and sph2xyz (instead of CoordinateTransform) so spherical coordinates are in longitude, latitude, radius order *)

xyz2sph[x_,y_,z_] = {ArcTan[x,y], ArcTan[Sqrt[x^2+y^2],z], Norm[{x,y,z}]};
sph2xyz[th_,ph_,r_] = r*{Cos[th]*Cos[ph], Sin[th]*Cos[ph], Sin[ph]};

xyz2sph[l_] := Apply[xyz2sph,l];
sph2xyz[l_] := Apply[sph2xyz,l];

(* all variables are angles, longitudes limited to +-Pi and latitudues
limited to +-Pi/2 *)

conds = {-Pi < tha, tha < Pi, -Pi < thb, thb < Pi, -Pi/2 < pha, pha < Pi/2, -Pi/2 < phb, phb < Pi/2};

(* convert longitude, latitude to Cartesian *)

ptA = sph2xyz[tha, pha, 1];
ptB = sph2xyz[thb, phb, 1];

(* find the cross product, a (not necessarily unit) vector perpendicular to both *)

(* Simplify at each step to make things easier, although some of these
Simplify's will have no effect *)

perp = Simplify[Cross[ptA, ptB], conds];

(* find the vector perpendicular to perp and ptA, which lies in the
same plane as A and B; the order of the cross product is chosen so
that the parametrization I will do later goes from A to B *)

planar = Simplify[Cross[perp, ptA], conds];

(* normalize the planar vector *)

planarNorm= Simplify[planar/Norm[planar], conds];

(* distance-preserving parameterization of the great circle from A to B *)

point[t_] = Simplify[ptA*Cos[t] + planarNorm*Sin[t], conds];

(* the angular distance between A and B *)

angDist = Simplify[VectorAngle[ptA, ptB], conds];

(* the Cartesian vector "r" of the way between A and B on the great circle *)

rWayCartesian[r_] = Simplify[point[r*angDist], conds];

(* and the spherical equivalent *)

rWaySpherical[r_] = Simplify[xyz2sph[rWayCartesian[r]], conds];

(* which yields the nasty formula above *)

*)

Notes:

The above is also at https://www.wolframcloud.com/obj/a1b56ad7-4deb-4afd-aa53-436ec91b4c48
I suppressed most of the output above (except for $Version), because I wasn't sure if it would be useful.
I realize most of my terms above are functions of tha, pha, thb, phb, but I've written them as non-functions for simplicity.
I know that Equally spaced points on a great circle arc between 2 points asks a similar question, but I don't think a closed form is given explicitly.

@Somos I don't think so. I did see an answer on that as well, but I don't think my problem is related, and I think my problem is simpler. — , Jan 15 '20 at 03:43
To make sure I understand what you want. You have two not antipodal points on a sphere. The unique great circle they determine can be considered a new equator with zero longitude given by one of the points. Now you want to find the Cartesian coordinates of the point on this new equator at a given longitude. Is this what you want? — Somos, Jan 15 '20 at 23:26
@Somos Not exactly. I want to find the spherical coordinates of the great circle in the old coordinate system. In other words, if you fly from Albuquerque to Tokyo (for example), what function f(t) returns the spherical coordinates of Albuquerque at f(0), the spherical coordinates of Tokyo at f(1), and the spherical coordinates 'r'-way along the great circle at f(r). Note that f should preserve distance: eg, f(0.1) should be 1/10th of the distance along the great circle. Actually, I've already found such a formula, I just want the simplest version — , Jan 16 '20 at 15:58

Simpler closed form for spherical coordinates of points on great circle?

0 Answers0