用Python动态模拟偏振光实验：从代码理解马吕斯定律

偏振光现象在光学实验中常常让人感到抽象难懂——那些旋转的偏振片、相位延迟的波片，以及随之变化的光强曲线，在传统实验报告中往往只能通过静态图表和公式呈现。但如果你能用代码让这些物理过程"活"起来呢？本文将带你用Python的Matplotlib库构建一个 交互式偏振光模拟器 ，通过动态可视化深入理解马吕斯定律和波片工作原理。

1. 环境准备与基础概念

在开始编码前，我们需要配置合适的Python环境。推荐使用Anaconda创建专属虚拟环境：

conda create -n polarization python=3.9
conda activate polarization
pip install numpy matplotlib ipywidgets

偏振光模拟的核心是理解几个关键物理量：

• 线偏振光 ：电矢量仅沿单一方向振动的光波
• 马吕斯定律 ：描述偏振光通过检偏器后光强变化的规律，数学表达式为 $I = I_0 \cos^2\theta$
• 波片 ：通过双折射效应产生相位延迟的光学元件，常见的有λ/4波片（产生π/2相位差）和λ/2波片（产生π相位差）

下面这个对比表能帮助快速理解不同偏振元件的作用：

元件类型	作用效果	相位改变	输出偏振态
偏振片	过滤特定方向振动	无	线偏振光
λ/4波片	产生π/2相位差	90°	椭圆/圆偏振光
λ/2波片	产生π相位差	180°	线偏振光（方向旋转）

提示：在模拟中，我们将把光波表示为电矢量的二维投影，x和y分量分别对应快轴和慢轴方向的振动。

2. 构建偏振光模拟框架

让我们先建立一个基础的光波模型。线偏振光可以表示为沿特定方向振动的电磁波：

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

def linear_polarized_wave(amplitude, angle, time):
    """生成线偏振光波"""
    Ex = amplitude * np.cos(angle) * np.cos(time)
    Ey = amplitude * np.sin(angle) * np.cos(time)
    return Ex, Ey

接下来实现偏振片的数学模型。根据马吕斯定律，偏振片会衰减与透光轴垂直的振动分量：

def polarizer(input_Ex, input_Ey, polarizer_angle):
    """模拟偏振片效应"""
    # 将输入光分解到偏振片的透光轴和消光轴方向
    E_parallel = input_Ex * np.cos(polarizer_angle) + input_Ey * np.sin(polarizer_angle)
    E_perpendicular = -input_Ex * np.sin(polarizer_angle) + input_Ey * np.cos(polarizer_angle)
    
    # 消光轴分量被完全阻挡
    output_Ex = E_parallel * np.cos(polarizer_angle)
    output_Ey = E_parallel * np.sin(polarizer_angle)
    return output_Ex, output_Ey

为了直观展示偏振片的作用，我们可以创建一个动态演示：

def animate_polarizer_rotation():
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
    time = np.linspace(0, 2*np.pi, 100)
    
    def update(frame):
        ax1.clear(); ax2.clear()
        angle = frame * np.pi / 180  # 转换为弧度
        
        # 原始光波（水平偏振）
        Ex, Ey = linear_polarized_wave(1, 0, time)
        ax1.plot(time, Ex, 'r', label='Ex')
        ax1.plot(time, Ey, 'b', label='Ey')
        
        # 通过旋转偏振片后的光波
        Ex_p, Ey_p = polarizer(Ex, Ey, angle)
        ax2.plot(time, Ex_p, 'r--', label='Ex after polarizer')
        ax2.plot(time, Ey_p, 'b--', label='Ey after polarizer')
        
        # 设置图表属性
        for ax in [ax1, ax2]:
            ax.set_ylim(-1.1, 1.1)
            ax.legend()
        ax1.set_title('Original Polarized Light')
        ax2.set_title(f'After Polarizer (θ={frame}°)')
    
    anim = FuncAnimation(fig, update, frames=np.arange(0, 180, 2), interval=100)
    plt.close()
    return anim

这段代码会生成一个动画，展示当偏振片旋转时，透射光波各分量的变化情况。你可以明显看到当偏振片透光轴与偏振方向垂直时，光强变为零。

3. 波片效应的动态模拟

波片的模拟比偏振片更复杂，因为它会引入相位延迟。我们需要分别处理快轴和慢轴分量：

