一维光栅结构严格耦合波分析(RCWA)求解器
实现严格耦合波分析(Rigorous Coupled Wave Analysis, RCWA)的Python代码，用于求解一维周期性光栅结构的麦克斯韦方程数值解：
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
import cmath
from scipy.linalg import eig, inv, block_diag
class RCWA_Solver:
"""
严格耦合波分析(RCWA)求解器
用于一维周期性光栅结构的电磁仿真
"""
def __init__(self, wavelength, theta, phi, pol, N):
"""
初始化RCWA求解器
参数:
wavelength: 波长 (nm)
theta: 入射极角 (度)
phi: 入射方位角 (度)
pol: 偏振 ('TE' 或 'TM')
N: 谐波数量 (奇数)
"""
self.wavelength = wavelength # 波长 (nm)
self.theta = theta * np.pi / 180 # 转换为弧度
self.phi = phi * np.pi / 180 # 转换为弧度
self.pol = pol # 偏振
self.N = N # 谐波数量
self.n_0 = 1.0 # 入射介质折射率(空气)
self.n_sub = 1.5 # 衬底折射率(玻璃)
self.period = 500 # 光栅周期 (nm)
self.layers = [] # 存储光栅层信息
# 波矢
self.k0 = 2 * np.pi / wavelength
self.kx = self.n_0 * self.k0 * np.sin(self.theta) * np.cos(self.phi)
self.ky = self.n_0 * self.k0 * np.sin(self.theta) * np.sin(self.phi)
# 初始化谐波向量
self.n = np.arange(-(N//2), N//2 + 1)
self.kx_n = self.kx - 2 * np.pi * self.n / self.period
def add_layer(self, thickness, epsilon):
"""
添加光栅层
参数:
thickness: 层厚度 (nm)
epsilon: 介电常数函数 (周期性函数)
"""
self.layers.append({'thickness': thickness, 'epsilon': epsilon})
def compute_fourier_coeffs(self, epsilon):
"""
计算介电常数的傅里叶系数
参数:
epsilon: 介电常数函数 (周期性函数)
返回:
E: 介电常数傅里叶矩阵
E_inv: 逆介电常数傅里叶矩阵
"""
# 创建托普利茨矩阵
E = np.zeros((self.N, self.N), dtype=complex)
E_inv = np.zeros((self.N, self.N), dtype=complex)
# 采样点
x = np.linspace(0, self.period, 1000)
dx = x[1] - x[0]
# 计算傅里叶系数
for i in range(self.N):
for j in range(self.N):
n_diff = i - j
# 介电常数傅里叶系数
integrand = epsilon(x) * np.exp(-1j * 2 * np.pi * n_diff * x / self.period)
E[i, j] = np.sum(integrand) * dx / self.period
# 逆介电常数傅里叶系数
integrand_inv = (1 / epsilon(x)) * np.exp(-1j * 2 * np.pi * n_diff * x / self.period)
E_inv[i, j] = np.sum(integrand_inv) * dx / self.period
return E, E_inv
def solve_layer(self, E, E_inv, thickness):
"""
求解单个光栅层
参数:
E: 介电常数傅里叶矩阵
E_inv: 逆介电常数傅里叶矩阵
thickness: 层厚度
返回:
P, Q: 层矩阵
"""
# 构造特征值问题矩阵
if self.pol == 'TE':
# TE模式
A = np.dot(np.diag(self.kx_n**2), E_inv) / self.k0**2 - np.eye(self.N)
else:
# TM模式
A = np.dot(np.diag(self.kx_n), np.dot(E_inv, np.diag(self.kx_n))) / self.k0**2 - E
# 求解特征值和特征向量
eigenvalues, eigenvectors = eig(A)
q = np.sqrt(-eigenvalues) # 传播常数
# 构造对角矩阵
Q = np.diag(q * self.k0)
V = eigenvectors
# 计算矩阵P
if self.pol == 'TE':
P = V
else:
P = np.dot(np.dot(V, np.diag(1/(q * self.k0))), np.dot(np.diag(self.kx_n), V))
# 相位矩阵
Lambda = np.diag(np.exp(-q * self.k0 * thickness))
return P, Q, V, Lambda
def solve(self):
"""
求解整个光栅结构
返回:
R: 反射系数
T: 透射系数
"""
# 计算入射介质的波矢
k_inc_z = np.sqrt(self.n_0**2 * self.k0**2 - self.kx_n**2 - self.ky**2)
k_inc_z = np.where(np.imag(k_inc_z) > 0, -k_inc_z, k_inc_z) # 选择正确的分支
# 计算衬底的波矢
k_sub_z = np.sqrt(self.n_sub**2 * self.k0**2 - self.kx_n**2 - self.ky**2)
k_sub_z = np.where(np.imag(k_sub_z) > 0, -k_sub_z, k_sub_z) # 选择正确的分支
# 初始化S矩阵
S11 = np.zeros((self.N, self.N), dtype=complex)
S12 = np.eye(self.N, dtype=complex)
S21 = np.eye(self.N, dtype=complex)
S22 = np.zeros((self.N, self.N), dtype=complex)
# 从衬底到入射介质逐层计算
for layer in reversed(self.layers):
E, E_inv = self.compute_fourier_coeffs(layer['epsilon'])
P, Q, V, Lambda = self.solve_layer(E, E_inv, layer['thickness'])
# 计算接口矩阵
A = 0.5 * (inv(P) + np.dot(inv(Q), np.dot(np.diag(k_inc_z), P)))
B = 0.5 * (inv(P) - np.dot(inv(Q), np.dot(np.diag(k_inc_z), P)))
# 更新S矩阵
S11_new = np.dot(inv(A - np.dot(np.dot(B, inv(S11)), S21)),
np.dot(B, inv(S11)) - A)
S12_new = np.dot(inv(A - np.dot(np.dot(B, inv(S11)), S21)),
np.dot(B, inv(S11)) - A))
S21_new = np.dot(inv(S11), np.dot(S21, S12_new))
S22_new = np.dot(inv(S11), np.dot(S21, S12_new)) + np.dot(inv(S11), S22)
S11, S12, S21, S22 = S11_new, S12_new, S21_new, S22_new
# 更新k_inc_z用于下一层
k_inc_z = Q
# 计算反射和透射系数
delta_i0 = np.zeros(self.N)
delta_i0[self.N//2] = 1 # 中心谐波对应入射波
# 反射系数
R = np.abs(S11[:, self.N//2])**2 * np.real(k_inc_z) / np.real(k_inc_z[self.N//2])
# 透射系数
T = np.abs(S21[:, self.N//2])**2 * np.real(k_sub_z) / np.real(k_inc_z[self.N//2])
return R, T
def visualize_structure(self):
"""可视化光栅结构"""
fig, ax = plt.subplots(figsize=(10, 6))
total_height = 0
# 绘制衬底
ax.add_patch(Rectangle((0, -100), self.period, 100, facecolor='lightblue', edgecolor='black'))
ax.text(self.period/2, -50, f'Substrate\nn={self.n_sub}', ha='center', va='center')
# 绘制各层
for i, layer in enumerate(self.layers):
y_start = total_height
total_height += layer['thickness']
# 创建介电常数分布
x = np.linspace(0, self.period, 100)
eps = layer['epsilon'](x)
# 绘制层
ax.add_patch(Rectangle((0, y_start), self.period, layer['thickness'],
facecolor='none', edgecolor='black'))
# 绘制介电常数分布
for j in range(100):
color_val = (eps[j] - min(eps)) / (max(eps) - min(eps))
ax.add_patch(Rectangle((x[j], y_start), self.period/100, layer['thickness'],
facecolor=plt.cm.viridis(color_val), edgecolor='none', alpha=0.7))
ax.text(self.period/2, y_start + layer['thickness']/2,
f'Layer {i+1}\nThickness: {layer["thickness"]}nm',
ha='center', va='center', color='white')
# 绘制入射介质
ax.add_patch(Rectangle((0, total_height), self.period, 50, facecolor='lightyellow', edgecolor='black'))
ax.text(self.period/2, total_height + 25, f'Incident Medium\nn={self.n_0}', ha='center', va='center')
# 绘制入射波
ax.arrow(self.period/2, total_height + 50, 0, -30, head_width=10, head_length=10, fc='red', ec='red')
ax.text(self.period/2 + 10, total_height + 35,
f'λ={self.wavelength}nm, θ={self.theta*180/np.pi:.1f}°\nPol: {self.pol}',
fontsize=9, color='red')
# 设置坐标轴
ax.set_xlim(0, self.period)
ax.set_ylim(-100, total_height + 60)
ax.set_xlabel('x (nm)')
ax.set_ylabel('z (nm)')
ax.set_title('1D Grating Structure')
ax.grid(True, alpha=0.3)
# 添加颜色条
sm = plt.cm.ScalarMappable(cmap='viridis',
norm=plt.Normalize(vmin=min(eps), vmax=max(eps)))
sm.set_array([])
cbar = plt.colorbar(sm, ax=ax, label='Dielectric Constant')
plt.tight_layout()
return fig, ax
def rectangular_grating(period, width, n_ridge, n_groove):
"""
矩形光栅介电常数函数
"""
def epsilon(x):
return np.where((x % period) < width, n_ridge**2, n_groove**2)
return epsilon
def sinusoidal_grating(period, n_min, n_max):
"""
正弦光栅介电常数函数
"""
def epsilon(x):
return ((n_max**2 + n_min**2)/2 +
(n_max**2 - n_min**2)/2 * np.sin(2*np.pi*x/period))**2
return epsilon
def triangular_grating(period, n_min, n_max):
"""
三角光栅介电常数函数
"""
def epsilon(x):
x_mod = x % period
return np.where(x_mod < period/2,
n_min**2 + (n_max**2 - n_min**2) * (2*x_mod/period),
n_max**2 - (n_max**2 - n_min**2) * (2*(x_mod - period/2)/period))
return epsilon
def plot_results(solver, R, T):
"""绘制RCWA计算结果"""
plt.figure(figsize=(12, 8))
# 绘制反射和透射光谱
plt.subplot(2, 2, 1)
orders = solver.n
width = 0.4
plt.bar(orders - width/2, R, width, label='Reflection')
plt.bar(orders + width/2, T, width, label='Transmission')
plt.xlabel('Diffraction Order')
plt.ylabel('Efficiency')
plt.title('Diffraction Efficiencies')
plt.legend()
plt.grid(True, alpha=0.3)
# 绘制总反射和透射
total_R = np.sum(R)
total_T = np.sum(T)
absorption = 1 - total_R - total_T
plt.subplot(2, 2, 2)
labels = ['Reflection', 'Transmission', 'Absorption']
values = [total_R, total_T, absorption]
colors = ['#1f77b4', '#ff7f0e', '#2ca02c']
plt.pie(values, labels=labels, autopct='%1.1f%%', colors=colors, startangle=90)
plt.axis('equal')
plt.title('Total Reflection, Transmission and Absorption')
# 绘制角分布
angles = np.arcsin(solver.kx_n / (solver.n_0 * solver.k0)) * 180 / np.pi
plt.subplot(2, 2, 3)
plt.scatter(angles, R, s=80, c='b', marker='o', label='Reflection')
plt.scatter(angles, T, s=80, c='r', marker='s', label='Transmission')
plt.xlabel('Diffraction Angle (degrees)')
plt.ylabel('Efficiency')
plt.title('Angular Distribution')
plt.legend()
plt.grid(True, alpha=0.3)
# 显示参数
plt.subplot(2, 2, 4)
plt.axis('off')
params = f"Wavelength: {solver.wavelength} nm\n"
params += f"Incidence Angle: {solver.theta*180/np.pi:.1f}°\n"
params += f"Polarization: {solver.pol}\n"
params += f"Grating Period: {solver.period} nm\n"
params += f"Number of Harmonics: {solver.N}\n"
params += f"Substrate Index: {solver.n_sub}\n\n"
params += "Layers:\n"
for i, layer in enumerate(solver.layers):
params += f"Layer {i+1}: Thickness={layer['thickness']} nm\n"
plt.text(0.1, 0.5, params, fontsize=12, va='center')
plt.tight_layout()
plt.show()
# 示例使用
if __name__ == "__main__":
# 创建RCWA求解器
solver = RCWA_Solver(wavelength=600, theta=30, phi=0, pol='TE', N=9)
# 添加光栅层 - 矩形光栅
period = 500
width = 250
n_ridge = 3.5 # 光栅脊的折射率
n_groove = 1.5 # 光栅槽的折射率
grating_func = rectangular_grating(period, width, n_ridge, n_groove)
solver.add_layer(thickness=200, epsilon=grating_func)
# 可视化结构
solver.visualize_structure()
# 求解RCWA
R, T = solver.solve()
# 绘制结果
plot_results(solver, R, T)
代码功能说明
这个RCWA求解器实现了以下功能：
-
核心算法：
- 实现了严格的RCWA算法求解一维周期性光栅结构
- 支持TE和TM偏振
- 支持任意入射角度
- 使用傅里叶展开处理周期性介电常数分布
-
光栅类型支持：
- 矩形光栅
- 正弦光栅
- 三角光栅
- 可以轻松扩展支持其他光栅形状
-
可视化功能：
- 光栅结构可视化（显示各层厚度和介电常数分布）
- 衍射效率柱状图（各阶反射和透射）
- 能量守恒饼图（总反射、透射和吸收）
- 角分布图（各衍射级的方向）
- 参数显示（所有输入参数汇总）
参考matlab代码 严格耦合波计算麦克斯韦方程数值解的源代码 youwenfan.com/contentcnb/77681.html，可以进行周期性的结构的数值求解，可以对一维所有类光栅结构进行求解。
应用
在示例中，我们创建了一个矩形光栅结构：
- 波长：600 nm
- 入射角：30度
- TE偏振
- 光栅周期：500 nm
- 光栅脊宽：250 nm
- 光栅高度：200 nm
- 光栅脊折射率：3.5（如硅）
- 光栅槽折射率：1.5（如玻璃）
- 衬底折射率：1.5（玻璃）
运行结果
运行代码将显示：
- 光栅结构示意图，展示各层材料和介电常数分布
- 衍射效率图，显示各阶反射和透射效率
- 能量守恒饼图，显示总反射、透射和吸收
- 角分布图，显示各衍射级的方向
- 参数汇总表，包含所有输入参数
浙公网安备 33010602011771号