梯度折射率透镜的光程计算方法

周末在arxiv上闲逛时，看到一篇关于GRIN（Gradient Index）透镜的论文。 梯度折射率这个概念挺有意思的，虽然之前没接触过，但看起来用编程来理解应该不难。 这次尝试整理了一下学到的知识，顺便用Python算了几个例子。

1. 什么是GRIN透镜？

普通透镜是靠表面的曲率来聚焦光线，折射率在整个透镜内部是均匀的。 而GRIN透镜恰恰相反：表面可以是平的，但内部折射率是连续变化的。

这种设计有几个好处：

• 可以减少表面反射损耗
• 能够更灵活地控制光路
• 在某些情况下可以做得更紧凑

典型的应用比如光纤耦合、内窥镜、复印机的成像系统等。

2. 折射率分布模型

最常见的是抛物线型折射率分布（也叫selfoc透镜）。 在径向上，折射率按照下面的规律变化：

其中：
• n₀ 是轴上（r=0）的折射率
• r 是到光轴的距离
• A 是梯度常数，决定了折射率变化的速度

这种分布会让光线在透镜内部沿正弦曲线传播，挺有意思的。

3. 费马原理与光程

在均匀介质中，光走直线，光程就是 n×d（折射率×几何距离）。 但在变折射率介质中，光的路径本身就是曲线，怎么算光程呢？

对于任意路径，光程可以表示为积分：

这里的积分是沿着光线路径的线积分。在GRIN介质中，需要先求出光线路径，然后再积分计算光程。

4. 用Python求解光线路径

在抛物线型GRIN中，光线方程可以解析求解，但更通用的方法是用数值方法。 我用scipy的odeint来求解光线传播的微分方程。

光线方程

根据几何光学，光线的传播可以用二阶微分方程描述：

把它改写成一阶方程组，就可以用常微分方程求解器来算了：

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# 折射率分布函数
def n(r, n0=1.5, A=0.1):
    return n0 * (1 - A * r**2 / 2)

# 折射率梯度
def dn_dr(r, n0=1.5, A=0.1):
    return -n0 * A * r

# 光线方程（一阶形式）
def ray_eq(y, z, n0=1.5, A=0.1):
    r, dr_dz = y
    d2r_dz2 = (1 / n(r, n0, A)) * dn_dr(r, n0, A)
    return [dr_dz, d2r_dz2]

# 初始条件：r0=0.5mm, 初始角度2度
r0 = 0.5e-3  # 起始位置
angle = np.radians(2)  # 入射角
dr_dz0 = np.tan(angle)  # 初始斜率

# 求解
z = np.linspace(0, 10e-3, 1000)  # 传播10mm
solution = odeint(ray_eq, [r0, dr_dz0], z)

# 绘制光线轨迹
plt.figure(figsize=(10, 6))
plt.plot(z * 1e3, solution[:, 0] * 1e3, linewidth=2)
plt.xlabel('传播距离 z (mm)')
plt.ylabel('径向位置 r (mm)')
plt.title('GRIN透镜中的光线轨迹')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

运行后可以看到光线在GRIN透镜中的正弦曲线轨迹，周期性地靠近和远离光轴。

5. 计算光程

有了光线轨迹r(z)之后，就可以计算光程了。微元段的光程是：

对整个路径积分：

# 计算光程
def calculate_opl(z, r, dr_dz, n0=1.5, A=0.1):
    # 计算每段的光程
    ds = np.sqrt(1 + dr_dz**2)  # 路径微元
    n_values = n(r, n0, A)      # 折射率
    opl = np.trapz(n_values * ds, z)  # 数值积分
    return opl

# 从求解结果中提取
r_path = solution[:, 0]
dr_dz_path = solution[:, 1]

opl = calculate_opl(z, r_path, dr_dz_path)
print(f"总光程: {opl*1e3:.6f} mm")
print(f"几何距离: {z[-1]*1e3:.6f} mm")
print(f"等效折射率: {opl/z[-1]:.6f}")

运行结果示例：

可以看到，光程大于几何距离，这是因为光在传播过程中经过了不同折射率的区域。

6. 可视化折射率分布

为了更直观，我还画了折射率的2D分布图：

# 创建2D网格
r_grid = np.linspace(-2e-3, 2e-3, 200)
z_grid = np.linspace(0, 10e-3, 400)
R, Z = np.meshgrid(r_grid, z_grid)

# 计算折射率分布
N = n(np.abs(R), n0=1.5, A=0.1)

# 绘制
plt.figure(figsize=(12, 6))
plt.contourf(Z*1e3, R*1e3, N, levels=50, cmap='viridis')
plt.colorbar(label='折射率 n')

# 叠加光线轨迹
plt.plot(z*1e3, r_path*1e3, 'r-', linewidth=2, label='光线轨迹')
plt.plot(z*1e3, -r_path*1e3, 'r-', linewidth=2)

plt.xlabel('传播距离 z (mm)')
plt.ylabel('径向位置 r (mm)')
plt.title('GRIN透镜中的折射率分布与光线轨迹')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

这张图很清楚地展示了光线是如何在折射率梯度的作用下弯曲的。 轴心区域（r=0）折射率最高，边缘区域逐渐降低。

7. 实际应用中的考虑

虽然理论上很漂亮，但实际制造GRIN透镜还是有难度的：

制造工艺

常用的方法包括离子交换、化学气相沉积等。要精确控制折射率分布并不容易， 而且只能做成特定的形状（比如圆柱形）。

像差问题

虽然GRIN可以减少某些像差，但也会引入新的像差。 实际设计时需要仔细优化折射率分布曲线。

成本

GRIN透镜通常比普通透镜贵不少，只在性能要求高或空间受限的场合才用。

8. 总结与思考

通过这次学习，有几点收获：

• 变折射率介质中的光程计算需要数值积分
• Python的科学计算库（numpy、scipy）处理这类问题很方便
• 可视化能帮助理解抽象的概念
• 理论和实际应用之间还是有不小的gap

下一步想看看能不能把这个方法推广到更复杂的折射率分布， 比如非轴对称的情况。可能需要用到有限元或光线追迹软件了。