:Mathematica Version: 4.0

:Copyright: public domain

:Author: John Hendrickson (debguy@sourceforge.net)

:Summary:

1) list of  geocentric <> geodetic <> XYZ  related equations*

2) quick plots using them (3d poly globe)

3) wikipedia article summarizing the above

4) one equation is by author of this paper.  shows that using (new) geometry observation, Heiskanen and Moritz or Ferrari equations were not necessary as the problem is not inheirently recursive as they suggested. (see article)

*gleened from many onlne sources

geodesy equations

key

BELOW:

λ	longitude, swept from x to y on x-y plane
φ	angle Rn makes with x (phi)
Φ	angle r makes with x (Phi)
ϕ	angle Rn makes with x (phi)
φλh lat lon geodetic, height above x,y
ΦΛr lat lon spherical, radius

ϕ φ Φ	there are three \phi i've seen, not always used the same way
		they don't cut paste play across apps or web well
			mm TEX unicode ASCII mathml html not same as "new web"
			req's hand edit per site / per app

l	λ	longitude, swept from x to y on x-y plane
L	Λ	longitude, swept from x to y on x-y plane
		same as λ unless ellipsoid has mult. e's
w	φ	latitude,  angle Rn makes with x
W	Φ	latitude,  angle r makes with x
h	h	height above sphere or ellipse from earth reference point
		oddly books use same symbol however earth is defined
		i assume h is part of r. for large data use (r+h)
wlh		XYZg	lat lon geodetic, height above x,y
WLr		XYZe	lat lon spherical, radius

See section: "same equations, with random notes about arriving at for extra details" and "Wikipedia article".  there are many papers online and on wikipedia about geodesy equations, they do not all agree on naming and form (some appear incorrect).

general interest: try free software "stellarium".  there are navigation softwares free.  wikipedia has meny geodesy articles. mathematica 9 has built in solar simulation. http://www.fourmilab.ch/yoursky/ is a good example of java interactive astronomy.

basic data and equations  definitions w/o notes

Clear[approxPolarDist, approxPolarDist2, GeodeticToGeocentric, 
      GeocentricToGeodetic, GeocentricToGeodetic2, WLr, wlh, wlhNewt, XYZe, 
      XYZg, secondEccentricityFun, firstEccentricitySquaredFun]

mag[
          r_] := (Plus @@ (r^2))^(1/
              2);
unit[r_] := 
        Module[{mr}, mr = mag[r]; 
          If[mr === 0, 0, r/mr]];
midpt[r_] := 
        Module[{mr}, lr = Length@r; 
          If[lr === 0, 0, Plus @@ (r)/lr]];

EarthEquatorialRadius  =  6378136.6 Meter;
EarthPolarRadius  =  6356751.9 Meter;

XYZ data is ECEF

(*WGS - 
        84 datum*) 
 (*a*) 
semiMajorAxis = 
        6378137.0 Meter;
 (*1/
        f*) 
reciprocalOfFlattening = 
        298.257223;
  
semiMinorAxis = 
        6356752  
          Meter;
  
firstEccentricitySquared = 
        6.69437999014 10^-3;
  
secondEccentricitySquared = 6.73949674228 10^-3;

note e is same by elliptic e eq. found in math books

e = firstEccentricitySquared^(1/2);
a = semiMajorAxis;
b = semiMinorAxis;

h is height above surface

flattening = 1/298257223563;

b[a_,f_]:=a(1-f);

firstEccentricitySquaredFun[a_, b_] := 
        1 - (b^2/a^2);

wgse = 
        firstEccentricitySquared^(1/2);
fefe = firstEccentricitySquaredFun[a, b]^(1/2);

secondEccentricityFun[a_, b_] := (a^2/b^2) - 
          1;

rn[a_, e_, w_] := a/(1 - e^2 Sin[w]^2)^(1/2);

rn[a_, b_, w_, "b"] := 
        a^2/(a^2 (Cos[w]^2)  + (b^2 Sin[w]^2))^(1/
                2);

XYZg[w_, l_, h_, a_, e_] := 
        Module[{Rn}, 
          Rn = rn[a, e, w]; 
{(Rn + h) Cos[w] 
              Cos[l], 
(Rn + h) Cos[w] 
              Sin[l], 
((1 - e^2) Rn + h) 
              Sin[w]}];

instructions: start with w == geocentric

wlh[X_, Y_, Z_, a_, e_, w_] := 
        Module[{Rn, rc}, 
rc = mag[{X, Y}]; 
          Rn = rn[a, e, 
              w]; 
{ArcTan[(Z + e^2 Rn Sin[w])/
                rc], 
Sign[Y] 
              ArcCos[X/rc], 
(rc/Cos[w]) - 
              Rn}];

wlhNewt[X_, Y_, Z_, a_, e_, 
        count__:  4] := 
Module[{r, rc, rc2, w, W, L, h, i, 
          W2}, 
r = mag[{X, Y, Z}]; 
        rc = mag[{X, Y}]; 
w = ArcTan[rc/Z];   
Do[
h = 
            rc/Cos[w] - rn[a, e, w]; 
w = 
            ArcTan[(Z/
                    rc) (1 - 
                      e^2 (rn[a, e, 
                              w]/(rn[a, e, w] + 
                                h)))^-1], 
{i, 1, 
            count}]; 
h = 
          rc/Cos[w] - rn[a, e, w]; 
{W, L, r} = 
          WLr[X, Y, Z]; 
rc2 = (rn[a, e, w] + h) 
            Cos[W]; 
Zn = (rn[a, e, w] + h) Sin[
              w]; 
RnZn = 
          e^2 rn[a, e, w];   
{w, 
          L, h}]

wlh[X_, Y_, Z_, a_, e_, w_] := 
        Module[{Rn, rc}, 
rc = mag[{X, Y}]; 
          Rn = rn[a, e, 
              w]; 
{ArcTan[(Z + e^2 Rn Sin[w])/
                rc], 
Sign[Y] 
              ArcCos[X/rc], 
(rc/Cos[w]) - 
              Rn}];

XYZe[W_, L_, 
          r_] := 
{r Cos[W] Cos[
              L], 
r Cos[W] Sin[
              L], 
r Sin[W]};

XYZe[W_, L_, r_, 
          h_] := 
{(r + h) Cos[W] Cos[
              L], 
(r + h) Cos[W] Sin[
              L], 
(r + h) Sin[W]};

