Auxiliary Latitude Formulas

Finding the coefficients numerically and symbolically

Craig Rollins

National Geospatial-Intelligence Agency
Office of GEOINT Sciences
St. Louis, MO

Ellipse

Off[General::spell1, General::spell];
a = 5; b = 4; c = (Sqrt[a^2 - b^2]);
ϵ = c/a;
deg = Pi/180;
θ1 = 50  deg // N;
x1 = a Cos[θ1];
z1 = b Sin[θ1];
ϕ1 = ArcTan[a/b*Tan[θ1]];
ψ1 = ArcTan[b/a*Tan[θ1]];
ω1 = (π/2)*
          EllipticPi[(ϵ^2), ϕ1, (ϵ^2)]/EllipticPi[(ϵ^2), Pi/2, (ϵ^2)];
R2 = 
        a (Sqrt[ (1/2) + ((1 - (ϵ^2))/(4  ϵ)) Log[((1 + ϵ)/(1 - ϵ))]]);
β1 = 
        ArcSin[((((b/R2)))^2) (Sin[ϕ1]/(2 (1 - (ϵ^2) Sin[ϕ1]^2)) + (1/(4  ϵ)) 
                  Log[((1 + ϵ Sin[ϕ1])/(1 - ϵ Sin[ϕ1]))])];
χ1 = (π/2) - 
          2 ArcTan[
              Tan[(π/4) - ϕ(1/2)] (((((1 + ϵ Sin[ϕ1])/(1 - ϵ Sin[ϕ1]))))^(ϵ/2))];
<< Graphics`ImplicitPlot`;
gra1 = 
        ImplicitPlot[
          x^2/a^2 + z^2/b^2 == 1, {x, -a - 1, 
            a + 1}, {z, -b - 1, b + 1}, AspectRatio -> Automatic, 
          AxesOrigin -> {0, 0}, PlotStyle -> Thickness[0.01], 
          Ticks -> None, 
          TextStyle -> {FontFamily -> "Arial", FontSize -> 20}, 
          AxesLabel -> {"x", "z"},  
          DisplayFunction -> Identity];
txta = 
        Graphics[
          Text[StyleForm["a", FontWeight -> "Plain", 
              FontColor -> Red,  FontSize -> 20], {2.0, 
              0.4}]];
txtb = 
        Graphics[
          Text[
            StyleForm["b", FontWeight -> "Plain", 
              FontColor -> Red,  FontSize -> 20], {0.4, 
              2.0}]];
Show[gra1, txta, txtb, 
        DisplayFunction -> $DisplayFunction];
On[General::spell1, General::spell];

ϵ = (Sqrt[1 - ((((b/a)))^2)])

(Geodetic) Latitude ϕ

gra2 = 
        ImplicitPlot[
          z - z1 == Tan[ϕ1]*(x - x1), {x, -a, a}, {z, 0, 
            b}, PlotStyle -> Red, AspectRatio -> Automatic, 
          DisplayFunction -> Identity];
gra3 = 
        Graphics[
          Text[StyleForm["ϕ", FontWeight -> "Bold", 
              FontColor -> Red,  FontSize -> 20], {2.4, 
              0.5}]];
Show[gra1, gra2, gra3, 
        DisplayFunction -> $DisplayFunction];

Geocentric Latitude ψ

gra4 = 
        ImplicitPlot[
          z - z1 == Tan[ψ1]*(x - x1), {x, -a, a}, {z, 0, 
            b}, PlotStyle -> Red, AspectRatio -> Automatic, 
          DisplayFunction -> Identity];
gra5 = 
        Graphics[
          Text[StyleForm["ψ", FontWeight -> "Bold", 
              FontColor -> Red,  FontSize -> 20], {1.4, 
              0.5}]];
Show[gra1, gra4, gra5, 
        DisplayFunction -> $DisplayFunction];

Parametric Latitude θ

gra6 = 
        ListPlot[{{a Cos[θ1], b Sin[θ1]}}, 
          PlotStyle -> {PointSize[0.03], Red}, 
          DisplayFunction -> Identity];
gra7 = 
        Graphics[
          Text[StyleForm["x = a Cos(θ) \nz = b Sin(θ)", 
              FontWeight -> "Bold", FontColor -> Red,  
              FontSize -> 20,  FontSlant -> "Italic"], {6.5, 
              3.0}]];
