Python光学仿真：基于菲涅耳公式的波动模型实现

苏见山

发布时间：2026-10-14| 4人浏览过

原创

本文详解如何使用Python进行光学仿真，通过引入相位项构建波动模型。重点解析菲涅耳公式在计算光波透过率中的应用，区分p分量与s分量。提供完整Python代码示例，绘制不同入射角下的振幅反射率与透过率变化曲线，实现从几何光学到波动光学的仿真升级。

从物理机制来看，波动模型比光线模型更贴近光的本质；但从物理学发展脉络看，波动光学旨在解决几何光学无法处理的问题，可视作光线模型的进阶版本。在编程实现层面，波动光学可简化理解为在光线模型基础上增加一个相位项。

波动模型

通常情况下，频率、振幅和相位三个特征足以确定空间中的波场，因此光波场可表示为：


import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
z = np.arange(15,200)*10    #单位为nm
x = np.arange(15,200)*10
x,z = np.meshgrid(x,z)      #创建坐标系
E = 1/np.sqrt(x**2+z**2)*np.cos(2*np.pi*np.sqrt(x**2+z**2)/(532*1e-9))
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_surface(x,z,E)
plt.show()

其结果如图所示

菲涅耳公式

几何光学依据费马原理可推导折射定律，但无法获取光波的透过率。菲涅耳公式在几何光学基础上解决了这一缺失。

由于光是横波集合，可根据电矢量振动方向将其分解为平行于入射面（p分量）和垂直于入射面（s分量）的两个分量。当光束在折射率为n1和n2的两种介质界面发生折射时，对于p光，其电矢量平行于入射面，磁矢量垂直于入射面（仅含s分量）；对于s光则相反，具体如图所示。

因此，对于p光而言，其关系式为：

针对磁矢量，存在如下关系：

我们可以利用Python绘制不同入射角下振幅反射率和透过率的变化情况。


import matplotlib.pyplot as plt
import numpy as np
def fresnel(theta, n1, n2):
    theta = theta*np.pi/180
    xTheta = np.cos(theta)
    mid = np.sqrt(1-(n1/n2*np.sin(theta))**2)          #中间变量
    rp = (n2*xTheta-n1*mid)/(n2*xTheta+n1*mid)  #p分量振幅反射率
    rs = (n1*xTheta-n2*mid)/(n1*xTheta+n2*mid)
    tp = 2*n1*xTheta/(n2*xTheta+n1*mid)
    ts = 2*n1*xTheta/(n1*xTheta+n2*mid)
    return rp, rs, tp, ts
def testFres(n1=1,n2=1.45):         #默认n2为1.45
    theta = np.arange(0,90,0.1)+0j
    a = theta*np.pi/180
    rp,rs,tp,ts = fresnel(theta,n1,n2)
    fig = plt.figure(1)
    plt.subplot(1,2,1)
    plt.plot(theta,rp,'-',label='rp')
    plt.plot(theta,rs,'-.',label='rs')
    plt.plot(theta,np.abs(rp),'--',label='|rp|')
    plt.plot(theta,np.abs(rs),':',label='|rs|')
    plt.legend()
    plt.subplot(1,2,2)
    plt.plot(theta,tp,'-',label='tp')
    plt.plot(theta,ts,'-.',label='ts')
    plt.plot(theta,np.abs(tp),'--',label='|tp|')
    plt.plot(theta,np.abs(ts),':',label='|ts|')
    plt.legend()
    plt.show()
if __init__=="__main__":
    testFres()

所得图像如下所示：

通过Python绘图，只需将上述程序中的testFres替换为以下代码即可实现。


def testFres(n1=1,n2=1.45):
    theta = np.arange(0,90,0.1)+0j
    a = theta*np.pi/180
    rp,rs,tp,ts = fml.fresnel(theta,n1,n2)
    Rp = np.abs(rp)**2
    Rs = np.abs(rs)**2
    Rn = (Rp+Rs)/2
    Tp = n2*np.sqrt(1-(n1/n2*np.sin(a))**2)/(n1*np.cos(a))*np.abs(tp)**2
    Ts = n2*np.sqrt(1-(n1/n2*np.sin(a))**2)/(n1*np.cos(a))*np.abs(ts)**2
    Tn = (Tp+Ts)/2
    fig = plt.figure(2)
    plt.subplot(1,2,1)
    plt.plot(theta,Rp,'-',label='R_p')
    plt.plot(theta,Rs,'-.',label='R_s')
    plt.plot(theta,Rn,'-',label='R_n')
    plt.legend()
    plt.subplot(1,2,2)
    plt.plot(theta,Tp,'-',label='T_p')
    plt.plot(theta,Ts,'-.',label='T_s')
    plt.plot(theta,Tn,'--',label='T_n')
    plt.legend()
    plt.show()

最终得出：