RCWA计算一维光栅衍射效率

clc;clear all; %% parameter c0 = 3e8; % light speed theta = 0 * pi /180; % incident angle lamda = 0.5; % wavelength w = 2 * pi * c0./ lamda; % angle frequency order = 1; order_max = order; % diffraction order max order_min = -order; % diffraction order min d = 0.5; % thickness of each layer from front to back, [um] N = length(d); % number of layers perior = 1; % perior k0 = 2 * pi ./ lamda; nrd = 1.5; % Refractive index of the convex parts of the grating ngr = 1; % Refractive index of the concave parts of the grating i = order_min:order_max; %%% e的表达式 e: permittivity of incident e(1) and substrate media e(2) e(1) = 1; % incident medium permitivity e(2) = 1.5.^2; % out medium permitivity kxi = k0 .* (sqrt(e(1)).* sin(theta) - i * (lamda /perior)); % f1 = [0,0]; f = 0.5; %% grating field permitivity fourier expansion % e0 = nrd^2 * f + ngr^2 * (1 - f); % h = 1:order_max - order_min; % eh = (em - ed) .* sin(pi * h * f)./(pi * h); % p = order_min - order_max:-1; % ep = (em - ed) .* sin(pi * h * f)./(pi * h); for h = 1:length(i) if i(h) == 0 e_grating(h) = nrd^2 * f + ngr^2 * (1 - f); else e_grating(h) = (nrd ^ 2 - ngr ^ 2) * sin(f * pi * i(h))/(i(h) * pi); end end % e_g = [ep e0 eh]; E = zeros(length(i),length(i)); range = max(i) - min(i); mod = mod(range,2); if mod == 0 mod = mod + 2; else mod = mod + 1; end for s = 1:length(i) for p = 1:length(i) %下标计算 if s - p > order_max || s - p < order_min E(s,p) = 0; else E(s,p) = e_grating(s-p+((range + mod)/2)); end end end %% wavevector in z direction if k0 .* sqrt(e(1)) > kxi kz1 = sqrt(e(1) .* k0.^2 - kxi.^2); else kz1 = -1j .*sqrt((kxi).^2 - e(1) .* k0.^2); end if k0 .* sqrt(e(1)) > kxi kz2 = sqrt(e(2) .* k0.^2 - kxi.^2); else kz2 = -1j .*sqrt((kxi).^2 - e(2) .* k0.^2); end kzz = [kz1;kz2]; %% Eigenvalues and eigenvectors of coupled wave equations I = eye(order * 2 +1); Kx = diag(kxi ./ k0); A = Kx.^2 - E; [W , Q] = eig(A); %Q:eigenvalue, W:corresponding vector Qm = sqrt(Q); V = W * Qm; %% The amplitude coefficient of propagation under the interface is solved delta = zeros(2 * order + 1, 1); delta(1 + order) = 1; X = zeros(2 * order +1, 2 * order +1); for x = 1:length(X) X(x,x) = exp(-k0 .* Qm(x,x) .* d); end % Yinc = diag(kz1./k0); %First medium: air (or dielectric) % Ysub = diag(kz2./k0); %Last medium: substrate % Yinc = kz1./k0; %First medium: air (or dielectric) % Ysub = kz2./k0; %Last medium: substrate Y1 = zeros(2 * order +1, 2 * order +1); for x = 1:length(Y1) Y1(x,x) = kz1(x)./k0; % First medium: air (or dielectric) end Y2 = zeros(2 * order +1, 2 * order +1); for x = 1:length(Y2) Y2(x,x) = kz2(x)./k0; % First medium: air (or dielectric) end M_1_1 = 1j * Y1 * W + V; M_1_2 = (1j * Y1 * W - V) * X; M_2_1 = (V - 1j * Y2 * W) * X; M_2_2 = -V - 1j * Y2 * W; matrix_pro = [M_1_1,M_1_2;M_2_1,M_2_2]; solver_up = [1j * Y1 * delta + 1j * sqrt(e(1)) * cos(theta) * delta]; solver_down = zeros(2*order+1,1); solver = [solver_up;solver_down]; %[3*3 3*3] [3*3] % *[C1;C2] = %[3*3 3*3] [3*3] 求解得到的C为6*1，分别对应-1：1：1的C+和C- C = inv(matrix_pro) * solver; %% Normalized amplitudes of class i reflected and transmitted waves Ti = ( W * X * C(1:2*order+1,1) + W * C(2*order+2:2*(2*order+1),1)) ; Ri = ( W * C(1:2*order+1,1) + W * X * C(2*order+2:2*(2*order+1),1)) - delta; %% Diffraction efficiency of grating series = real((kz1')./(k0.*sqrt(e(1)).*cos(theta))); DEr = Ri.*conj(Ri).* series; DEt = Ti.*conj(Ti).*real((kz2')./(k0.*sqrt(e(1)).*cos(theta))); RR = Ri.*conj(Ri); R = sum(DEr); T = sum(DEt); diffraction_all = sum(DEr + DEt); %% Normalization is used for lumerical comparison diffraction_T = abs(DEt) ./sum(abs(DEt)); diffraction_R = abs(DEr) ./sum(abs(DEr)); RR_obj = real((kz1')./(k0.*sqrt(e(1)).*cos(theta))); % Diff_lumerical_grating_T=[ % 0.284915 % 0.100292 % 0.614792]; %lumerical_grating_R=[0.459171 0.0816573 0.459171]