[Python知识库]【Python光学仿真】光线在线性分层的变折射率介质中传播的规律

光线模型默许了一条基本的原则：即道路千万条，光只走最短的那条。

为了让问题充分简单，先考虑折射率分层分布的情况，即假设空间中折射率分布为

$$n=n_0+kx$$

取光线传播方向为坐标平面，建立坐标系$xOy$，记光束入射设介质折射率分布为$n(x,y)$，光束与$x$轴夹角为$\theta$。

考虑折射率只在$x$轴方向变化，即$n(x)$。设介质中某点$(x,y)$入射光线$(x,y,\theta )$，其临近的出射位置为$(x+\text{d}x,y+\text{d}y)$，有

$$n(x)\sin \theta (x)=n(x+\text{d}x)\sin \theta (x+\text{d}x)$$

记外界折射率为$n_0$，入射角为${\theta }_0$，则对于任意一点，有

$$n(x)\sin \theta (x)=n_0\sin {\theta }_0$$

则

$$\begin{array}{rl}\text{d}y& =\text{d}x\tan \theta (x)\\ & =\text{d}x\frac{\sin \theta (x)}{\sqrt{1?{\sin }^2\theta (x)}}\\ (\frac{\text{d}y}{\text{d}x}{)}^2& =\frac{n_0^2{\sin }^2{\theta }_0}{n(x{)}^2?n_0^2{\sin }^2{\theta }_0}\end{array}$$

此即在折射率分层介质中，光的传播规律。

当$n(x)=n_0+kx$时，记$N_0=n_0\sin {\theta }_0$，上式写为

$$\begin{array}{rl}(\frac{\text{d}y}{\text{d}x}{)}^2& =\frac{N_0^2}{(n_0+kx{)}^2?N_0^2}\\ \frac{\text{d}y}{\text{d}x}& =\frac{1}{\sqrt{(\frac{n_0+kx}{N_0}{)}^2?1}}\end{array}$$

令$X=\frac{n_0+kx}{N_0}$，即$x=\frac{N_{0}X?n_0}{k}$，上式变为

$$\begin{array}{rl}\frac{\text{d}y}{\text{d}X}\frac{k}{2N_0}& =\frac{1}{\sqrt{X^2?1}}\\ \to y& =\frac{2N_0}{k}\int \frac{1}{\sqrt{X^2?1}}\text{d}x\\ \stackrel{X=\frac{1}{\cos \alpha }}{}& =\frac{2N_0}{k}\int \frac{\text{d}\alpha }{\cos \alpha }\\ & =\frac{2N_0}{k}\ln |\frac{1+\sin \alpha }{\cos \alpha }|\\ \stackrel{X=\frac{1}{\cos \alpha }}{}& =\frac{2N_0}{k}\ln |X+\sqrt{X^2?1})|\end{array}$$

解得

$$y=\frac{2N_0}{k}\ln |\frac{n_0+kx}{N_0}+\sqrt{(\frac{n_0+kx}{N_0}{)}^2?1})|+C$$

设入射位置为$(0,0)$，则上式变为

$$y=\frac{2N_0}{k}\ln \frac{n_0+kx+\sqrt{(n_0+kx{)}^2?N_0^2}}{n_0+\sqrt{n_0^2?N_0^2}}$$

代码为

def linearRefractive(n0=1,k=0.01,theta0=0.5):
    x = np.arange(400)
    nx = n0 + k*x
    N0 = n0*np.sin(theta0)
    y = 2*N0/k*np.log(
        (nx+np.sqrt(nx**2-N0**2))/(n0+np.sqrt(n0**2-N0**2)))
    Y = 1/np.tan(theta0)*x+y[0]

    #折射率场
    N,_ = np.meshgrid(x,np.arange(np.max(Y)))
    N = n0 + k*N

    plt.imshow(N.T)
    plt.plot(y,x)
    plt.plot(Y,x)

    plt.show()

得到光线的传播图像为

图中两条线分别表示光线轨迹和光线的入射方向。

所谓分层介质，当然不必处处折射率梯度恒定，即$n(x)$可以选取更加复杂的变化规律，而此时下列方程不一定存在解析解

$$(\frac{\text{d}y}{\text{d}x}{)}^2=\frac{n_0^2{\sin }^2{\theta }_0}{n(x{)}^2?n_0^2{\sin }^2{\theta }_0}$$

故选取其差分方式

$$\text{d}y=\frac{n_0\sin {\theta }_0}{\sqrt{n(x{)}^2?n_0^2{\sin }^2{\theta }_0}}\text{d}x$$

当折射率周期变化时，有

def triangleRefractive(n0=1,k=0.5,theta0=0.5):
    x = np.arange(400)
    nx = n0 + 1 + k*np.sin(x/20)
    N0 = n0*np.sin(theta0)
    y = x*1.0
    for i in x[1:]:
        y[i] = N0/np.sqrt(nx[i]**2-N0**2)+y[i-1]
    Y = 1/np.tan(theta0)*x+y[0]

    #折射率场
    N,_ = np.meshgrid(x,np.arange(np.max(y)))
    N = n0 + k*np.sin(N/20)

    plt.imshow(N)
    plt.plot(x,y)
    plt.plot(x[Y<np.max(y)],Y[Y<np.max(y)])

    plt.show()

如图所示