Show[gra1, gra6,  gra7, DisplayFunction -> $DisplayFunction, 
        PlotRange -> All];

Rectifying Latitude ω

gra8 = 
        ParametricPlot[{{a Cos[θ], b Sin[θ]} // 
              Evaluate}, {θ, 0, θ1}, 
          AspectRatio -> Automatic, Axes -> False, 
          PlotStyle -> {{Thickness[0.012], Red}}, 
          DisplayFunction -> Identity];
gra9 = 
        Graphics[
          Text[StyleForm["ω", FontWeight -> "Bold", 
              FontColor -> Red,  FontSize -> 20], {5.0, 
              2.5}]];
Show[gra1, gra8, gra9,  
        DisplayFunction -> $DisplayFunction];

1 deg =(1/90) Arclength Equator-to-Pole

Authalic Latitude β
Conformal Latitude χ

gr10 = 
        ParametricPlot3D[{a Cos[θ] Cos[λ], 
            a Cos[θ] Sin[λ], 
            b Sin[θ]}, {θ, -Pi/2, 
            Pi/2}, {λ, -Pi, Pi}, AspectRatio -> Automatic, 
          Axes -> False, Boxed -> False, 
          PlotPoints -> {18, 36}, 
          DisplayFunction -> Identity];
gr11 = 
        ParametricPlot3D[{a Cos[θ] Cos[λ], 
            a Cos[θ] Sin[λ], b Sin[θ]}, {θ, 
            0, θ1}, {λ, -Pi, Pi}, 
          AspectRatio -> Automatic, Axes -> False, 
          Boxed -> False, PlotPoints -> {30, 45}, 
          DisplayFunction -> Identity];
gr12 = 
        Show[gr10, gr11, 
          DisplayFunction -> Identity];
gr13 = 
        ParametricPlot3D[{R2 Cos[β] Cos[λ], 
            R2 Cos[β] Sin[λ], 
            R2 Sin[β]}, {β, -Pi/2, 
            Pi/2}, {λ, -Pi, Pi}, AspectRatio -> Automatic, 
          Axes -> False, Boxed -> False, 
          PlotPoints -> {18, 36}, 
          DisplayFunction -> Identity];
gr14 = 
        ParametricPlot3D[{R2 Cos[β] Cos[λ], 
            R2 Cos[β] Sin[λ], R2 Sin[β]}, {β, 
            0, β1}, {λ, -Pi, Pi}, 
          AspectRatio -> Automatic, Axes -> False, 
          Boxed -> False, PlotPoints -> {30, 45}, 
          DisplayFunction -> Identity];
gr15 = 
        Show[gr13, gr14, 
          DisplayFunction -> Identity];
Show[GraphicsArray[{{gr12, gr15}}, GraphicsSpacing -> 0], 
        DisplayFunction -> $DisplayFunction];

Comparisons (ϵ = 3/5)

({{ϕ, Geodetic, ϕ1/deg}, {θ, Parametric, θ1/deg}, {β, Authalic, β1/deg}, {ω, Rectifying, ω1/deg}, {χ, Conformal, χ1/deg}, {ψ, Geocentric, ψ1/deg}})//TableForm//(NumberForm[#, {8, 6}, 
            NumberPadding -> {"", "0"}] &)

Commonalities

General series method

ω = χ + (Subscript[a, 2]) Sin(
                  2  χ) + (Subscript[a, 4]) Sin(4  χ) + (Subscript[a, 6]) Sin(
                  6  χ) + (Subscript[a, 8]) Sin(8  χ) +  ...

Oscar Adams (1921)

From χ to ω

The coefficients are given by:

(Subscript[a, k]) = (2/(k m)) (∫_( 0)^( π/2))Cos[k χ]/(w^3) \[DifferentialD]ϕ

where m is a constant, and w and χ are functions of ϕ given by:

m = (∫_( 0)^( π/2))(1/(w^3)) \[DifferentialD]ϕ
w = (Sqrt[1 - (ϵ^2) Sin[ϕ]^2])
χ = (π/2) - 
        2  ArcTan[
            Tan[(π/4) - (ϕ/2)] (((((1 + ϵ Sin[ϕ])/(1 - ϵ Sin[ϕ]))))^(ϵ/2))]

