用PythonMatplotlib复现偏振光实验从马吕斯定律到波片可视化附完整代码偏振光是光学中一个既基础又重要的概念它不仅揭示了光的横波本质还在3D显示、光学测量、量子通信等领域有着广泛应用。传统物理实验往往受限于设备条件和环境因素而借助Python的数据可视化能力我们可以突破这些限制在代码中构建一个理想实验室。本文将带你用NumPy和Matplotlib从零开始完整复现偏振光实验的核心过程。1. 环境准备与基础理论在开始编码前我们需要配置合适的Python环境并理解关键的光学原理。推荐使用Anaconda创建独立环境conda create -n polarization python3.9 conda activate polarization pip install numpy matplotlib ipython偏振光的数学表示是模拟的基础。线偏振光可以表示为E_x E_0 * cos(ωt - kz) * cosθ E_y E_0 * cos(ωt - kz) * sinθ其中θ是偏振方向与x轴的夹角。当这束光通过检偏器另一个偏振片时根据马吕斯定律输出光强I与入射光强I₀的关系为I I₀ * cos²(φ)φ是两个偏振片透振方向的夹角。波片的工作原理则基于双折射效应。λ/4波片会引入π/2的相位差将线偏振光转换为椭圆偏振光λ/2波片引入π相位差使偏振方向旋转2θ。这些效应可以通过琼斯矩阵严格描述光学元件琼斯矩阵线偏振片(θ0)[[1,0],[0,0]]λ/4波片(快轴x)exp(iπ/4)[[1,0],[0,-i]]λ/2波片(快轴x)i[[1,0],[0,-1]]2. 马吕斯定律的可视化验证让我们首先实现马吕斯定律的验证。关键步骤包括生成偏振光、模拟检偏过程以及绘制光强随角度变化的曲线。import numpy as np import matplotlib.pyplot as plt def malus_law(I0, phi): 计算马吕斯定律预测的光强 return I0 * (np.cos(np.radians(phi)) ** 2) # 参数设置 I0 1.0 # 初始光强 angles np.linspace(0, 360, 361) # 0到360度 intensities malus_law(I0, angles) # 绘制结果 plt.figure(figsize(10,6)) plt.plot(angles, intensities, label理论曲线) plt.xlabel(检偏器旋转角度(度), fontsize12) plt.ylabel(相对光强, fontsize12) plt.title(马吕斯定律验证, fontsize15) plt.grid(True, linestyle--, alpha0.7) plt.legend() plt.show()运行这段代码会得到一条完美的余弦平方曲线。为了增加实验真实感我们可以添加一些实验噪声# 添加高斯噪声模拟真实实验 noise np.random.normal(0, 0.02, len(angles)) experimental_intensities intensities noise plt.figure(figsize(10,6)) plt.scatter(angles[::10], experimental_intensities[::10], colorred, label模拟实验数据) plt.plot(angles, intensities, label理论曲线) plt.legend() plt.show()数据拟合可以进一步验证马吕斯定律。使用SciPy的curve_fit函数from scipy.optimize import curve_fit def fit_func(phi, I0): return I0 * (np.cos(np.radians(phi)) ** 2) popt, pcov curve_fit(fit_func, angles, experimental_intensities) print(f拟合得到I0{popt[0]:.3f} (理论值1.0))3. 波片效应的动态演示波片对偏振光的改变是偏振光学中最迷人的现象之一。我们将创建交互式可视化来展示λ/4和λ/2波片的效果。3.1 λ/4波片产生椭圆偏振光λ/4波片的关键在于引入π/2的相位差。我们可以用以下函数模拟这个过程def quarter_wave_plate(E_in, theta): 模拟λ/4波片效应 E_in: 入射光电矢量 [Ex, Ey] theta: 波片快轴与x轴夹角(度) theta_rad np.radians(theta) # 琼斯矩阵 J np.exp(1j*np.pi/4) * np.array([ [np.cos(theta_rad)**2 1j*np.sin(theta_rad)**2, (1-1j)*np.sin(theta_rad)*np.cos(theta_rad)], [(1-1j)*np.sin(theta_rad)*np.cos(theta_rad), np.sin(theta_rad)**2 1j*np.cos(theta_rad)**2] ]) return J E_in为了直观展示效果我们创建偏振态随时间变化的动画from matplotlib.animation import FuncAnimation def create_polarization_animation(waveplate_typequarter): fig plt.figure(figsize(8,8)) ax fig.add_subplot(111, projection3d) # 初始化线偏振光 (沿x轴) E_in np.array([1, 0]) # 设置波片角度 wp_angle 45 # 波片旋转角度 # 生成时间序列 t np.linspace(0, 2*np.pi, 100) def update(frame): ax.cla() wp_angle frame % 180 if waveplate_type quarter: E_out quarter_wave_plate(E_in, wp_angle) else: E_out half_wave_plate(E_in, wp_angle) # 计算电场随时间变化 Ex np.real(E_out[0] * np.exp(1j*t)) Ey np.real(E_out[1] * np.exp(1j*t)) # 绘制偏振态 ax.plot(Ex, Ey, t, lw2) ax.set_xlim(-1,1) ax.set_ylim(-1,1) ax.set_zlim(0,2*np.pi) ax.set_title(f{waveplate_type} waveplate at {wp_angle}°) anim FuncAnimation(fig, update, frames180, interval50) plt.close() return anim3.2 λ/2波片旋转偏振方向λ/2波片会使偏振方向旋转两倍于波片旋转角度。其琼斯矩阵表示为def half_wave_plate(E_in, theta): 模拟λ/2波片效应 theta_rad np.radians(theta) J 1j * np.array([ [np.cos(2*theta_rad), np.sin(2*theta_rad)], [np.sin(2*theta_rad), -np.cos(2*theta_rad)] ]) return J E_in我们可以创建一个可视化函数来对比输入和输出偏振态def visualize_waveplate_effect(waveplate_typehalf, angle30): t np.linspace(0, 2*np.pi, 100) E_in np.array([1, 0]) # 沿x轴的线偏振光 if waveplate_type half: E_out half_wave_plate(E_in, angle) else: E_out quarter_wave_plate(E_in, angle) # 计算电场 Ex_in np.real(E_in[0] * np.exp(1j*t)) Ey_in np.real(E_in[1] * np.exp(1j*t)) Ex_out np.real(E_out[0] * np.exp(1j*t)) Ey_out np.real(E_out[1] * np.exp(1j*t)) # 绘图 fig, (ax1, ax2) plt.subplots(1, 2, figsize(12,5)) ax1.plot(Ex_in, Ey_in) ax1.set_title(Input Polarization) ax1.set_xlim(-1.2,1.2) ax1.set_ylim(-1.2,1.2) ax2.plot(Ex_out, Ey_out) ax2.set_title(fAfter {waveplate_type} waveplate at {angle}°) ax2.set_xlim(-1.2,1.2) ax2.set_ylim(-1.2,1.2) plt.tight_layout() plt.show()4. 综合实验平台构建将上述模块整合成一个交互式偏振光实验平台我们可以使用IPython的交互式控件from ipywidgets import interact, FloatSlider, Dropdown def polarization_experiment(waveplate_typequarter, waveplate_angle0, polarizer_angle0): # 生成线偏振光 E_in np.array([1, 0]) # 波片效应 if waveplate_type quarter: E_after_wp quarter_wave_plate(E_in, waveplate_angle) elif waveplate_type half: E_after_wp half_wave_plate(E_in, waveplate_angle) else: E_after_wp E_in # 检偏器效应 phi np.radians(polarizer_angle) J_pol np.array([ [np.cos(phi)**2, np.cos(phi)*np.sin(phi)], [np.cos(phi)*np.sin(phi), np.sin(phi)**2] ]) E_final J_pol E_after_wp # 计算光强 I np.sum(np.abs(E_final)**2) # 可视化 t np.linspace(0, 2*np.pi, 100) Ex np.real(E_after_wp[0] * np.exp(1j*t)) Ey np.real(E_after_wp[1] * np.exp(1j*t)) plt.figure(figsize(8,8)) plt.plot(Ex, Ey, lw2) plt.title(fPolarization State\nWaveplate: {waveplate_type} at {waveplate_angle}°\n fPolarizer: {polarizer_angle}°\nIntensity: {I:.3f}) plt.xlim(-1.2,1.2) plt.ylim(-1.2,1.2) plt.grid(True) plt.show() # 创建交互界面 interact(polarization_experiment, waveplate_typeDropdown(options[none, quarter, half]), waveplate_angleFloatSlider(min0, max180, step5, value0), polarizer_angleFloatSlider(min0, max180, step5, value0))这个交互式工具允许你选择波片类型无、λ/4、λ/2调整波片角度0-180度旋转检偏器角度实时观察输出偏振态和光强5. 进阶应用与扩展掌握了基础模拟后我们可以探索更复杂的偏振光学现象5.1 椭圆偏振光的表征椭圆偏振光可以用斯托克斯参数完整描述def stokes_parameters(E): 计算斯托克斯参数 Ex, Ey E S0 np.abs(Ex)**2 np.abs(Ey)**2 S1 np.abs(Ex)**2 - np.abs(Ey)**2 S2 2 * np.real(Ex * np.conj(Ey)) S3 2 * np.imag(Ex * np.conj(Ey)) return np.array([S0, S1, S2, S3]) # 示例计算不同偏振态的斯托克斯参数 linear_x np.array([1, 0]) # 沿x轴线偏振 linear_y np.array([0, 1]) # 沿y轴线偏振 circular_r np.array([1, -1j])/np.sqrt(2) # 右旋圆偏振 print(Linear X:, stokes_parameters(linear_x)) print(Linear Y:, stokes_parameters(linear_y)) print(Circular R:, stokes_parameters(circular_r))5.2 偏振光的干涉当两束偏振光相遇时它们的干涉模式取决于偏振状态def polarization_interference(E1, E2, phase_diff): 模拟两束偏振光的干涉 E_total E1 E2 * np.exp(1j*phase_diff) return np.sum(np.abs(E_total)**2) # 创建干涉图样 phase_diffs np.linspace(0, 2*np.pi, 100) intensities [polarization_interference( np.array([1,0]), # 沿x轴线偏振 np.array([0,1]), # 沿y轴线偏振 pd) for pd in phase_diffs] plt.plot(phase_diffs, intensities) plt.xlabel(Phase Difference (rad)) plt.ylabel(Interference Intensity) plt.title(Interference of Orthogonal Polarizations) plt.show()5.3 偏振成像应用偏振信息在计算机视觉中有重要应用。我们可以模拟一个简单的偏振成像场景def simulate_polarization_image(size256): 模拟包含不同偏振特性的场景 image np.zeros((size, size, 3)) # 区域1漫反射表面部分偏振 image[50:150, 50:150] [0.8, 0.3, 0.3] # 区域2金属反射强偏振 image[50:150, 150:250] [0.2, 0.7, 0.7] # 区域3玻璃透射改变偏振态 image[150:250, 50:150] [0.5, 0.5, 0.9] # 区域4各向同性散射非偏振 image[150:250, 150:250] [0.9, 0.6, 0.2] return image def capture_with_polarizer(image, angle): 模拟通过偏振片拍摄 theta np.radians(angle) intensity (image[:,:,0]*np.cos(theta)**2 image[:,:,1]*np.sin(theta)**2 image[:,:,2]*np.sin(2*theta)/2) return intensity # 模拟不同偏振角度拍摄 angles [0, 45, 90, 135] images [capture_with_polarizer(simulate_polarization_image(), a) for a in angles] # 显示结果 fig, axes plt.subplots(2, 2, figsize(10,10)) for ax, img, angle in zip(axes.flat, images, angles): ax.imshow(img, cmapgray) ax.set_title(fPolarizer at {angle}°) ax.axis(off) plt.tight_layout() plt.show()