def wave_plate(input_Ex, input_Ey, plate_angle, phase_retardation):
    """模拟波片效应"""
    # 将输入光分解到波片的快轴和慢轴方向
    E_fast = input_Ex * np.cos(plate_angle) + input_Ey * np.sin(plate_angle)
    E_slow = -input_Ex * np.sin(plate_angle) + input_Ey * np.cos(plate_angle)
    
    # 应用相位延迟（慢轴分量相位滞后）
    E_slow_retarded = E_slow * np.exp(-1j * phase_retardation)
    
    # 转换回x-y坐标系
    output_Ex = E_fast * np.cos(plate_angle) - E_slow_retarded * np.sin(plate_angle)
    output_Ey = E_fast * np.sin(plate_angle) + E_slow_retarded * np.cos(plate_angle)
    return output_Ex.real, output_Ey.real

特别有趣的是λ/4波片的行为。当入射线偏振光与波片快轴成45°角时，会产生圆偏振光：

def demo_quarter_wave_plate():
    time = np.linspace(0, 4*np.pi, 200)
    Ex, Ey = linear_polarized_wave(1, np.pi/4, time)  # 45°偏振光
    
    # 通过λ/4波片
    Ex_qwp, Ey_qwp = wave_plate(Ex, Ey, 0, np.pi/2)  # 快轴沿x方向
    
    # 绘制结果
    plt.figure(figsize=(10, 5))
    plt.plot(time, Ex_qwp, 'r', label='Ex after QWP')
    plt.plot(time, Ey_qwp, 'b', label='Ey after QWP')
    plt.title('Circular Polarization Generated by λ/4 Wave Plate')
    plt.legend()
    plt.grid(True)

运行这段代码，你会看到x和y分量产生了π/2的相位差，这正是圆偏振光的特征。通过调整入射线偏振光与波片快轴的夹角，你还能观察到不同椭圆度的椭圆偏振光。

4. 完整实验系统模拟

现在我们把所有组件组合起来，构建一个完整的偏振光实验模拟系统：

class PolarizationExperiment:
    def __init__(self):
        self.polarizer_angle = 0      # 起偏器角度
        self.plate_angle = 0          # 波片角度
        self.plate_type = 'none'      # 波片类型：'none', 'qwp', 'hwp'
        self.analyzer_angle = 0       # 检偏器角度
    
    def set_components(self, polar_ang, plate_ang, plate_type, analyzer_ang):
        self.polarizer_angle = polar_ang
        self.plate_angle = plate_ang
        self.plate_type = plate_type
        self.analyzer_angle = analyzer_ang
    
    def run_experiment(self, time_points):
        # 生成自然光（随机偏振方向）
        angle = np.random.uniform(0, 2*np.pi)
        Ex, Ey = linear_polarized_wave(1, angle, time_points)
        
        # 通过起偏器
        Ex, Ey = polarizer(Ex, Ey, self.polarizer_angle)
        
        # 通过波片
        if self.plate_type == 'qwp':
            Ex, Ey = wave_plate(Ex, Ey, self.plate_angle, np.pi/2)
        elif self.plate_type == 'hwp':
            Ex, Ey = wave_plate(Ex, Ey, self.plate_angle, np.pi)
        
        # 通过检偏器
        Ex, Ey = polarizer(Ex, Ey, self.analyzer_angle)
        
        # 计算光强
        intensity = Ex**2 + Ey**2
        return intensity.mean()  # 平均光强

这个类可以模拟完整的偏振光实验流程。我们可以用它来验证马吕斯定律：

def verify_malus_law():
    experiment = PolarizationExperiment()
    angles = np.linspace(0, 180, 36)
    intensities = []
    
    for angle in angles:
        experiment.set_components(0, 0, 'none', angle)
        intensity = experiment.run_experiment(np.linspace(0, 2*np.pi, 100))
        intensities.append(intensity)
    
    # 绘制结果
    plt.figure(figsize=(10, 5))
    plt.plot(angles, intensities, 'bo-', label='Simulated Data')
    plt.plot(angles, np.cos(np.radians(angles))**2, 'r--', label='cos²θ')
    plt.xlabel('Analyzer Angle (degrees)')
    plt.ylabel('Relative Intensity')
    plt.title('Verification of Malus Law')
    plt.legend()
    plt.grid(True)

运行这段代码，你会看到模拟结果与马吕斯定律的理论预测完美吻合。这种通过代码验证物理定律的方式，比单纯看实验数据表格要直观得多。

5. 交互式可视化工具

为了让实验更加生动，我们可以使用ipywidgets创建一个交互式控制面板：

from ipywidgets import interact, FloatSlider, Dropdown