WLr[X_, Y_, Z_] := 
        Module[{r}, 
r = 
            mag[{X, Y, Z}]; 
{If[r === 0, 0, 
              ArcSin[Z/r]], 
If[Y === 0, 
              If[X >= 0, 0, 
                Pi], 
Sign[
                  Y] (Pi - Sign[Y] 2 ArcTan[(X/Y)])/
                  2], 
r}];

WLr[X_, Y_, Z_, 
          r_] := 
{If[r === 0, 0, 
            ArcSin[Z/r]], 
If[Y === 0, 
            If[X >= 0, 0, 
              Pi], 
Sign[
                Y] (Pi - Sign[Y] 2 ArcTan[(X/Y)])/
                2], 
r};

GeodeticToGeocentric[w_, e_, rn_, h_] := 
        If[Abs[w] === Pi/2, w, 
          ArcTan[(1 - e^2 (rn/(rn + h))) Tan[w]]];

GeocentricToGeodetic[W_, e_, rn_, h_] := 
        If[Abs[W] === Pi/2, W, 
          ArcTan[(1 - e^2 (rn/(rn + h)))^-1 
              Tan[W]]];

GeocentricToGeodetic2[{X_, Y_, Z_}, a_, b_] := 
        Module[{gradFP}, 
gradFP = {2  X/a^2, 2  Y/a^2, 
              2  Z/b^2};   
Pi/2 - 
            ArcCos[unit[gradFP] . {0, 0, 1}]];

approxPolarDist[a_, e_, ul_, ll_] := 
      Module[{w}, 
        a Integrate[(1 - e^2  Sin[w]^2)^(1/2), {w, ll, 
              ul}]]

approxPolarDist2[a_, e_, ul_, ll_] := 
      Module[{w}, (ul - 
              ll) ((a^(3/2) + b^(3/2))/2)^(2/
                3)]

Clear[wlhEDm3, wlhEDm1, wlhm1, wlhEDm2, ellipseDilation]

ellipseDilation[{X_, Y_, Z_}, a_, 
          b_] := 
((b^2 X^2 + b^2 Y^2 + 
                  a^2 Z^2)^(1/2))/(a b);

wlhEDm1[X_, Y_, Z_, a_, b_, 
        e_] := 
Module[{P, We, Le, re, w, l, h, 
          ed}, 
P = {X, Y, Z}; 
{We, Le, 
            re} = WLr[X, Y, Z]; 
w = 
          GeocentricToGeodetic[We, e, 1, 0]; 
ed = 
          ellipseDilation[P, a, b]; 
h = 
          If[ed === 0, 0, mag[P] - mag[P/ed]]; 
{w, Le, 
          h}]

wlhEDm2[X_, Y_, Z_, a_, b_, 
        e_] := 
Module[{P, We, Le, re, w, l, h, 
          ed}, 
P = {X, Y, Z}; 
ed = 
          ellipseDilation[P, a, b]; 
{We, Le, re} = 
          WLr[X/ed, Y/ed, Z/ed]; 
w = 
          GeocentricToGeodetic[We, e, 1, 0]; 
h = 
          If[ed === 0, 0, mag[P] - mag[P/ed]]; 
{w, Le, 
          h}]

wlhEDm3[X_, Y_, Z_, a_, b_, 
        e_] := 
Module[{P, Le, Pe, w, l, h, gradFP, ed, 
          mxy}, 
P = {X, Y, Z}; 
        mxy = mag[{X, Y}]; 
ed = 
          ellipseDilation[P, a, b]; 
w = 
          GeocentricToGeodetic2[P, a, b]; 
Le = 
          If[mxy === 0, 0, Sign[Y] ArcCos[X/mxy]]; 
h = 
          If[ed === 0, 0, mag[P] - mag[P/ed]]; 
{w, Le, 
          h}]

make world (ww) using eqs ^^ and plot it

lat lon 0,0  is ang 0

details

