引言
如果需要計算蜃景，則需計算不同折射率下的光路軌跡，這裡考慮層狀線性變化的折射率對光路的影響。
模型
在水平面建立坐標系。假設折射率分層變化$\eta(y) = \eta_0 + ky$，和水平距離的尺度而言，顯然梯度$k$是一個小量，這也是線性近似處理的原因。考慮從地面$(x_1, 0)$射出的光線最終回到地面$(x_2, 0)$，下面計算其軌跡。
計算
考慮作用量
$$
S = \int_{x_1}^{x_2} \mathcal{L}(y, y’) dx
$$
其中
$$
\mathcal{L}(y, y’) = \eta(y)\sqrt{1 + (y’)^2}
$$
考慮$\delta S = 0$，得到歐拉方程
$$
\partial_y \mathcal{L} - \frac{d}{dx} (\partial_{y’} \mathcal{L}) = 0
$$
即
$$
k \sqrt{1 + (y’)^2} = \frac{d}{dx} \left[ \frac{(n_0 + ky)y’}{\sqrt{1 + (y’)^2}} \right]
$$
這個方程不易解，不過有微擾的辦法。
微擾
考慮$k = 0$的情況，這樣得到零階的解
$$
y_0(x) = \pm \frac{c}{\sqrt{n_0^2 - c^2}} x + d = mx + d
$$
其中$m$，$c$和$d$均為常數。此時光路為直線，這也與實際物理意義一致。假設最終解可以寫成級數展開的形式
$$
y(x) = y_0(x) + k y_1(x) + k^2 y_2(x) + \cdots
$$
現在考慮方程的左邊，考慮近似
$$
(y’)^2 = (y_0’ + k y_1’ + \cdots)^2 \approx (y_0’)^2 + 2 k y_0’ y_1’ + O(k^2)
$$
這樣泰勒展開得到
$$
\text{LHS} = k \sqrt{1 + (y’)^2} \approx k \sqrt{1 + (y_0’)^2} + k^2 \frac{y_0’ y_1’}{\sqrt{1 + (y_0’)^2}} + O(k^3)
$$
接著考慮右邊函數的分子
$$
\begin{equation}
\begin{gathered}
& (n_0 + k y_0 + k^2 y_1 + \cdots)(y_0’ + k y_1’ + \cdots) \\
\approx& n_0 y_0’ + k (n y_1’ + y_0 y_0’) + k^2 (n_0 y_2’ + y_0 y_1’ + y_1 y_0’) + O(k^3)
\end{gathered}
\end{equation}
$$
以及分母
$$
\frac{1}{\sqrt{1 + (y’)^2}} \approx \frac{1}{\sqrt{1 + (y_0’)^2}} - \frac{k y_0’ y_1’}{(1 + (y_0’)^2)^{3/2}} + O(k^2)
$$
因此
$$
\begin{equation}
\begin{gathered}
\frac{(n_0 + ky)y’}{\sqrt{1 + (y’)^2}} \\
\approx \left[ n_0 y_0’ + k (n y_1’ + y_0 y_0’) + \cdots \right] \left[ \frac{1}{\sqrt{1 + (y_0’)^2}} - \frac{k y_0’ y_1’}{(1 + (y_0’)^2)^{3/2}} + \cdots \right]
\end{gathered}
\end{equation}
$$
比較系數
右邊函數的零階為
$$
\frac{n y_0’}{\sqrt{1 + (y_0’)^2}}
$$
一階為
$$
\frac{n y_1’ + y_0 y_0’}{\sqrt{1 + (y_0’)^2}} - \frac{n (y_0’)^2 y_1’}{(1 + (y_0’)^2)^{3/2}}
$$
現在計算右邊微分
$$
\frac{d}{dx} \left[ \frac{n_0 y_0’}{\sqrt{1 + (y_0’)^2}} + k \left( \frac{n_0 y_1’ + y_0 y_0’}{\sqrt{1 + (y_0’)^2}} - \frac{n_0 (y_0’)^2 y_1’}{(1 + (y_0’)^2)^{3/2}} \right) + O(k^2) \right]
$$
其中$y_0’ = m$，令$a = \sqrt{1 + m^2}$，因此右邊微分號裡面的一階系數為
$$
\begin{equation}
\begin{aligned}
& y_1’ \left( \frac{n_0}{a} - \frac{n_0 m^2}{a^3} \right) + \frac{y_0 m}{a} \\
=& y_1’ \cdot \frac{n_0 (1 + m^2 - m^2)}{(1 + m^2)^{3/2}} + \frac{y_0 m}{\sqrt{1 + m^2}} \\
=& \frac{n_0 y_1’}{(1 + m^2)^{3/2}} + \frac{y_0 m}{\sqrt{1 + m^2}}
\end{aligned}
\end{equation}
$$
因此
$$
\frac{d}{dx} \left[ \frac{n_0 y_1’}{(1 + m^2)^{3/2}} + \frac{y_0 m}{\sqrt{1 + m^2}} \right] = \frac{n_0 y_1’’}{(1 + m^2)^{3/2}} + \frac{m^2}{\sqrt{1 + m^2}}
$$
一階微擾方程
比較左右兩邊的一階系數得到方程
$$
\frac{n_0 y_1’’}{(1 + m^2)^{3/2}} + \frac{m^2}{\sqrt{1 + m^2}} = \sqrt{1 + m^2}
$$
這是一個簡單的二階微分方程，就不作展開了。最後通解為
$$
y(x) \approx y_0 + k y_1 = (m x + d) + k \left( \frac{1 + m^2}{2 n_0} x^2 + c_1 x + c_2 \right)
$$
其中$c_1$與$c_2$也均為常數，由初始條件決定。根據文中條件可知，軌跡近似為
$$
y(x) \approx k\frac{1 + m^2}{2n_0}(x - x_1)(x - x_2)
$$
過兩點之拋物線，須知$m$即為初始斜率，而$x_2$由$x_1$和$m$決定
$$
x_2 \approx x_1 - \frac{2mn_0}{k(1 + m^2)}
$$
注意這裡的初始斜率是一個小量。
討論
有兩點結論：
- 地表光線經過大氣層折射回到地面的軌跡是扁平的拋物線。
- 梯度必須為負，即隨高度上升折射率下降。
這個條件可以發生在傍晚的天空：地表由於比熱較小的緣故迅速散熱而導致溫度下降，同時由於上層空氣溫度依舊較高，從而讓遠處光線折射回到地面成為可能。同樣討論也適用於海面，海水由於比熱較大的緣故導致接近海面的氣溫不高。後續會討論出現海市蜃樓的條件。
我是在返程的火車上思考這個問題，當時我看著窗外的夕陽在紙上寫下上述結果。
补充
這裡以拋物線近似估計仰角。簡單起見，考慮$m \approx 0$以及$n_0 \approx 1$，假設光源位於原點。
首先在不考慮偏折的情況下，距離原點$L$的豎直方向長度為$\Delta l$的瞳孔收集到的臨界光束為
$$
y_1 = 0, \quad y_2 = \frac{\Delta l}{L}x
$$
相應臨界角為$\theta_1 = 0$和$\theta_2 = \Delta l / L$，因而角度差為$\Delta \theta = \Delta l / L$。
接下來考慮有偏折的情況，臨界光束分別為
$$
y_1 = \frac{k}{2}x(x-x_1), \quad y_2 = \frac{k}{2}x(x-x_2)
$$
這裡$x_2 > x_1 = L$。臨界角為$\theta_1 = -kx_1 / 2$和$\theta_2 = -kx_2 / 2$，角度差為$\Delta \theta = -k(x_2 - x_1) / 2$，同時注意到$\Delta l = -\frac{k}{2}x_1(x_2 - x_1)$，因此依舊有$\Delta \theta = \Delta l / L$，角度差完全一致。假設光線密度正比於角度差，那麼光線偏折並不會讓進入瞳孔的光線變多，而是和無偏折的情況一致。
兩者情況的差異在於進入瞳孔的光線的角度，後者可以繞開原先水平方向阻擋的障礙物。這個仰角即為$-kL / 2$，這可以作為一個測量梯度$k$的辦法。
下面計算像的位置，進入瞳孔的臨界光線的切線方程分別為
$$
y = \frac{kx_1}{2}(x - x_1), \quad y - \Delta l = \frac{k}{2}(2x_1 - x_2)(x - x_1)
$$
解出位置
$$
x = x_1 + \frac{2\Delta l}{k(x_2 - x_1)} = L - \frac{\Delta l}{\Delta \theta} = 0
$$
位置不變。以及高度
$$
y = -\frac{kx_1^2}{2} = -\frac{kL^2}{2}
$$
可見偏折等效於提升光源的高度。
倍率
現在考慮放大倍率，假設物體有一定的尺寸不能忽略，假設高度為$h$，把頂端放在原點處，類似有臨界光線
$$
y_1 = \frac{k}{2}x(x-x_1), \quad y_2 = \frac{k}{2}x(x-x_2)
$$
注意此時$x_1 \neq L$。代入數據得到
$$
-h = \frac{k}{2}L(L - x_1), \quad -h + \Delta l = \frac{k}{2}L(L - x_2)
$$
這樣$\Delta l = -\frac{k}{2}L(x_2 - x_1)$，依舊有$\Delta \theta = \Delta l / L$。現在計算像的位置，類似兩條臨界光線的切線方程分別為
$$
y + h = \frac{k}{2}(2L - x_1)(x - L), \quad y + h - \Delta l = \frac{k}{2}(2L - x_2)(x - L)
$$
類似可以解出位置
$$
x = L + \frac{2\Delta l}{k(x_2 - x_1)} = L - \frac{\Delta l}{\Delta \theta} = 0
$$
維持不變，而求解高度則需要解出未知參數$x_1$
$$
x_1 = L + \frac{2h}{kL}
$$
代入得到高度
$$
y = -h - \frac{k}{2}L(2L - x_1) = -h - \frac{kL^2}{2} + h = - \frac{kL^2}{2}
$$
表明物體整體被提高了$-kL^2 / 2$，而大小不變。