GPS

$ % macros % MathJax \newcommand{\bigtimes}{\mathop{\vcenter{\hbox{$\Huge\times\normalsize$}}}} % prefix with \! and postfix with \!\! % sized grouping symbols \renewcommand {\aa} [1] {\left\langle {#1} \right\rangle} % <> angle brackets \newcommand {\bb} [1] {\left[ {#1} \right]} % [] brackets \newcommand {\cc} [1] {\left\{ {#1} \right\}} % {} curly braces \newcommand {\mm} [1] {\lVert {#1} \rVert} % || || double norm bars \newcommand {\nn} [1] {\lvert {#1} \rvert} % || norm bars \newcommand {\pp} [1] {\left( {#1} \right)} % () parentheses % unit \newcommand {\unit} [1] {\bb{\mathrm{#1}}} % measurement unit % math \newcommand {\fn} [1] {\mathrm{#1}} % function name % sets \newcommand {\setZ} {\mathbb{Z}} \newcommand {\setQ} {\mathbb{Q}} \newcommand {\setR} {\mathbb{R}} \newcommand {\setC} {\mathbb{C}} % arithmetic \newcommand {\q} [2] {\frac{#1}{#2}} % quotient, because fuck \frac % trig \newcommand {\acos} {\mathrm{acos}} % \mathrm{???}^{-1} is misleading \newcommand {\asin} {\mathrm{asin}} \newcommand {\atan} {\mathrm{atan}} \newcommand {\atantwo} {\mathrm{atan2}} % at angle = atan2(y, x) \newcommand {\asec} {\mathrm{asec}} \newcommand {\acsc} {\mathrm{acsc}} \newcommand {\acot} {\mathrm{acot}} % complex numbers \newcommand {\z} [1] {\tilde{#1}} \newcommand {\conj} [1] {{#1}^\ast} \renewcommand {\Re} {\mathfrak{Re}} % real part \renewcommand {\Im} {\mathrm{I}\mathfrak{m}} % imaginary part \newcommand {\abs} [1] {\nn{#1}} % quaternions \newcommand {\quat} [1] {\tilde{\mathbf{#1}}} % quaternion symbol \newcommand {\uquat} [1] {\check{\mathbf{#1}}} % versor symbol \newcommand {\gquat} [1] {\tilde{\boldsymbol{#1}}} % greek quaternion symbol \newcommand {\guquat}[1] {\check{\boldsymbol{#1}}} % greek versor symbol % vectors \renewcommand {\vec} [1] {\mathbf{#1}} % vector symbol \newcommand {\uvec} [1] {\hat{\mathbf{#1}}} % unit vector symbol \newcommand {\gvec} [1] {\boldsymbol{#1}} % greek vector symbol \newcommand {\guvec} [1] {\hat{\boldsymbol{#1}}} % greek unit vector symbol % special math vectors \renewcommand {\r} {\vec{r}} % r vector [m] \newcommand {\R} {\vec{R}} % R = r - r' difference vector [m] \newcommand {\ur} {\uvec{r}} % r unit vector [#] \newcommand {\uR} {\uvec{R}} % R unit vector [#] \newcommand {\ux} {\uvec{x}} % x unit vector [#] \newcommand {\uy} {\uvec{y}} % y unit vector [#] \newcommand {\uz} {\uvec{z}} % z unit vector [#] \newcommand {\urho} {\guvec{\rho}} % rho unit vector [#] \newcommand {\utheta} {\guvec{\theta}} % theta unit vector [#] \newcommand {\uphi} {\guvec{\phi}} % phi unit vector [#] \newcommand {\un} {\uvec{n}} % unit normal vector [#] % vector operations \newcommand {\inner} [2] {\left\langle {#1} , {#2} \right\rangle} % <,> \newcommand {\outer} [2] {{#1} \otimes {#2}} \newcommand {\norm} [1] {\mm{#1}} \renewcommand {\dot} {\cdot} % dot product \newcommand {\cross} {\times} % cross product % matrices \newcommand {\mat} [1] {\mathbf{#1}} % matrix symbol \newcommand {\gmat} [1] {\boldsymbol{#1}} % greek matrix symbol % ordinary derivatives \newcommand {\od} [2] {\q{d #1}{d #2}} % ordinary derivative \newcommand {\odn} [3] {\q{d^{#3}{#1}}{d{#2}^{#3}}} % nth order od \newcommand {\odt} [1] {\q{d{#1}}{dt}} % time od % partial derivatives \newcommand {\de} {\partial} % partial symbol \newcommand {\pd} [2] {\q{\de{#1}}{\de{#2}}} % partial derivative \newcommand {\pdn} [3] {\q{\de^{#3}{#1}}{\de{#2}^{#3}}} % nth order pd \newcommand {\pdd} [3] {\q{\de^2{#1}}{\de{#2}\de{#3}}} % 2nd order mixed pd \newcommand {\pdt} [1] {\q{\de{#1}}{\de{t}}} % time pd % vector derivatives \newcommand {\del} {\nabla} % del \newcommand {\grad} {\del} % gradient \renewcommand {\div} {\del\dot} % divergence \newcommand {\curl} {\del\cross} % curl % differential vectors \newcommand {\dL} {d\vec{L}} % differential vector length [m] \newcommand {\dS} {d\vec{S}} % differential vector surface area [m^2] % special functions \newcommand {\Hn} [2] {H^{(#1)}_{#2}} % nth order Hankel function \newcommand {\hn} [2] {H^{(#1)}_{#2}} % nth order spherical Hankel function % transforms \newcommand {\FT} {\mathcal{F}} % fourier transform \newcommand {\IFT} {\FT^{-1}} % inverse fourier transform % signal processing \newcommand {\conv} [2] {{#1}\ast{#2}} % convolution \newcommand {\corr} [2] {{#1}\star{#2}} % correlation % abstract algebra \newcommand {\lie} [1] {\mathfrak{#1}} % lie algebra % other \renewcommand {\d} {\delta} % optimization %\DeclareMathOperator* {\argmin} {arg\,min} %\DeclareMathOperator* {\argmax} {arg\,max} \newcommand {\argmin} {\fn{arg\,min}} \newcommand {\argmax} {\fn{arg\,max}} % waves \renewcommand {\l} {\lambda} % wavelength [m] \renewcommand {\k} {\vec{k}} % wavevector [rad/m] \newcommand {\uk} {\uvec{k}} % unit wavevector [#] \newcommand {\w} {\omega} % angular frequency [rad/s] \renewcommand {\TH} {e^{j \w t}} % engineering time-harmonic function [#] % classical mechanics \newcommand {\F} {\vec{F}} % force [N] \newcommand {\p} {\vec{p}} % momentum [kg m/s] % \r % position [m], aliased \renewcommand {\v} {\vec{v}} % velocity vector [m/s] \renewcommand {\a} {\vec{a}} % acceleration [m/s^2] \newcommand {\vGamma} {\gvec{\Gamma}} % torque [N m] \renewcommand {\L} {\vec{L}} % angular momentum [kg m^2 / s] \newcommand {\mI} {\mat{I}} % moment of inertia tensor [kg m^2/rad] \newcommand {\vw} {\gvec{\omega}} % angular velocity vector [rad/s] \newcommand {\valpha} {\gvec{\alpha}} % angular acceleration vector [rad/s^2] % electromagnetics % fields \newcommand {\E} {\vec{E}} % electric field intensity [V/m] \renewcommand {\H} {\vec{H}} % magnetic field intensity [A/m] \newcommand {\D} {\vec{D}} % electric flux density [C/m^2] \newcommand {\B} {\vec{B}} % magnetic flux density [Wb/m^2] % potentials \newcommand {\A} {\vec{A}} % vector potential [Wb/m], [C/m] % \F % electric vector potential [C/m], aliased % sources \newcommand {\I} {\vec{I}} % line current density [A] , [V] \newcommand {\J} {\vec{J}} % surface current density [A/m] , [V/m] \newcommand {\K} {\vec{K}} % volume current density [A/m^2], [V/m^2] % \M % magnetic current [V/m^2], aliased, obsolete % materials \newcommand {\ep} {\epsilon} % permittivity [F/m] % \mu % permeability [H/m], aliased \renewcommand {\P} {\vec{P}} % polarization [C/m^2], [Wb/m^2] % \p % electric dipole moment [C m], aliased \newcommand {\M} {\vec{M}} % magnetization [A/m], [V/m] \newcommand {\m} {\vec{m}} % magnetic dipole moment [A m^2] % power \renewcommand {\S} {\vec{S}} % poynting vector [W/m^2] \newcommand {\Sa} {\aa{\vec{S}}_t} % time averaged poynting vector [W/m^2] % quantum mechanics \newcommand {\bra} [1] {\left\langle {#1} \right|} % <| \newcommand {\ket} [1] {\left| {#1} \right\rangle} % |> \newcommand {\braket} [2] {\left\langle {#1} \middle| {#2} \right\rangle} $

My unpublished GPS notes.

Geometry

Conic Sections Dandelin spheres show that an ellipse is the locus of points with the property that the sum of distances from two foci is constant. $$ L = \norm{\r - \vec{f}_-} + \norm{\r - \vec{f}_+} $$ let $\vec{f}_{\pm} = \pm \ux c$ semi-major axis $a$ semi-minor axis $b$ linear eccentricity $c$ $$ L = c + a + [a - c] = 2 a $$ $$ \frac{x^2}{a^2} + \frac{y^2}{b^2} = 1 $$ $$ b^2 = a^2 - c^2 $$
$$ 2 a = \sqrt{[x + c]^2 + y^2} + \sqrt{[x - c]^2 + y^2} $$ square this to get $$ \begin{split} 4 a^2 &= [x + c]^2 + y^2 + 2 \sqrt{[x + c]^2 + y^2} \sqrt{[x - c]^2 + y^2} + [x - c]^2 + y^2 \\ 4 a^2 &= 2 x^2 + 2 y^2 + 2 c^2 + 2 \sqrt{[[x + c]^2 + y^2] [[x - c]^2 + y^2]} \\ 2 a^2 &= x^2 + y^2 + c^2 + \sqrt{[x^2 + 2 x c + c^2 + y^2] [x^2 - 2 x c + c^2 + y^2]} \\ 2 a^2 &= x^2 + y^2 + c^2 + \sqrt{x^4 - 2 x^3 c + x^2 c^2 + x^2 y^2 + 2 x^3 c - 4 x^2 c^2 + 2 x c^3 + 2 x y^2 c + x^2 c^2 - 2 x c^3 + c^4 + y^2 c^2 + x^2 y^2 - 2 x y^2 c + y^2 c^2 + y^4} \\ 2 a^2 &- x^2 - y^2 - c^2 = \sqrt{x^4 + y^4 + c^4 - 2 x^2 c^2 + 2 x^2 y^2 + 2 y^2 c^2} \end{split} $$ square this to get $$ \begin{split} 4 a^4 &- 4 a^2 x^2 - 4 a^2 y^2 - 4 a^2 c^2 + x^4 + 2 x^2 y^2 + 2 x^2 c^2 + y^4 + 2 y^2 c^2 + c^4 = x^4 + y^4 + c^4 - 2 x^2 c^2 + 2 x^2 y^2 + 2 y^2 c^2 \\ 4 a^4 &- 4 a^2 x^2 - 4 a^2 y^2 - 4 a^2 c^2 + 4 x^2 c^2 = 0 \\ a^2 [a^2 - c^2] &= [a^2 - c^2] x^2 + a^2 y^2 \\ a^2 b^2 &= b^2 x^2 + a^2 y^2 \\ 1 &= \frac{x^2}{a^2} + \frac{y^2}{b^2} \end{split} $$
$$ \begin{bmatrix} x \\ y \end{bmatrix} = \begin{bmatrix} a \cos \theta \\ b \sin \theta \end{bmatrix} $$ eccentricity $e$ $$ e \equiv \frac{c}{a} = \sqrt{1 - \frac{b^2}{a^2}} $$ semi-latus rectum $p$ Let $x = \pm c = \pm a e$. Then $$ \frac{[\pm a e]^2}{a^2} + \frac{y^2}{b^2} = 1 $$ thus $$ y = \pm b \sqrt{1 - e^2} $$ $$ p = b \sqrt{1 - e^2} = a [1 - e^2] = \frac{b^2}{a} $$

Polar Forms

Around the center of an ellipse $$ \boxed{r = a [1 - e \cos \theta]} $$

$$ \begin{split} r &= \sqrt{[x - c]^2 + [y - 0]^2} \\ &= \sqrt{[a \cos \theta - c]^2 + b^2 \sin^2 \theta} \\ &= \sqrt{[a \cos \theta - a e]^2 + a^2 [1 - e^2] \sin^2 \theta} \\ &= \sqrt{a^2 \cos^2 \theta - 2 a^2 e \cos \theta + a^2 e^2 + a^2 \sin^2 \theta - a^2 e^2 \sin^2 \theta} \\ &= a \sqrt{1 - 2 e \cos \theta + e^2 [1 - \sin^2 \theta]} \\ &= a \sqrt{[1 - e \cos \theta]^2} \\ &= a [1 - e \cos \theta] \end{split} $$

Around a focus $$ \boxed{r = \frac{p}{1 + e \cos \theta}} $$

$$ 2 a = \norm{\r} + \norm{\r - \ux 2 c} $$ $$ 2 a = r + \sqrt{[r \cos \theta - 2 c]^2 + [r \sin \theta]^2} $$ $$ [2 a - r]^2 = [r \cos \theta - 2 c]^2 + [r \sin \theta]^2 $$ $$ 4 a^2 - 4 a r + r^2 = r^2 \cos^2 \theta - 4 c r \cos \theta + 4 c^2 + r^2 \sin^2 \theta $$ $$ a^2 - c^2 = r [a - c \cos \theta] $$ $$ r = \frac{b^2}{a - c \cos \theta} = \frac{\frac{b^2}{a}}{1 - \frac{c}{a} \cos \theta} = \frac{p}{1 - e \cos \theta} $$

Orbital Mechanics

eccentric anomaly $E$ mean motion $n$ mean anomaly $M$ Kepler's Laws
  1. Planet orbits trace out ellipses with sun at one focus
  2. Planet orbits sweep out equal area in equal time
  3. $T^2 \propto a^3$
Newtonian mechanics: $$ \F = \odt{\p} $$ $$ \F = m \a $$ $$ \F = -G \frac{m m'}{R^2} \uR $$ where $G = 6.672e-11 [N m^2/kg^2]$ is the gravitational constant and $\R \equiv \r - \r'$ points from source to field position $$\F = -\grad U$$ $$ U(\r) = -G \frac{m'}{R} $$

Kepler Problem

$$ \begin{cases} \F_{12} = m_1 \a_1 = -G \frac{m_1 m_2}{R_{12}^2} \uR_{12} \\ \F_{21} = m_2 \a_2 = -G \frac{m_2 m_1}{R_{21}^2} \uR_{21} \end{cases} $$ $$\F_{12} = -\F_{21}$$ $$ \r_{CM} = \frac{m_1 \r_1 + m_2 \r_2}{m_1 + m_2} $$ $M = m_1 + m_2$ $$ m_1 \a_1 + m_2 \a_2 = \vec{0} $$ $$ \a_{CM} = \vec{0} $$ $$ \v_{CM} = \v_{CM,0} $$ $$ \L = \r_1 \cross \p_1 + \r_2 \cross \p_2 \\ $$ $$ \mathcal{L} = K - U $$ $$ \mathcal{L}(t, \r_1, \v_1, \r_2, \v_2) = \frac{1}{2} m_1 v_1^2 + \frac{1}{2} m_2 v_2^2 - U(\r_1, \r_2) $$ $$ \mu = \frac{m_1 m_2}{m_1 + m_2} $$ $$ \r_{21} = \r_2 - \r_1 $$ $$ \r_1 = \r_{CM} - \frac{m_2}{M} \r_{21} $$ $$ \r_2 = \r_{CM} + \frac{m_1}{M} \r_{21} $$ $$ \mathcal{L}(t, \r_{CM}, \v_{CM}, \r_{21}, \v_{21}) = \frac{1}{2} M v_{CM}^2 + \frac{1}{2} \mu v_{21}^2 - U(r_{21}) $$ $$ \begin{split} K &=\frac{1}{2} M v_{CM}^2 + \frac{1}{2} \mu v_{21}^2 \\ &= \frac{1}{2} M \frac{\sum_i [m_1 v_{1,i} + m_2 v_{2,i}]^2}{M^2} + \frac{1}{2} \frac{m_1 m_2}{M} \sum_i [v_{2,i} - v_{1,i}]^2 \\ &= \frac{1}{2 M} \sum_i [m_1^2 v_{1,i}^2 + 2 m_1 m_2 v_{1,i} v_{2,i} + m_2^2 v_{2,i}^2 + m_1 m_2 v_{2,i}^2 - 2 m_1 m_2 v_{1,i} v_{2,i} + m_1 m_2 v_{1,i}^2 ] \\ &= \frac{1}{2 M} \sum_i [[m_1 + m_2] m_1 v_{1,i}^2 + [m_1 + m_2] m_2 v_{2,i}^2] \\ &= \frac{1}{2} \sum_i [m_1 v_{1,i}^2 + m_2 v_{2,i}^2] \\ &= \frac{1}{2} m_1 v_1^2 + \frac{1}{2} m_2 v_2^2 \end{split} $$ $$\mathcal{L} = \mathcal{L}_{CM} + \mathcal{L}_{rel}$$ As shown, the CM is inertial, so choose to locate the origin at the CM which makes the CM kinect energy term zero $$ \mathcal{L}(t, \vec{0}, \vec{0}, \r_{21}, \v_{21}) = \frac{1}{2} \mu v_{21}^2 - U(r_{21}) $$ This is Lagrangian is identical in form to a single particle in an inverse square-law field.

Using the TLE

The TLE gives mean motion, epoch, and mean anomaly at epoch, while the GPS measurement yields time $$ \boxed{n \equiv \frac{2 \pi}{T}} $$ compute mean anomaly as $$ \boxed{M = M_0 + n [t - t_0]} $$ Compute semi-major axis given mean motion and Earth's gravitational constant using Kepler's Law $$ \boxed{a = \bb{\frac{G [m + M_\oplus]}{n^2}}^{\frac{1}{3}} \approx \bb{\frac{G M_\oplus}{n^2}}^{\frac{1}{3}}} $$ Note, the earth gravitational constant is known to more significant figures than either factor. Compute semi-minor axis given semi-major axis and eccentricity $$ \boxed{b = a \sqrt{1 - e^2}} $$ Compute linear eccentricity given semi-major axis and eccentricity $$ \boxed{c = a e} $$ Compute eccentric anomaly given mean anomaly and eccentricity using Kepler's equation $$ \boxed{M = E - e \sin E} $$ This equation can't be solved analytically. Instead, solve it iteratively. For instance, use Newton's method. Another way proposed by [1] is to iterate over $$ E_{i+1} = E_i + \Delta E $$ $$ \Delta M = M - [E_i - e \sin E_i] $$ $$ \Delta E = \frac{\Delta M}{1 - e \cos E_i} $$ $E_0 = 0.0$ $\Delta E = M - e \sin M$ $\epsilon = 1.0e-8$ [rad] test $\abs{\Delta E} ≤ \epsilon$ test smaller than max iterations Finally, compute satellite position given semi-major axis, semi-minor axis, linear eccentricity, and eccentric anomaly $$ \boxed{\r = \begin{bmatrix} a \cos E - c \\ b \sin E \\ 0 \\ 1 \end{bmatrix}} $$

Geodetic Stuff

parametric latitude $\beta$ geocentric latitude $\theta$ geodetic latitude $\phi$ geodetic longitude $\lambda$ geodetic height $h$ flattening $f$ $$ f \equiv 1 - \sqrt{1 - e^2} $$

ecef <-> geodetic

$(\phi, \lambda, h) \rightarrow (x, y, z)$ $$ \boxed{ \begin{cases} x = [K a^2 + h] \cos \phi \cos \lambda \\ y = [K a^2 + h] \cos \phi \sin \lambda \\ z = [K b^2 + h] \sin \phi \end{cases} \quad\text{where}\quad K = \frac{1}{\sqrt{a^2 \cos^2 \phi + b^2 \sin^2 \phi}} } $$
$$ \frac{x^2}{a^2} + \frac{z^2}{b^2} = 1 $$ $$ \frac{2 x}{a^2} + \frac{2 z}{b^2} \pd{z}{x} = 0 $$ $$ \pd{z}{x} = - \frac{b^2}{a^2} \frac{x}{z} $$ $$ \pd{z}{x} = -\tan \pp{\frac{\pi}{2} - \phi} = -\cot \phi $$ $$ a^2 z \cos \phi = b^2 x \sin \phi $$ $$ a^4 z^2 \cos^2 \phi = b^4 x^2 \sin^2 \phi $$ $$ b^2 x^2 + a^2 z^2 = a^2 b^2 $$ $$ \begin{cases} b^4 x^2 \sin^2 \phi + a^2 b^2 z^2 \sin^2 \phi = a^2 b^4 \sin^2 \phi \\ a^2 b^2 x^2 \cos^2 \phi + a^4 z^2 \cos^2 \phi = a^4 b^2 \cos^2 \phi \end{cases} $$ $$ \begin{cases} a^2 z^2 [a^2 \cos^2 \phi + b^2 \sin^2 \phi] = a^2 b^4 \sin^2 \phi \\ b^2 x^2 [a^2 \cos^2 \phi + b^2 \sin^2 \phi] = a^4 b^2 \cos^2 \phi \end{cases} $$ $$ \begin{cases} z \sqrt{a^2 \cos^2 \phi + b^2 \sin^2 \phi} = b^2 \sin \phi \\ x \sqrt{a^2 \cos^2 \phi + b^2 \sin^2 \phi} = a^2 \cos \phi \end{cases} $$ $$ \begin{cases} z = K b^2 \sin \phi \\ x = K a^2 \cos \phi \end{cases} $$ where $K = [a^2 \cos^2 \phi + b^2 \sin^2 \phi]^{-\frac{1}{2}}$ $$ \begin{cases} x = K a^2 \cos \phi \cos \lambda \\ y = K a^2 \cos \phi \sin \lambda \\ z = K b^2 \sin \phi \end{cases} $$ $$ \begin{cases} x = [K a^2 + h] \cos \phi \cos \lambda \\ y = [K a^2 + h] \cos \phi \sin \lambda \\ z = [K b^2 + h] \sin \phi \end{cases} $$
$(x, y, z) \rightarrow (\phi, \lambda, h)$ $$ \boxed{ \begin{cases} \phi = \atan \pp{ \frac{z}{\sqrt{x^2 + y^2}} \frac{[K a^2 + h]}{[K b^2 + h]} } \\ \lambda = \atantwo(y, x) \\ h = \frac{\sqrt{x^2 + y^2}}{\cos \phi} - K a^2 \end{cases} \quad\text{where}\quad K = \frac{1}{\sqrt{a^2 \cos^2 \phi + b^2 \sin^2 \phi}} } $$ These equations are coupled and must be solved iteratively.
$$ \lambda = \atantwo(y, x) $$ $$ \sqrt{x^2 + y^2} = \sqrt{[K a^2 + h]^2 \cos^2 \phi \cos^2 \lambda + [K a^2 + h]^2 \cos^2 \phi \sin^2 \lambda} = [K a^2 + h] \cos \phi $$ $$ h = \frac{\sqrt{x^2 + y^2}}{\cos \phi} - K a^2 $$ $$ \frac{z}{\sqrt{x^2 + y^2}} = \frac{[K b^2 + h] \sin \phi}{[K a^2 + h] \cos \phi} $$ $$ \phi = \atan \pp{ \frac{z}{\sqrt{x^2 + y^2}} \frac{[K a^2 + h]}{[K b^2 + h]} } $$
THE FOLLOWING IS IN SPHERICAL COORDINATES spherical to rectangular matrix $$ \boxed{ \ur = \begin{bmatrix} \sin \theta \cos \phi \\ \sin \theta \sin \phi \\ \cos \theta \\ 0 \end{bmatrix} \quad \utheta = \begin{bmatrix} \cos \theta \cos \phi \\ \cos \theta \sin \phi \\ -\sin \theta \\ 0 \end{bmatrix} \quad \uphi = \begin{bmatrix} -\sin \phi \\ \cos \phi \\ 0 \\ 0 \end{bmatrix} \quad \r = \begin{bmatrix} r \sin \theta \cos \phi \\ r \sin \theta \sin \phi \\ r \cos \theta \\ 1 \end{bmatrix} } $$ $$ \mat{M}_{rs} = \begin{bmatrix} \ur | \utheta | \uphi | \r \end{bmatrix} $$ unit sphere to spheroid matrix $$ \mat{M}_{spheroid} = \fn{diag}(a, a, b, 1) $$ $$ \mat{N} = [\mat{M}^{-1}]^T $$ $$ \mat{N}_{spheroid} = \fn{diag}(1/a, 1/a, 1/b, 1) $$ spheroid to rectangular matrix $$ \boxed{ \un = \begin{bmatrix} K b \sin \theta \cos \phi \\ K b \sin \theta \sin \phi \\ K a \cos \theta \\ 0 \end{bmatrix} \quad \utheta = \begin{bmatrix} K a \cos \theta \cos \phi \\ K a \cos \theta \sin \phi \\ -K b \sin \theta \\ 0 \end{bmatrix} \quad \uphi = \begin{bmatrix} -\sin \phi \\ \cos \phi \\ 0 \\ 0 \end{bmatrix} \quad \r = \begin{bmatrix} [K b h + a] \sin \theta \cos \phi \\ [K b h + a] \sin \theta \sin \phi \\ [K a h + b] \cos \theta \\ 1 \end{bmatrix} } $$ $$ \mat{M}_{rs} = \begin{bmatrix} \un | \utheta | \uphi | \r \end{bmatrix} $$
$$ \vec{n} = \mat{N}_{spheroid} \ur = \begin{bmatrix} \frac{1}{a} \sin \theta \cos \phi \\ \frac{1}{a} \sin \theta \sin \phi \\ \frac{1}{b} \cos \theta \\ 0 \end{bmatrix} $$ $$ \norm{\vec{n}} = \sqrt{\frac{1}{a^2} \sin^2 \theta \cos^2 \phi + \frac{1}{a^2} \sin^2 \theta \sin^2 \phi + \frac{1}{b^2} \cos^2 \theta} = \frac{1}{ab} \sqrt{a^2 \cos^2 \theta + b^2 \sin^2 \theta} = \frac{1}{Kab} $$ where $K = 1 / \sqrt{a^2 \cos^2 \theta + b^2 \sin^2 \theta}$ $$ \un = \frac{\vec{n}}{\norm{\vec{n}}} = \begin{bmatrix} K b \sin \theta \cos \phi \\ K b \sin \theta \sin \phi \\ K a \cos \theta \\ 0 \end{bmatrix} $$ $$ \vec{\theta} = \mat{M}_{spheroid} \utheta \begin{bmatrix} a \cos \theta \cos \phi \\ a \cos \theta \sin \phi \\ -b \sin \theta \\ 0 \end{bmatrix} $$ $$ \norm{\vec{\theta}} = \sqrt{a^2 \cos^2 \theta \cos^2 \phi + a^2 \cos^2 \theta \sin^2 \phi + b^2 \sin^2 \theta} = \sqrt{a^2 \cos^2 \theta + b^2 \sin^2 \theta} = \frac{1}{K} $$ $$ \utheta = \frac{\vec{\theta}}{\norm{\vec{\theta}}} = \begin{bmatrix} K a \cos \theta \cos \phi \\ K a \cos \theta \sin \phi \\ -K b \sin \theta \\ 0 \end{bmatrix} $$ $$ \vec{\phi} = \mat{M}_{spheroid} \uphi \begin{bmatrix} -a \sin \phi \\ a \cos \phi \\ 0 \\ 0 \end{bmatrix} $$ $$ \norm{\vec{\phi}} = \sqrt{a^2 \sin^2 \phi + a^2 \cos^2 \phi} = a $$ $$ \uphi = \frac{\vec{\phi}}{\norm{\vec{\phi}}} = \begin{bmatrix} -\sin \phi \\ \cos \phi \\ 0 \\ 0 \end{bmatrix} $$ $$ \r = \mat{M}_{spheroid} \r + \un h = \begin{bmatrix} [K b h + a] \sin \theta \cos \phi \\ [K b h + a] \sin \theta \sin \phi \\ [K a h + b] \cos \theta \\ 1 \end{bmatrix} $$
$$ \boxed{\tan \beta = \frac{1}{\sqrt{1 - e^2}} \tan \theta = \sqrt{1 - e^2} \tan \phi} $$
$$ \tan \theta = \frac{b \sin \beta}{a \cos \beta} = \frac{b}{a} \tan \beta $$ $$ \tan \phi = \frac{b \sin \beta - z_0}{a \cos \beta} $$ $$ \r + t \un = \vec{z}_0 $$ $$ \begin{bmatrix} 0 \\ z_0 \end{bmatrix} = \begin{bmatrix} a \cos \beta \\ b \sin \beta \end{bmatrix} + \begin{bmatrix} K b \cos \beta \\ K a \sin \beta \end{bmatrix} t $$ $$ t = -\frac{a}{K b} $$ $$ z_0 = b \sin \beta \bb{1 - \frac{a^2}{b^2}} $$ $$ \tan \phi = \frac{b \sin \beta - b \sin \beta \bb{1 - \frac{a^2}{b^2}}}{a \cos \beta} = \frac{a}{b} \tan \beta $$ $$ \frac{b}{a} = \sqrt{1 - e^2} $$

References