Needs just WorldData`, if mm 9 does not have it you can find WorldData.m on Wolfram Library.

Needs["Miscellaneous`WorldData`"]

make ww

the below makes ww, rough earth polys to plot , pretty bad but will suffice

Clear[tessfunPoly4, areamax, junk1, junk2, junk3, wfg, wfg2, w, w0, ww, 
      www, x, y, z, w, v, i, l, splitLen, tessDia, hull]

Module[{i, l, v}, 
        i =  .3; 
w0 = 
          Map[(
v = WorldData[#]; 
                i += 1/300; 
l = 
                  Length@Flatten@v; 
If[l == 4, 
                  v[[1]] = Append[v[[1]], v[[1, 1]] + {10, 25}]; 
                  v[[1]] = 
                    Append[v[[1]], v[[1, 1]]]]; 
If[
                  l == 2, 
                  v[[1]] = Append[v[[1]], v[[1, 1]] + {10, 50}]; 
                  v[[1]] = Append[v[[1]], v[[1, 1]] - {50, 10}]; 
                  v[[1]] = 
                    Append[v[[1]], v[[1, 1]]]]; 
If[
                  Length@Flatten@v === 
                    2 Length@v, 
v = Polygon[v[[1]]], 
                  v = (Map[Polygon, 
                        v])]; 
 (* {GrayLevel[i], 
                    v} *) 
v
) &
, World]
];

splitLen = 100;
tessDia = 500;

splitLen = 100;
tessDia = 2000;
hull = False;

splitLen = 
        500;
tessDia = 1000;
hull = False;

\(tessfunPoly4[vl_List] := 
      Module[{inc, v, center, len, x, vv, xx, 
          m}, 

If[
          hull =!= True, 
v = 
            vl
, 
v = 
            ConvexHull[vl]; 
v = 
            Map[\((vl[\([#]\)])\) &, v]; 
v = 
            Append[v, 
              First@v];
]; 

 (*\ fill\ in\ incomplete\ *) 
If[
          Last[v] =!= First[v], 
          AppendTo[v, First[v]]]; 
len = 
          Length[v]; 
 (*\ 
          find\ center\ make\ pie\ slices\ ; \ 
          is\ cc\ if\ poly\ is\ *) 
center = 
          Apply[Plus, v]/len; 
 (*\ 
          break\ 2\ longest\ sides\ *) 
If[len < 4, 
          Return[v]]; 

Do[
m = 
            mag[v[\([inc - 1]\)] - v[\([inc]\)]]; 
If[
            m > splitLen, 
 (*\ \(Print["m=", 
                  m];\)\ *) 
v = 
              Insert[v, \((v[\([inc - 1]\)] + v[\([inc]\)])\)/2, 
                inc]
], {inc, 2, 
            len}
]; 
If[once, 
          Return[v]]; 

xx = \(Table[
                mag[v[\([inc]\)] - center] // N
, {inc, 1, 
                  len}] // Sort\) // Reverse; 
xx = 
          xx[\([1]\)]; 
If[xx < tessDia, 
          Return[v]]; 
len = 
          Length[v]; 
Table[{v[\([inc - 1]\)], 
            v[\([inc]\)], center}
, {inc, 2, 
            len}]
]\)

junk3[v_List]:=If[areamax[v][[1]]>500,Print["hit"];True,False]

junk3[v_List]:=True

junk2[v_List]:=If[junk3[v],junk1[tessfunPoly4[v]],Polygon[v]]

junk2 =. ;

junk2[v_List] := 
      Module[{x, y}, 
x = 
          tessfunPoly4[v // N]; 
If[
          NumericQ[x[[1, 1]]] === True, 
Polygon[
            x // N], 
Map[Polygon, 
            x // N]
]
]

junk2[v_List]:=If[junk3[v],Polygon[tessfunPoly4[v]],Polygon[v]]

ww =. ; v =. ;

w0 /. Polygon[v_] -> v // Flatten // Length

ww = w0 /. Polygon[v_] -> junk2[v];

ww /. Polygon[v_] -> v // Flatten // Length

plot world

Clear[we, ws, wg, x, y, z];

Map[Polygon, {{{1, 2}, {3, 4}}, {{1, 2}, {3, 4}}}]

ang = 0;

wfg2 = 
        WireFrame[
          Graphics3D[
Sphere[(a/Meter)  .95 // N, 20, 
              15]]];

world = ww;

triangle that marks (0,0) which shows 0,0 in data is
	+x is towards us (it's either that or -x is)
	and y is L/R
earth-centered earth-fixed (ECEF or ECF) on wiki says this is good
	because w does this and wg is same, assume same
TODO is still to check w and wg for acc. and prec.

world=Prepend[world,{Polygon[10{{10,0},{-10,-10},{-10,10},{10,0}}]}];

---------------------------
plot geocentric data  (must be fit where it intersects ellipse)

double converting will work but poor example, extra steps
	XYZg, φλh, GeodeticToGeocentric

ellipseDilation[{X_,Y_,Z_},a_,b_]:=
    ((b^2 X^2+b^2 Y^2+a^2 Z^2)^(1/2))/(a b);

Pe[P_] := P/ellipseDilation[P, a/Meter, b/Meter] // N;
we = 
        world /. {x_?NumericQ, y_?NumericQ} -> 
Pe[
              XYZe[x (Pi/10800), y (Pi/10800), a/Meter]];
we // Graphics3D // Show;

show use P/ellipse not P-ellipse
	Show[Graphics3D[{Line[{{0,0,0},{0,2,4}}],Line[{{0,1,0},1/5{0,4,8}}]}]]

there are two ways to get normal. however vector will be easier.  vector to what ??

---------------------------
plot geodetic data

wg = 
        world /. {x_?NumericQ, y_?NumericQ} -> 
XYZg[
                x (Pi/10800), y (Pi/10800), 0, a/Meter, e] // 
          N;
wg // Graphics3D // Show;

---------------------------
plot spheric data w/o height

ws = 
        world /. {x_?NumericQ, y_?NumericQ} -> 
XYZe[
                x (Pi/10800), y (Pi/10800), a/Meter] // 
          N;
ws // Graphics3D // Show;

--------------------------------------------------------------
quick-plots
	show lat lon 0,0  is ang 0

ang = -(Pi/2 + Pi/2 + Pi/2);
ang = 0;

www = 
      Show[
Graphics3D[
{GrayLevel[ .7], 
            we}, 
ViewPoint -> {2 Cos[ang], 
              2  Sin[ang], 0}, 
          ViewVertical -> {0, 0, 1}
, 
          LightSources -> {
{{-4, 1, 4}, 
                RGBColor[ .7,  .7,  .1]}, 
{{4, -1, 4}, 
                RGBColor[ .5,  .5,  .9]}}
]
, wfg2, Boxed -> False
]

www = 
      Show[
Graphics3D[
{GrayLevel[ .7], 
            wg}, 
ViewPoint -> {2 Cos[ang], 
              2  Sin[ang], 0}, 
          ViewVertical -> {0, 0, 1}
, 
          LightSources -> {
{{-4, 1, 4}, 
                RGBColor[ .7,  .7,  .1]}, 
{{4, -1, 4}, 
                RGBColor[ .5,  .5,  .9]}}
]
, wfg2, Boxed -> False
]

--------------------------------------------------------------
make .gif for wiki, rotated and lit diff.

ang = -(Pi/2 + Pi/2 + Pi/2);

g = 
      Show[
Graphics3D[
{GrayLevel[ .7], 
            wg}
, 
          ViewPoint -> {2 Cos[ang], 2  Sin[ang], 0}, 
          ViewVertical -> {0, 0, 1}
, 
          LightSources -> {
{{-4, 1, 4}, 
                RGBColor[ .7,  .7,  .1]}, 
{{4, -1, 4}, 
                RGBColor[ .5,  .5,  .9]}}
]
, wfg2, Boxed -> False
]

Display["geodeticquickplot.gif",g,"GIF",ImageSize->{240,240}];

rayshade`
(is set not evaluatable)

note: section set not evaluatable

rayshade` (with POVray !) is a package that can be found via Wolfram Library

Needs["rayshade`"]

Options[Rayshade]

NotebookOpen["~/rayshade/rayshade-help.nb"]

\(\(\(SetDirectory["/tmp/"]\)\(\ \)\( (*\ "/home/user/"\ *) \)\)\)

gr = 
        ContourPlot3D[
          x y z, {x, -1, 1}, {y, -1, 1}, {z, -1, 1}, 
          Contours -> { .1}];

Rayshade["file.ray", gr, RayPreView -> True, 
      RayFastLight -> True]

ang = 0;
ang = -(Pi/2 + Pi/2 + Pi/2);

g = 
      Show[
Graphics3D[
{GrayLevel[ .7], 
            wg}
, 
          ViewPoint -> {2 Cos[ang], 2  Sin[ang], 0}, 
          ViewVertical -> {0, 0, 1}
, 
          LightSources -> {
{{-4, 1, 4}, 
                RGBColor[ .7,  .7,  .1]}, 
{{4, -1, 4}, 
                RGBColor[ .5,  .5,  .9]}}
]
, wfg2, Boxed -> False
]

Rayshade["file.ray", g, RayFastLight -> True, 
      Boxed -> False, RayTubeLines -> True]

NOTE: need to tesselate BEFORE ^^ mapping, RayTesselate->True