def interactive_polarization_simulator():
    def update(polarizer_angle=0, plate_angle=0, 
               plate_type='none', analyzer_angle=0):
        experiment = PolarizationExperiment()
        experiment.set_components(
            np.radians(polarizer_angle),
            np.radians(plate_angle),
            plate_type,
            np.radians(analyzer_angle)
        )
        
        # 模拟实时光波
        time = np.linspace(0, 2*np.pi, 100)
        Ex, Ey = linear_polarized_wave(1, np.radians(polarizer_angle), time)
        
        # 应用各组件
        if plate_type == 'qwp':
            Ex, Ey = wave_plate(Ex, Ey, np.radians(plate_angle), np.pi/2)
        elif plate_type == 'hwp':
            Ex, Ey = wave_plate(Ex, Ey, np.radians(plate_angle), np.pi)
        
        Ex, Ey = polarizer(Ex, Ey, np.radians(analyzer_angle))
        
        # 绘制结果
        plt.figure(figsize=(12, 4))
        plt.subplot(121)
        plt.plot(time, Ex, 'r', label='Ex')
        plt.plot(time, Ey, 'b', label='Ey')
        plt.ylim(-1.1, 1.1)
        plt.legend()
        plt.title('Electric Field Components')
        
        plt.subplot(122)
        plt.plot(Ex, Ey)
        plt.xlim(-1.1, 1.1)
        plt.ylim(-1.1, 1.1)
        plt.gca().set_aspect('equal')
        plt.title('Polarization State')
        plt.show()
    
    interact(update,
             polarizer_angle=FloatSlider(min=0, max=180, step=5, value=0),
             plate_angle=FloatSlider(min=0, max=180, step=5, value=0),
             plate_type=Dropdown(options=['none', 'qwp', 'hwp'], value='none'),
             analyzer_angle=FloatSlider(min=0, max=180, step=5, value=0))

这个交互工具让你可以实时调整各个光学元件的角度，立即看到偏振态的变化。比如你可以：

1. 设置起偏器为0°，检偏器为90°，观察消光现象
2. 在中间插入λ/4波片，旋转波片角度，观察如何产生椭圆偏振光
3. 使用λ/2波片，验证它会使偏振方向旋转两倍角度

6. 高级应用与扩展

掌握了基础模拟后，我们可以进一步扩展模型，模拟更复杂的偏振现象：

偏振光的干涉 ：当两束相干偏振光叠加时，干涉图样取决于它们的偏振状态。我们可以扩展模型来模拟这一现象：

def polarization_interference():
    # 创建两束相干线偏振光
    time = np.linspace(0, 4*np.pi, 500)
    Ex1, Ey1 = linear_polarized_wave(1, np.radians(0), time)   # 水平偏振
    Ex2, Ey2 = linear_polarized_wave(1, np.radians(60), time)  # 60°偏振
    
    # 添加相位差模拟光程差
    Ex2 = Ex2 * np.cos(0.5*time)
    Ey2 = Ey2 * np.cos(0.5*time)
    
    # 叠加两束光
    Ex_total = Ex1 + Ex2
    Ey_total = Ey1 + Ey2
    
    # 绘制结果
    plt.figure(figsize=(12, 5))
    plt.subplot(121)
    plt.plot(time, Ex_total, 'r', label='Ex')
    plt.plot(time, Ey_total, 'b', label='Ey')
    plt.title('Superposition of Two Polarized Waves')
    plt.legend()
    
    plt.subplot(122)
    plt.plot(Ex_total, Ey_total)
    plt.gca().set_aspect('equal')
    plt.title('Resultant Polarization State')
    plt.show()

琼斯矩阵形式化 ：对于更复杂的偏振光学系统，可以使用琼斯矩阵 formalism 来简化计算：

def jones_matrix_simulation():
    # 定义琼斯矩阵
    def polarizer_jones(angle):
        c, s = np.cos(angle), np.sin(angle)
        return np.array([[c*c, c*s], [c*s, s*s]])
    
    def waveplate_jones(angle, phase_retardation):
        c, s = np.cos(angle), np.sin(angle)
        return np.array([
            [c*c + s*s*np.exp(-1j*phase_retardation), c*s*(1 - np.exp(-1j*phase_retardation))],
            [c*s*(1 - np.exp(-1j*phase_retardation)), s*s + c*c*np.exp(-1j*phase_retardation)]
        ])
    
    # 初始光波（水平偏振）
    E_in = np.array([1, 0])
    
    # 构建光学系统：偏振片(0°) → λ/4波片(45°) → 偏振片(90°)
    P1 = polarizer_jones(0)
    QWP = waveplate_jones(np.pi/4, np.pi/2)
    P2 = polarizer_jones(np.pi/2)
    
    # 计算输出光波
    E_out = P2 @ QWP @ P1 @ E_in
    
    print(f"Output field: {E_out}")
    print(f"Output intensity: {np.abs(E_out[0])**2 + np.abs(E_out[1])**2}")

这种方法特别适合模拟由多个偏振元件组成的复杂光学系统。你可以轻松地添加或重新排列各种光学元件，而无需重写整个模拟逻辑。