GISの分野では位置座標の表現に緯度・経度を用いることが多いのですが、それでもやはり3D表現や宇宙規模の話（衛星測位等）が絡む場面では $(x, y, z)$ の直交座標系が用いられます。
3Dグラフィクスなどの世界にいる方にとっては、経緯度＋高さ よりも $(x, y, z)$ の直交座標のほうがずっと馴染み深いでしょう。
この記事では、世界測地系 WGS 84 の 緯度・経度・高さ（楕円体高）で表される点を、WGS 84 の地心座標系 (x, y, z) に変換する計算をむりやり導いてみます（答え合わせもします）。
技術的にいえば、以下のような変換をします：
-
From: WGS 84 Geographic 3D (lat, lng, height)
EPSG:4979 -
To: WGS 84 Geocentric (x, y, z)
EPSG:4978
地球が真の球体であれば計算は簡単なのですが、残念ながらそうではなく、WGS 84 もそうですが一般に地球は楕円体でモデリングされるので、計算がややこしくなります。
PROJなどのライブラリが容易に利用できる場合はそれに任せてもいいかもしれません。一方、用途によってはPROJなどの巨大な依存を組み込まずに自分でさっさと実装したほうが良いでしょう。変換のコードはとても短いので。
導出の方法はいくつかあると思いますが、ここでは安直に考えていきます。
安直に導出してみる
まずサンプルコード (Python) を示します。解説はあとで。
# 基本の楕円体パラメータ
A = 6378137.0 # 長半径
INV_F = 298.257223563 # 偏平率の逆数 (exact)
# 楕円体パラメータから導ける値
F = 1 / INV_F # 偏平率
# B = A * (1 - F) # 短半径
E_SQ = F * (2 - F) # 離心率の2乗
# E = E_SQ ** 0.5 # 離心率
# 入力: WGS 84 Geographic 3D (EPSG:4979)
(lat_deg, lng_deg, height) = (34.290, 135.630, 100)
print(f"INPUT (WGS 84 Geographic): {lat_deg}, {lng_deg}, {height}")
(phi, lam) = (math.radians(lat_deg), math.radians(lng_deg))
# 下式は r = 1 / (1 / A ** 2 + (tan_psi / b) ** 2) ** 0.5 を整理したもの
r = A / (1 + (math.tan(phi) * math.tan(phi)) * (1 - E_SQ)) ** 0.5
# 下式は z = 1 / ((1 / (A * tan_psi) ** 2) + (1 / b ** 2)) ** 0.5 を整理したもの
tan_psi = (1 - E_SQ) * math.tan(phi)
z = A / (1 / (tan_psi * tan_psi) + 1 / (1 - E_SQ)) ** 0.5
x = r * math.cos(lam)
y = r * math.sin(lam)
# 高さを加味する
nx = math.cos(phi) * math.cos(lam)
ny = math.cos(phi) * math.sin(lam)
nz = math.sin(phi)
(x, y, z) = (x + nx * height, y + ny * height, z + nz * height)
print("OUTPUT (WGS 84 Geocentric):", x, y, z)
計算結果は、PROJ による結果と比べて数ナノメートル以下の差しかありませんでした。PROJの計算はさらに整理されているので（後述）、数値計算の誤差ということなのでしょう。
解説
Step 1. 地心緯度を求める
地心座標 (x, y, z) を計算する過程では、地理緯度 $\varphi$ ではなく地心緯度 $\psi$ を使います。
我々がふだん「緯度」というときに指すものは地理緯度です。地理緯度と地心座標の違いについては、Wikipediaなどを参照してください。
地心緯度は以下のようにして求められます：
$$
\tan \psi = (1-e^2) \tan \varphi
$$
Step 2. 楕円体上の点の計算
まずは高さを無視して、経緯度 (lat, lng) に対応する楕円体表面上の点 (x, y, z) を求めましょう。
（以下、図がないとつらいかも…？）
地心 $(0, 0, 0)$ から赤道面上の $(x, y, 0)$ までの距離を $r$ とすると、地心緯度 $\psi$ と、$z$、$r$ の関係は：
$$
\tan \psi = \frac{z}{r}
$$
$$
\frac{r^2}{a^2} + \frac{z^2}{b^2} = 1
$$
これらを簡単に整理すれば、緯度をもとに 𝑧 と 𝑟 が求められます：
$$
\begin{align}
z &= 1 / \big( \frac{1}{(a \tan \psi) ^ 2} + \frac{1}{b ^ 2} \big) \\
r &= 1 / \big( \frac{1}{a ^ 2} + (\frac{\tan \psi}{b})^2 \big)
\end{align}
$$
x と y は以下のように求まります。 $\theta$ は経度です。
$$
\begin{align}x &= r \cdot \cos(\theta) \\
y &= r \cdot \sin(\theta)\end{align}
$$
Step 3. 高さを加算する
ここまでで楕円体表面上の点が求まりましたが、楕円体表面からの高さ（楕円体高）height が考慮されていません。
楕円体高を反映するには、当該の点での楕円体の法線ベクトルを使って高さを加味するだけでOKです。楕円体上の法線ベクトルは「地理緯度」から簡単に求まります。
サンプルコードの最後でこの処理をしています。
計算をさらに整理する
さきほどの計算（コード）を、三角関数の性質などを利用してさらに整理すると以下のようになります：
A = 6378137.0 # 長半径
INV_F = 298.257223563 # 偏平率の逆数 (exact)
# 楕円体パラメータから導ける値
F = 1 / INV_F # 偏平率
# B = A * (1 - F) # 短半径
E_SQ = F * (2 - F) # 離心率の2乗
# E = E_SQ ** 0.5 # 離心率
# 入力: WGS 84 Geographic 3D (EPSG:4979)
(lat_deg, lng_deg, height) = (34.290, 135.630, 100)
print(f"INPUT (WGS 84 Geographic): {lat_deg}, {lng_deg}, {height}")
(phi, lam) = (math.radians(lat_deg), math.radians(lng_deg))
n = A / (1 - E_SQ * math.sin(phi) * math.sin(phi)) ** 0.5
x = (n + height) * math.cos(phi) * math.cos(lam)
y = (n + height) * math.cos(phi) * math.sin(lam)
z = (n * (1 - E_SQ) + height) * math.sin(phi)
print("OUTPUT (WGS 84 Geocentric):", x, y, z)
ここまで整理すると、PROJの中で使われている計算と同じになります。https://github.com/OSGeo/PROJ/blob/master/src/conversions/cart.cpp
※実は、幾何学的にうまく考えると、いきなりこの式を導くことができるようです（Wikipediaの記事を参照）。
おわり
逆変換はとても難しい
逆方向の (x, y, z) から (lng, lat, height) への変換も考えてみたいところですが、そちらは難解なだけでなく、解の安定性などにもややこしさがあり、色々な研究がなされているトピックだそうです。
最近では、Vermeille の手法というものがよく使われているように見えます。
楕円体高 vs. 標高
多くの3Dデータ（例えば JGD2011 などの日本の座標参照系が使われているデータ）では、高さが「楕円体高」ではなく「標高」で表わされていることがあります。それらの座標系から変換する場合は、標高と楕円体高を相互に変換する処理も必要になります。
具体的には、国土地理院さんが公開しているようなジオイドモデルから得られたジオイド高を楕円体高に加算する必要があります。楕円体で表現された「理想」の地球モデルと、「現実」の地球の等ポテンシャル面（仮想平均海面）とを変換する必要があるということでしょうか。
今回の計算はWGS84同士での変換ですので、ジオイドについては考慮する必要がなく、考慮していません。