一维光栅结构严格耦合波分析(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_diagclass 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 / wavelengthself.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.perioddef 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.periodreturn E, E_invdef 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# 计算矩阵Pif self.pol == 'TE':P = Velse: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, Lambdadef 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, Tdef 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_heighttotal_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, axdef rectangular_grating(period, width, n_ridge, n_groove):"""矩形光栅介电常数函数"""def epsilon(x):return np.where((x % period) < width, n_ridge**2, n_groove**2)return epsilondef 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))**2return epsilondef triangular_grating(period, n_min, n_max):"""三角光栅介电常数函数"""def epsilon(x):x_mod = x % periodreturn 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 epsilondef plot_results(solver, R, T):"""绘制RCWA计算结果"""plt.figure(figsize=(12, 8))# 绘制反射和透射光谱plt.subplot(2, 2, 1)orders = solver.nwidth = 0.4plt.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_Tplt.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.piplt.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 = 500width = 250n_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()# 求解RCWAR, T = solver.solve()# 绘制结果plot_results(solver, R, T)

代码功能说明

这个RCWA求解器实现了以下功能：

1. 核心算法：

   • 实现了严格的RCWA算法求解一维周期性光栅结构
   • 支持TE和TM偏振
   • 支持任意入射角度
   • 使用傅里叶展开处理周期性介电常数分布

2. 光栅类型支持：

   • 矩形光栅
   • 正弦光栅
   • 三角光栅
   • 可以轻松扩展支持其他光栅形状

3. 可视化功能：

   • 光栅结构可视化（显示各层厚度和介电常数分布）
   • 衍射效率柱状图（各阶反射和透射）
   • 能量守恒饼图（总反射、透射和吸收）
   • 角分布图（各衍射级的方向）
   • 参数显示（所有输入参数汇总）

   参考matlab代码 严格耦合波计算麦克斯韦方程数值解的源代码 youwenfan.com/contentcnb/77681.html，可以进行周期性的结构的数值求解，可以对一维所有类光栅结构进行求解。

应用

在示例中，我们创建了一个矩形光栅结构：

• 波长：600 nm
• 入射角：30度
• TE偏振
• 光栅周期：500 nm
• 光栅脊宽：250 nm
• 光栅高度：200 nm
• 光栅脊折射率：3.5（如硅）
• 光栅槽折射率：1.5（如玻璃）
• 衬底折射率：1.5（玻璃）

运行结果

运行代码将显示：

1. 光栅结构示意图，展示各层材料和介电常数分布
2. 衍射效率图，显示各阶反射和透射效率
3. 能量守恒饼图，显示总反射、透射和吸收
4. 角分布图，显示各衍射级的方向
5. 参数汇总表，包含所有输入参数