NOTE: just learned pov may be mroe problematic to render at distances
			maybe for earth dist there's some magic to find don't care ...
	  luckly using mathmatica no adjustment is impossible ...

wfg3 = 
        WireFrame[
          Graphics3D[
Sphere[ .95 // N, 20, 
              15]]];

wp = 
        wg /. {x_?NumericQ, y_?NumericQ, 
              z_?NumericQ} -> 
{x/(a/Meter), 
              y/(a/Meter), z/(b/Meter)};

gp = 
      Show[
Graphics3D[
{GrayLevel[ .7], 
            wp}
, 
          ViewPoint -> {2 Cos[ang], 2  Sin[ang], 0}, 
          ViewVertical -> {0, 0, 1}
, 
          LightSources -> {
{{-4, 1, 4}, 
                RGBColor[ .7,  .7,  .1]}, 
{{4, -1, 4}, 
                RGBColor[ .5,  .5,  .9]}}
]
, wfg3, Boxed -> False
]

Rayshade["file2.pov", gp, RayPOV -> True, Boxed -> False, 
      RayTesselate -> True, RayTubeLines -> True]

was for no sphere unit of 1

Rayshade["file2.pov",{Graphics3D[{Sphere[{0,0,0},.95]}],g},RayPOV->True,RayBoxLines->False,RayTesselate->True]

due to tesselation and (?order of? ?need tracer?) looks like crap

Rayshade["file1.pov",www,RayTesselate->True,RayPreView->True,RayFastLight->True,RayPOV->True,RayGlass->True]

clear plot data

Clear[ws, we, wg, wp]

Clear[world, ww, g]

Wikipedia article
(text, for the above equations and plots)

==Quick Plots: geodetic, spheric, geocentric==

