I am trying to understand conversion from geodetic coordinates to ENU coordinates. As per this wikipedia page, its two step process: geodetic to ECEF coordinates and then ECEF coordinates to ENU coordinates. Wikipedia gives the conversion from geodetic to ECEF coordinates as follows:
$$X=(N(\phi)+h)\cos\phi \cos\lambda$$ $$Y=(N(\phi)+h)\cos\phi \sin\lambda$$ $$Z=\left(\frac{b^2}{a^2}N(\phi)+h\right)\sin\phi$$ $$=\left((1-e^2)N(\phi)+h\right)\sin\phi$$ $$=\left((1-f)^2N(\phi)+h\right)\sin\phi$$
where
$$N(\phi)=\frac{a^2}{\sqrt{a^2\cos^2\phi+b^2\sin^2\phi}}$$
$a$ is the equatorial radius (semi-major axis)
$b$ is the polar radius (semi-minor axis)
$f=1-\frac{b}{a}$ is the flattening of the ellipsoid
However, I am not satisfied with the direct formula. I tried deriving cartesian coordinates from spherical coordinates:
But this does not directly map to above problem and hence really does not match with what wikipedia states.
Can someone please help me with the proof of these equations? If not the proofs, the link / source / any resource discussing the same?
PS:
Below diagram may help visualize different quantities ref:
Update
I am able to derive equations for $X$ and $Y$, and a part of $Z$. Will love if someone help me with equation of $Z$ and $N(\phi)$.
This is how I derived it:
$$\cos \lambda=\frac{x}{l}$$ $$\text{(1) }X=l\cos\lambda$$ $$\cos \phi=\frac{l}{N+h}$$ $$\text{(2) }l=(N+h)\cos \phi$$ $$\text{(1),(2)}\rightarrow \text{(3) }X=(N+h)\cos \phi\cos \lambda$$ $$\sin\lambda=\frac{y}{l}$$ $$\text{(4) }y=l\sin\lambda$$ $$\text{(2),(4)}\rightarrow\text{(5) }Y=(N+h)\cos\phi\sin\lambda$$
Sorry for non-conventional numbering $(x)$