Numerical (Subscript[a, k])

Oscillating integrand is problematic

ϵ = 0.6; k = 12;
Print["ϵ = ", ϵ // N,  "   k = ", 
        k];
Plot[
χ = (π/2) - 
            2  ArcTan[
                Tan[(π/4) - (ϕ/2)] (((((1 + ϵ Sin[ϕ])/(1 - ϵ Sin[ϕ]))))^(ϵ/2))]; 
w = (Sqrt[1 - (ϵ^2) Sin[ϕ]^2]); 
Cos[k χ]/(w^3), {ϕ, 0, π/2}];

ϵ = 0.081; k = 12;
Print["ϵ = ", ϵ // N,  "   k = ", 
        k];
Plot[
χ = (π/2) - 
            2  ArcTan[
                Tan[(π/4) - (ϕ/2)] (((((1 + ϵ Sin[ϕ])/(1 - ϵ Sin[ϕ]))))^(ϵ/2))]; 
w = (Sqrt[1 - (ϵ^2) Sin[ϕ]^2]); 
Cos[k χ]/(w^3), {ϕ, 0, π/2}];

Numerical (Subscript[a, k])

Use Mathematica's extended arithmetic

$MinPrecision = 20;
invfla  =  300.8017000`20;  (* Everest ellipsoid (s) *)
fla = 1/invfla;
ϵsq = 
        fla (2 - fla);
ϵ = 
        Sqrt[ϵsq];
k = 
        16;
m = 
        NIntegrate[1/(w^3), {ϕ, 0, π/2}, WorkingPrecision -> 50, 
          AccuracyGoal -> 30];
 acoeff=(2/(k m))NIntegrate[
\(χ = π\/2 - 
                    2  ArcTan[
                        Tan[π\/4 - ϕ\/2] \((\(1 + ϵ\ Sin[ϕ]\)\/\(1 - ϵ\ Sin[ϕ]\))\)\^\(ϵ/2\)]; 
w = \@\(1 - \(ϵ\^2\) Sin[ϕ]\^2\); 
Cos[k\ χ]\/w\^3\),{ϕ, 0, π/2},
(* \(EDIT\ FOLLOWING\ TO\ CHANGE\ INTERNAL\ PRECISION\) *)
WorkingPrecision -> 50,AccuracyGoal -> 30] ;
Print[
        "\[Epsilon] = ", ϵ // N,  "   k = ", k,  
        "   coefficient = ", acoeff];

Numerical (Subscript[a, k])

Internal Precision not high enough causes these problems:

Symbolic (Subscript[a, k])

Clear[ϵ, ϵsq];  (* Clear numerical work *)

sim[x_] := 
      Simplify[x, 
        Assumptions -> {ϵ > 0, ϵ < 
              1, ϕ ≥ 
              0, ϕ ≤ π/2}]

k = 8;     (* Identify coefficient (Subscript[a, k]) *)
high = 12;    (* High power of ϵ  *)

m=Series[
        EllipticPi[(ϵ^2), (π/2), (ϵ^2)], {ϵ, 0, 
          high}]

(* Following is the function to be integrated *)
\(\(χ = π\/2 - 
            2  ArcTan[
                Tan[π\/4 - ϕ\/2] \((\(1 + ϵ\ Sin[ϕ]\)\/\(1 - ϵ\ Sin[ϕ]\))\)\^\(ϵ/2\)];\)

      \(w = \@\(1 - \(ϵ\^2\) Sin[ϕ]\^2\);\)

      \(integrand = Cos[k\ χ]\/w\^3;\)\)

Symbolic (Subscript[a, k])

(* Expand the integrand as a series in powers of ϵ *)
ser1 = 
        Series[integrand, {ϵ, 0, high}] // sim

(* Convert powers of Sin and Cos to multiple - 
            angle Sin and Cos *)
ser2 = Map[TrigReduce, ser1]

(* Integrate term by term *)
ser3 = 
        Integrate[ser2, {ϕ, 0, π/2}]

(* Apply constant factor -- non - trivial *)
ser4 = (2/(k m)) ser3

The results should be:
((49561 (ϵ^8))/41287680)+((28081 (ϵ^10))/20643840)+((27404441 (ϵ^12))/29727129600)+O[ϵ]^13