==Quick Plots: geodetic, spheric, geocentric==
Before going through tons of equations note:  (1) a capt would use methods for naut. almanac (not "value added" eq's) and ship software get updates (3) a nice map takes immense data, layers, and planning, and many avail online free.  NOTE the following is not verified:

  ArcMinuteToRadian[am_]:=am Pi/10800;
  e=firstEccentricitySquared^(1/2);
  a=semiMajorAxis;b=semiMinorAxis;

convert geodetic (lat,lon) in world to (x,y,z) and plot (most common usage)

  (* world={Polygon[{{lat,lon},...}], Polygon ...};*)
  rn[a_,e_,w_]:=a/(1-e^2 Sin[w]^2)^(1/2);
  XYZg[w_,l_,h_,a_,e_]:={
  (rn[a,e,w]+h)Cos[w]Cos[l],
  (rn[a,e,w]+h)Cos[w]Sin[l],
  ((1-e^2)rn[a,e,w]+h)Sin[w] }
  wg=w/.{x_?NumericQ,y_?NumericQ};
  ->XYZg[x(Pi/10800),y(Pi/10800),0,a/Meter,e]//N;
  wg//Graphics3D//Show;

[[File:geodeticquickplot.gif]]

convert spheric (lat,lon,r) in world to (x,y,z) and plot

  (* world={Polygon[{{lat,lon,r},...}], Polygon ...};*)
  XYZe[W_,L_,r_]:=
  {r Cos[W] Cos[L],r Cos[W] Sin[L],r Sin[W]};
  ws=w/.{x_?NumericQ,y_?NumericQ}->
  XYZe[x (Pi/10800),y (Pi/10800),r]//N;
  ws//Graphics3D//Show;

convert geocentric (lat,lon) in world to (x,y,z), and plot

  (* world={Polygon[{{lat,lon},...}], Polygon ...};*)
  ellipseDilation[{X_,Y_,Z_},a_,b_]:=
  ((b^2 X^2+b^2 Y^2+a^2 Z^2)^(1/2))/(a b);
  Pe[P_]:=P/ellipseDilation[P,a/Meter,b/Meter]//N;
  we=world/.{x_?NumericQ,y_?NumericQ}->
  Pe[XYZe[x(Pi/10800),y(Pi/10800),a/Meter]];
  we//Graphics3D//Show;

note: geocentric method: assume XYZe is on dilated ellipsoid (same e as earth), find dilation and contract to earth ellipse.

note: scale of 1:world is so "large" that geocentric v. geodtic can't be seen

== XYZ to Geodetic + all quad. eq. ==

== XYZ to Geodetic + all quad. eq. ==

Purpose: show pratical cut and paste equations that use map (lat,lon) throughout and work in all octants, for XYZ to geodetic.  note not fully tested !

  mag[r_]:=(Plus@@(r^2))^(1/2);

  unit[r_]:=Module[{mr},mr=mag[r];If[mr===0,0,r/mr]];

  GeodeticToGeocentric[w_,e_,rn_,h_]:=
  If[Abs[w]===Pi/2,w,
  ArcTan[(1-e^2(rn/(rn+h))) Tan[w]]];

  GeocentricToGeodetic[W_,e_,rn_,h_]:=
  If[Abs[W]===Pi/2,W,
  ArcTan[(1-e^2(rn/(rn+h)))^-1 Tan[W]]];

  XYZg[w_,l_,h_,a_,e_]:=Module[{Rn},Rn=rn[a,e,w];
  {(Rn+h) Cos[w] Cos[l],
  (Rn+h) Cos[w] Sin[l],
  ((1-e^2) Rn+h) Sin[w]}];

  XYZe[W_,L_,r_]:=
  {r Cos[W] Cos[L],
  r Cos[W] Sin[L],
  r Sin[W]};

  WLr[X_,Y_,Z_]:=Module[{r},r=mag[{X,Y,Z}];
  {If[r===0,0,ArcSin[Z/r]],
  If[Y===0,If[X>=0,0,Pi],
  Sign[Y](Pi-Sign[Y] 2 ArcTan[(X/Y)])/2],
  r}];

  rn[a_,e_,w_]:=a/(1-e^2 Sin[w]^2)^(1/2);

source:
  Heiskanen and Moritz (1967) revised
topic:
  XYZ to geodetic angle for WGS-84 earth ellipse
  (referred to on wiki pages as more difficult)
note:
  many papers sugg. to find wlh loop is req.
  (w and h cannot both be known by one calc.)
  WGS-84 inexact e: h .2 meter err at least
  i beleive obsoleted by arbitrary precision.
  (~ femto meter diff between methods, by arb)
  "converges in count (4 at most) to < 5 cm."
  "eq. for h diverges at poles." (fixup omited)
key:
  XYZ is ECEF, +X to you, +Y R, +Z U,(0,0) GMT
  h is geod. height above ellipse, if any
  rn line from pt to Z-axis usu. through y-axis
  rn rad. of curvature, line "normal" to surf
  Zn is height from rn origin to Z

  wlhNewt[X_,Y_,Z_,a_,e_,count__:4]:=
  Module[{r,rc,rc2,w,W,L,h,i,W2},
  r=mag[{X,Y,Z}];rc=mag[{X,Y}];
  w=ArcTan[rc/Z];(* geocentric *)
  Do[ h=rc/Cos[w]-rn[a,e,w];
  w=ArcTan[(Z/rc)
  (1-e^2(rn[a,e,w]/(rn[a,e,w]+h)))^-1],
  {i,1,count}];
  h=rc/Cos[w]-rn[a,e,w];
  {W,L,r}=WLr[X,Y,Z];
  rc2=(rn[a,e,w]+h) Cos[W];
  Zn=(rn[a,e,w]+h) Sin[w];
  RnZn=e^2 rn[a,e,w];(*Z+RnZn==Zn*)
  {w,L,h}];

  wlh[X_,Y_,Z_,a_,e_,w_]:=Module[{Rn,rc},
  rc=mag[{X,Y}];Rn=rn[a,e,w];
  {ArcTan[(Z+e^2 Rn Sin[w])/rc],
  Sign[Y] ArcCos[X/rc],
  (rc/Cos[w])-Rn}];

submitter: the below sugg. rn and h can be both known w/o newt. approx. because contour ellipses share properties. apparently same result as Heiskanen and Moritz (1967) using CAM

  ellipseDilation[{X_,Y_,Z_},a_,b_]:=
  ((b^2 X^2+b^2 Y^2+a^2 Z^2)^(1/2))/(a b);

  GeocentricToGeodetic2[{X_,Y_,Z_},a_,b_]:=
  Module[{gradFP},(*surf N vec at P*)
  gradFP={2X/a^2,2Y/a^2,2Z/b^2};
  Pi/2-ArcCos[unit[gradFP].{0,0,1}]];

  wlhEDm1[X_,Y_,Z_,a_,b_,e_]:=
  Module[{P,We,Le,re,w,l,h,ed},
  P={X,Y,Z};
  {We,Le,re}=WLr[X,Y,Z];
  w=GeocentricToGeodetic[We,e,1,0];
  ed=ellipseDilation[P,a,b];
  h=If[ed===0,-mag[P],mag[P]-mag[P/ed]];
  {w,Le,h}];

  wlhEDm3[X_,Y_,Z_,a_,b_,e_]:=
  Module[{P,Le,Pe,w,l,h,ed,mxy},
  P={X,Y,Z};mxy=mag[{X,Y}];
  ed=ellipseDilation[P,a,b];
  w=GeocentricToGeodetic2[P,a,b];
  Le=If[mxy===0,0,Sign[Y] ArcCos[X/mxy]];
  h=If[ed===0,-mag[P],mag[P]-mag[P/ed]];
  {w,Le,h}];

same equations, with random notes about arriving at
(is set non-evaluatable)

Need

limited notebook, does not require these

\(Map[
      Needs, {"Miscellaneous`Units`", "Geometry`Rotations`", "Miscellaneous`Geodesy`", "Miscellaneous`PhysicalConstants`", "Miscellaneous`WorldData`", "Miscellaneous`WorldNames`", "Miscellaneous`CityData`", "Miscellaneous`WorldPlot`", "Geometry`Rotations`", "Calculus`VectorAnalysis`"}]\)

this is wise if re-graphic allot , prevents kernel freeze

Unprotect[In, Out];

Clear[In, Out];

Protect[In, Out];

$HistoryLength = 5;

MemoryInUse[]

eq

Clear[approxPolarDist, approxPolarDist2, GeodeticToGeocentric, 
      GeocentricToGeodetic, GeocentricToGeodetic2, WLr, wlh, wlhNewt, XYZe, 
      XYZg, secondEccentricityFun, firstEccentricitySquaredFun]

mag[r_]:=(Plus@@(r^2))^(1/2);
unit[r_]:=Module[{mr},mr=mag[r];If[mr===0,0,r/mr]];
midpt[r_]:=Module[{mr},lr=Length@r;If[lr===0,0,Plus@@(r)/lr]];

from spheric			from ellipsoidal
(r+h) cos(ϕ) sin(λ)		r->Rn
(r+h) cos(ϕ) sin(λ)		r->Rn
(r+h) sin(ϕ)			r->(1-e^2) Rn , on z

EarthEquatorialRadius  =  6378136.6 Meter;
EarthPolarRadius  =  6356751.9 Meter;

XYZ data is ECEF

note it matters allot if one uses less accurate
	WGS-84 datum, however in some practices that is
	assumed and correction and nav tables apply after it

(*WGS - 
        84 datum*) 
 (*a*) 
semiMajorAxis = 
        6378137.0 Meter;
 (*1/
        f*) 
reciprocalOfFlattening = 
        298.257223;
  
semiMinorAxis = 
        6356752  
          Meter;
  
firstEccentricitySquared = 
        6.69437999014 10^-3;
  
secondEccentricitySquared = 6.73949674228 10^-3;

note e is same by elliptic e eq. found in math books

e = firstEccentricitySquared^(1/2);
a = semiMajorAxis;
b = semiMinorAxis;

i can't find any e v. e'.  apparently a blip on radar...

Module[{a = EarthEquatorialRadius, 
        b = EarthPolarRadius}, {1 - (b^2/a^2), (a^2/b^2) - 
          1}]

** "std formula"
	i have NO IDEA where this is from / what it's for

sr[r_,h_,w_]:=(C+h) Cos[w];

(* apparent ambiguity r b ?, see below*)
sz[r_,a_,b_,h_]:=(C b^2/a^2+h)Sin[w];

sC[a_,b_,w_]:=a^2/((a Cos[w])^2+(b Sin[w]^2))^(1/2)

the following refers to lattitude on ellipse
	from M (zenith merid) to
		C(rn,ψ)	geodetic		psi
		C(r,Ψ)	geocentric		Psi
	and cartesian coordinates

note p676 says
	(x-h)^2/a^2 + (y^2/b^2) == 1
	C(h,0) is center of ellipse, back to polar...
	h=-d e^2/(1-e^2)

h is height above surface

flattening = 1/298257223563;

(* semi major axis a *)
(* semi minor axis b *)
b[a_,f_]:=a(1-f);

oops.  2 diff formulas so far...
	a book: e=(a-(k^2)(b/a)^2)^(1/2)
	k=+1 ellipse, k=0 parabola, k=-1 hyperbola
	wiki: e 0<circle<ellipse<1<parabola<1<hyperbola
	angular eccentricity (same, for polar coord)
	e can be expressed in terms of the flattening factor

firstEccentricitySquaredFun[a_, b_] := 
        1 - (b^2/a^2);

wgse = 
        firstEccentricitySquared^(1/2);
fefe = firstEccentricitySquaredFun[a, b]^(1/2);

using p676 proof as a rule...
	e=c/a;
	e^2=c^2/a^2;
	c^2=a^2-b^2;
	c^2/a^2=(a^2-b^2)/a^2
	e^2=(1-(b^2/a^2))

?? local curvature of the ellipsoid along the first vertical
	*  radius of curvature in the prime vertical, N

secondEccentricityFun[a_, b_] := (a^2/b^2) - 
          1;

a book: e`=((1/k^2)(a/b)^2-1)^(1/2)

online papers totally disagree on lat,lon convention AND never no mind sector errors. so let's delineate that...

XYZ data is ECEF , looking at (0,0), +x is toward us, +y to right, z is up
r	is magnitude[{X,Y,Z}]
	const only if xyz is know to be on surface of sphere r
rc	is radius of xy plane given z, mag[{X,Y}]
map	lon data arranged for theta -180,180
	lat lon and (0,0) is {+a,0,0} lat

lat.	If[r===0,0,ArcSin[Z/r]]
lat.	If[rc===0,0,ArcTan[Z/rc]]

lon.	If[Y===0,If[X>=0,0,Pi],Sign[Y](Pi-Sign[Y] 2 ArcTan[(X/Y)])/2],
		(* avoid float call *)
lon.	If[rc===0,0,Sign[y] ArcCos[X/rc]]
lon.	If[mxy===0,0,Sign[Y] ArcCos[X/mag[{X,Y}]]]
		If[mxy===0,0,Sign[Y] ArcCos[{X,Y}.{1,0}/mag[{X,Y}]]]

Shperical	lat,lon data arranged for theta 0,360
			lat lon and (0,0) is {+a,0,0}
lat.		+ Pi/2 for 0,180
lon.		If[Sign[-y]≥0,2Pi,0] + (any of the lon. above)

parts of definitions ignored
	i think this only prevents div/0, though 2 gets full lon map range

(* y of sin is Z from r Sin[W] *)
Module[{x,y},Plot[x=(a/Meter)Cos[t];y=(a/Meter)Sin[t];ArcSin[y/mag[{x,y}]],{t,-Pi/2,Pi/2}]]
(* x will be rc *)
Module[{x,y},Plot[x=(a/Meter)Cos[t];y=(a/Meter)Sin[t];ArcTan[y/x],{t,-Pi/2,Pi/2}]]

Module[{x,y},
  Plot[x=(a/Meter)Cos[t];y=(a/Meter)Sin[t];Sign[y](Pi-Sign[y]2ArcTan[(x/y)])/2,{t,-Pi,Pi}]]
Module[{x,y},
  Plot[x=(a/Meter)Cos[t];y=(a/Meter)Sin[t];Sign[y]ArcCos[x/(a/Meter)],{t,-Pi,Pi}]]
Module[{x,y},
  Plot[x=(a/Meter)Cos[t];y=(a/Meter)Sin[t];Sign[y]ArcCos[{x,y}.{1,0}/mag[{x,y}]],{t,-Pi,Pi}]]
(* interesting.  this looks roughly like sunclock shading *)
Module[{x,y},
  Plot[x=(a/Meter)Cos[t];y=(a/Meter)Sin[t];ArcTan[{x,y}.{1,0}/mag[{x,y}]],{t,-Pi,Pi}]]

l[Yg_,Xg_]:=2 ArcTan[Y/(X+(X^2+Y^2)^(1/2))]
l[Yg_,Xg_]:=2 ArcTan[Y/X]
l[w_]:=ArcTan[(1-e^2)Tan[w]]

length of line from P to z normal to ellipse

rn[a_, e_, w_] := a/(1 - e^2 Sin[w]^2)^(1/2);

rn[a_, b_, w_, "b"] := 
        a^2/(a^2 (Cos[w]^2)  + (b^2 Sin[w]^2))^(1/
                2);

* all 3 papers agreed on this

XYZg[w_, l_, h_, a_, e_] := 
        Module[{Rn}, 
          Rn = rn[a, e, w]; 
{(Rn + h) Cos[w] 
              Cos[l], 
(Rn + h) Cos[w] 
              Sin[l], 
((1 - e^2) Rn + h) 
              Sin[w]}];

φλh	note: i think this estimates Δh by moving Z only "ok" if h is short
		  however my only test so far shows it's wrong a=1 b=1/4 if h==0

instructions: start with w == geocentric

wlh[X_, Y_, Z_, a_, e_, w_] := 
        Module[{Rn, rc}, 
rc = mag[{X, Y}]; 
          Rn = rn[a, e, 
              w]; 
{ArcTan[(Z + e^2 Rn Sin[w])/
                rc], 
Sign[Y] 
              ArcCos[X/rc], 
(rc/Cos[w]) - 
              Rn}];

wlhNewt[X_, Y_, Z_, a_, e_, 
        count__:  4] := 
Module[{r, rc, rc2, w, wi, W, L, 
          h, i, W2}, 
r = 
          mag[{X, Y, Z}]; 
rc = 
          mag[{X, Y}]; 
W2 = ArcTan[rc/Z];    
{W, L, 
            r} = WLr[X, Y, Z];    
w = 
          GeocentricToGeodetic[W, e, 1, 0]; 
w = 
          ArcTan[rc/Z];    
wi = 
          w; 
Do[
h = 
            rc/Cos[w] - rn[a, e, w]; 
w = 
            ArcTan[(Z/
                    rc) (1 - 
                      e^2 (rn[a, e, 
                              w]/(rn[a, e, w] + h)))^-1];   
, {i, 1, 
            count}]; 
h = 
          rc/Cos[w] - 
            rn[a, e, w]; 
rc2 = (rn[a, e, w] + h) 
            Cos[W]; 
Zn = (rn[a, e, w] + h) Sin[
              w]; 
RnZn = 
          e^2 rn[a, e, w];   
Print@{"r", r, "rc", 
            rc, "rn", rn[a, e, w], "WLr", {W, L, r}, "W-W2", 
            W - W2, "Zn", Zn, "RnZ", 
            RnZn, "Zn+RnZn", (Zn + RnZn), "rc2-rc", (rc2 - 
                rc), "init w", wi, "w-wi", w - wi, "w", 
            w (180/Pi), "w", w, "h", h}; 
{w, 
          L, h}
]

wlhNewt[X,Y,Z,a/Meter,e]

{"r", 6.359429057679795`*^6, "rc", 2.255011962018411`*^6, "rn", 
      6.396915483175567`*^6, 
      "WLr", {1.2083170864753618`, 0.7853981633974483`, 
        6.359429057679795`*^6}, "W-W2", 0.8458378461558269`, "Zn", 
      5.986271396347033`*^6, "RnZ", 42823.38300918727`, "Zn+RnZn", 
      6.02909477935622`*^6, "rc2-rc", 13292.34383018734`, "init w", 
      1.2105384889370825`, "w-wi", 9.560441327494118`*^-11, "w", 
      69.35874635971666`, "w", 1.210538489032687`, 
      "h", -0.27507965452969074`}

{1.21054,0.785398,-0.27508}

?? h should be 0, and even if iterated set to 0 
	h=rc/Cos[w]-rn[a,e,w] != 0

in this neighborhood look for about .13 deg diff

{P,GeocentricToGeodetic[69(Pi/180),e,1,0](180/Pi)}

{{1.59453425`*^6, 1.59453425`*^6, 5.946197019172506`*^6}, 
      69.12843598097585`}

Show[Graphics[
    {Circle[{0,0},{a/Meter,b/Meter}],
      Line[{{0,0},{mag[{X,Y}],Z}}],
      Line[{{0,0},{Z,Z}}]
      }
    ],AspectRatio->1]

⁃Graphics⁃

source:
	Heiskanen and Moritz (1967), submitter (2014), others
topic:
	XYZ to geodetic angle for WGS-84 earth ellipse
note:
	many papers sugg. to find wlh loop is req.
	WGS-84 inexact e causes h to be maybe .2 meter wrong
	(that w and h cannot both be known by one calc.)
	not tested.  likely obsoleted by arbitrary precision
	(~ femto meter diff between this(1 or 4 loops) and others shown, by arb)
	"This converges in a few iterations (4 at most) to a few centimeters."
	"This equation for h diverges at the poles." (fixup omited)
key:
	(XYZ is ECEF, +X toward you, +Y right, +Z up, (0,0,0) is GMT) 
	h is geod. height above ellipse, if any
	rn is line from pt on ellipse to Z axis (below z==0)
	(rn is radius of curvature)
	Zn is height from rn origin to Z

{W,L,r}=WLr[X,Y,Z]; 
w=GeocentricToGeodetic[W,e,1,0];

wlhNewt[X_, Y_, Z_, a_, e_, 
        count__:  4] := 
Module[{r, rc, rc2, w, W, L, h, i, 
          W2}, 
r = mag[{X, Y, Z}]; 
        rc = mag[{X, Y}]; 
w = ArcTan[rc/Z];   
Do[
h = 
            rc/Cos[w] - rn[a, e, w]; 
w = 
            ArcTan[(Z/
                    rc) (1 - 
                      e^2 (rn[a, e, 
                              w]/(rn[a, e, w] + 
                                h)))^-1], 
{i, 1, 
            count}]; 
h = 
          rc/Cos[w] - rn[a, e, w]; 
{W, L, r} = 
          WLr[X, Y, Z]; 
rc2 = (rn[a, e, w] + h) 
            Cos[W]; 
Zn = (rn[a, e, w] + h) Sin[
              w]; 
RnZn = 
          e^2 rn[a, e, w];   
{w, 
          L, h}]

wlh[X_, Y_, Z_, a_, e_, w_] := 
        Module[{Rn, rc}, 
rc = mag[{X, Y}]; 
          Rn = rn[a, e, 
              w]; 
{ArcTan[(Z + e^2 Rn Sin[w])/
                rc], 
Sign[Y] 
              ArcCos[X/rc], 
(rc/Cos[w]) - 
              Rn}];

*** geocentric to XYZ, see ΦΛr
	one book has params in reverse order !

XYZe[W_, L_, 
          r_] := 
{r Cos[W] Cos[
              L], 
r Cos[W] Sin[
              L], 
r Sin[W]};

XYZe[W_, L_, r_, 
          h_] := 
{(r + h) Cos[W] Cos[
              L], 
(r + h) Cos[W] Sin[
              L], 
(r + h) Sin[W]};

PointParametricPlot3D[{Cos[u] Cos[v],Cos[u] Sin[v],(1/4)Sin[u]},{u,-Pi/2,Pi/2},
  {v,-Pi,Pi},PlotPoints->50]
PointParametricPlot3D[{Cos[v] Sin[u],Sin[v] Sin[u],1/4 Cos[u]},{u,-Pi,Pi},{v,-Pi/2,Pi/2},PlotPoints->50]

WLr[X_, Y_, Z_] := 
        Module[{r}, 
r = 
            mag[{X, Y, Z}]; 
{If[r === 0, 0, 
              ArcSin[Z/r]], 
If[Y === 0, 
              If[X >= 0, 0, 
                Pi], 
Sign[
                  Y] (Pi - Sign[Y] 2 ArcTan[(X/Y)])/
                  2], 
r}];

WLr[X_, Y_, Z_, 
          r_] := 
{If[r === 0, 0, 
            ArcSin[Z/r]], 
If[Y === 0, 
            If[X >= 0, 0, 
              Pi], 
Sign[
                Y] (Pi - Sign[Y] 2 ArcTan[(X/Y)])/
                2], 
r};

note: if h is 0 rn doesn't matter

GeodeticToGeocentric[w_, e_, rn_, h_] := 
        If[Abs[w] === Pi/2, w, 
          ArcTan[(1 - e^2 (rn/(rn + h))) Tan[w]]];

GeocentricToGeodetic[W_, e_, rn_, h_] := 
        If[Abs[W] === Pi/2, W, 
          ArcTan[(1 - e^2 (rn/(rn + h)))^-1 
              Tan[W]]];

GeocentricToGeodetic2[{X_, Y_, Z_}, a_, b_] := 
        Module[{gradFP}, 
gradFP = {2  X/a^2, 2  Y/a^2, 
              2  Z/b^2};   
Pi/2 - 
            ArcCos[unit[gradFP] . {0, 0, 1}]];

there are many, many series (ie, binomial) to approx elliptic arc len
	UTM uses a besel curve based one

m (ψ) = 
      b (∫_0^(π/2))(Sqrt[1 - 
                  e^2  Sin[ϕ]^2]) dϕ

unsure if this is "better" or absolute.  apparently mm will take EllipticE anyhow

\(a?E \((e)\) = \(\(a\) \(∫\_0\%\(π/2\)\(\@\(1 - 
                  er2^2  Sin[β]^2\)\) dβ\)\(\ \)\)\)
er2 = e^2/(1 - e^2);  
    tan (β) = (1 - f) tan (ψ)

approxPolarDist[a_, e_, ul_, ll_] := 
      Module[{w}, 
        a Integrate[(1 - e^2  Sin[w]^2)^(1/2), {w, ll, 
              ul}]]

approxPolarDist2[a_, e_, ul_, ll_] := 
      Module[{w}, (ul - 
              ll) ((a^(3/2) + b^(3/2))/2)^(2/
                3)]

{2 approxPolarDist[a, e, Pi, 0], 
      2 approxPolarDist2[a, e, Pi, 0]}

(2 Pi a)

4 SpheroidalDistance[{0, 0}, {90, 0}, SemimajorAxis -> a/Meter, 
          Eccentricity -> e] // N

test point

a

used in quick def test just below

Module[{X, Y, Z, c, A, B, P}, X = a/(4 Meter); 
      Y = a/(4 Meter); c = 1; A = a/Meter; 
      B = b/Meter; 
P = {X, Y, 
          Z /. Solve[X^2/A^2 + Y^2/A^2 + Z^2/B^2 == 1, 
                Z][[1]]}; 
{X, Y, Z} = 
        P; 
Z = -P[[3]]; 
P = {X, 
          Y, Z}; 
Print@{{Z}, 
          X^2/A^2 + Y^2/A^2 + 
            Z^2/B^2}; 
P
]

P = %;
{X, Y, Z} = P;

φλh  m1-3  definitions

Clear[wlhEDm3, wlhEDm1, wlhm1, wlhEDm2, ellipseDilation]

note: ed is an intermed. mult. and can have any length acc justified

ellipseDilation[{X_, Y_, Z_}, a_, 
          b_] := 
((b^2 X^2 + b^2 Y^2 + 
                  a^2 Z^2)^(1/2))/(a b);

ellipseDilation[ .5 {X, Y, Z}, a/Meter, b/Meter]

{mag[P], mag[((P/1.2))], mag[P] - mag[((P/1.2))], 
      mag[(P - (P/1.2))]}

{mag[P], mag[((P/ .9))], mag[P] - mag[((P/ .8))], 
      mag[(P - (P/ .8))]}

wlhEDm1[X_, Y_, Z_, a_, b_, 
        e_] := 
Module[{P, We, Le, re, w, l, h, 
          ed}, 
P = {X, Y, Z}; 
{We, Le, 
            re} = WLr[X, Y, Z]; 
w = 
          GeocentricToGeodetic[We, e, 1, 0]; 
ed = 
          ellipseDilation[P, a, b]; 
h = 
          If[ed === 0, 0, mag[P] - mag[P/ed]]; 
{w, Le, 
          h}]

wlhEDm1[X, Y, Z, a/Meter, b/Meter, e]

% {(180/Pi), (180/Pi), 1}

wlhEDm2[X_, Y_, Z_, a_, b_, 
        e_] := 
Module[{P, We, Le, re, w, l, h, 
          ed}, 
P = {X, Y, Z}; 
ed = 
          ellipseDilation[P, a, b]; 
{We, Le, re} = 
          WLr[X/ed, Y/ed, Z/ed]; 
w = 
          GeocentricToGeodetic[We, e, 1, 0]; 
h = 
          If[ed === 0, 0, mag[P] - mag[P/ed]]; 
{w, Le, 
          h}]

wlhEDm2[X, Y, Z, a/Meter, b/Meter, e]

% {(180/Pi), (180/Pi), 1}

note for Le that ungfp does not share origin with {1,0,0}
	gradFP={2X/a^2,2Y/a^2,2Z/b^2};(*surf N vec at P*)
	w=Pi/2-ArcCos[unit[gradFP].{0,0,1}];

see chart in eq. def. above for Le

wlhEDm3[X_, Y_, Z_, a_, b_, 
        e_] := 
Module[{P, Le, Pe, w, l, h, gradFP, ed, 
          mxy}, 
P = {X, Y, Z}; 
        mxy = mag[{X, Y}]; 
ed = 
          ellipseDilation[P, a, b]; 
w = 
          GeocentricToGeodetic2[P, a, b]; 
Le = 
          If[mxy === 0, 0, Sign[Y] ArcCos[X/mxy]]; 
h = 
          If[ed === 0, 0, mag[P] - mag[P/ed]]; 
{w, Le, 
          h}]

wlhEDm3[X, Y, Z, a/Meter, b/Meter, e]

% {(180/Pi), (180/Pi), 1}

Module[{X, Y, Z, gradFP, 
        ungfp}, {X, Y, Z} = { .2 + Pi/4,  .2 + Pi/4, 
          0}; 
gradFP = {2  X/(a/Meter)^2, 
          2  Y/(a/Meter)^2, 
          2  Z/(b/
                    Meter)^2};   
ungfp = 
        gradFP/mag[
            gradFP]; 
ungfp; 
{Sign[
            Y] ArcCos[X/mag[{X, Y, Z}]], ArcCos[ungfp . {1, 0, 0}], 
        Sign[Y] ArcCos[X/mag[{X, Y, Z}]] - 
          ArcCos[ungfp . {1, 0, 0}]}
]

Fxyz[X_,Y_,Z_,a_,b_]:={2X/a^2,2Y/a^2,2Z/b^2};
w=Pi/2-ArcCos[gradFP.{0,0,1}/(mag[gradFP]mag[{1,0,0}])];
Le=ArcCos[gradFP.{1,0,0}/(mag[gradFP]mag[{0,1,0}])];
Le=Sign[Y] ArcCos[X/mag[{X,Y,Z}]];

------------------------------
show eD method pictorially

assume line rn is normal to inner ellipse earth (perpandicular, standing on), then it is normal on dilated one.  XYZ gives pt on dilated one so only dilation is solved for (see eD), then h ("the hard part") is just dist between two points.  doesn't seem recursion is needed.  this method is simple to learn and understand.  note this has NOT been not verified for satellite use

Show[
      Graphics[{Circle[{0, 0}, {7, 1}], Circle[{0, 0}, 1.2 {7, 1}], 
          Line[{{0, -3/4}, 1.5 {1, 1}}], Line[{{-7, 0}, {7, 0}}], 
          Line[{{0, -1}, {0, 1}}], RGBColor[1, 0, 0], PointSize[ .025], 
          Point[{1.2, 1}], Point[{1.3, 1.2}]}]]

--------------------------
many online papers show rn, h, ellipse for newton method pictorially, so search for that