✅作者简介:热爱科研的Matlab仿真开发者,修心和技术同步精进,代码获取、论文复现及科研仿真合作可私信。
🍎个人主页:Matlab科研工作室
🍊个人信条:格物致知。
更多Matlab完整代码及仿真定制内容点击👇
🔥 内容介绍
超材料是一种具有独特电磁特性的新型人工材料,其结构尺寸远小于入射光的波长。三维超材料由于其结构的复杂性,可以实现对光波的更精细控制,从而在光学领域具有广泛的应用前景。本文研究了三维超材料对光的衍射特性,分析了衍射引起的反射、透射和场分布的变化规律。
引言
光波在通过周期性结构时会发生衍射,产生多个衍射波。三维超材料的周期性结构可以有效地控制衍射波的传播方向和强度,从而实现对光波的调控。本文通过数值模拟的方法,研究了三维超材料对光的衍射特性,分析了衍射引起的反射、透射和场分布的变化规律。
数值模拟方法
本文采用有限元方法对三维超材料的衍射特性进行数值模拟。模拟模型为一个周期性的三维超材料结构,入射光为平面波。通过求解麦克斯韦方程组,可以得到超材料结构对入射光的衍射场分布。
结果与讨论
反射率和透射率
图1给出了三维超材料结构的反射率和透射率随入射角的变化曲线。可以看出,当入射角较小时,超材料结构的反射率较低,透射率较高。随着入射角的增加,反射率逐渐增加,透射率逐渐减小。在某些特定的入射角下,超材料结构可以实现对入射光的全反射或全透射。
场分布
图2给出了三维超材料结构在不同入射角下的场分布。可以看出,当入射角较小时,衍射波主要集中在超材料结构的表面附近。随着入射角的增加,衍射波逐渐向超材料结构内部传播。在某些特定的入射角下,衍射波可以在超材料结构内部形成驻波。
结论
本文通过数值模拟的方法,研究了三维超材料对光的衍射特性,分析了衍射引起的反射、透射和场分布的变化规律。研究结果表明,三维超材料可以有效地控制衍射波的传播方向和强度,从而实现对光波的调控。这些特性为三维超材料在光学器件、光通信和光计算等领域提供了广阔的应用前景。波导:**三维超材料可以作为光波导,引导和控制光在特定路径中传播,实现光信号的传输和处理。
-
**隐形斗篷:**通过调控超材料的衍射特性,可以实现对物体的光学隐形,使物体在特定波段下对光波透明。
结论
三维超材料光的衍射引起的反射、透射和场分布特性是影响其光学性能的关键因素。通过理解和控制这些衍射效应,可以设计和优化三维超材料的光学器件,实现各种光学功能,推动光学技术的发展。随着三维超材料研究的不断深入,其在光学领域的应用潜力将进一步得到挖掘和拓展。
📣 部分代码
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Inputs and Structure Parameters
y = 500; % Range of wavelength in nm
% Angles
th = 0; % Incidence Angle (deg)
phi = 0; % Angle between 'x' axis and plane of incidence (deg)
psi = 90; % Polarization angle (deg) 0 -> TM, 90 -> TE
% Unit cell parameters
mlt = 1; % Scaling factor for the hexagonal unit cell
D_x = 300; % Length of one period in 'x'
D_x = 2*round(D_x*mlt/2); % Scaled and rounded to nearest even number
D_y = 2*round(sqrt(3)*D_x/2); % Length of one period in 'y'
% Metal-dielectric metamaterials
D_m = 10; % Thickness of metal layer in 'z'
D_d = 30; % Thickness of dielectric layer in 'z'
NN_z = 6; % =10; for uncarved % Number of uncarved unit cells (after the moth-eyes)
% Moth-eye Structures
N_m = 4; % =0; for uncarved % Number of elementary cells carved into pillars
r_1 = 50*mlt; % Radius of pillars at top
r_2 = 145*mlt; % Radius of pillars at bottom
N_zd = 3; % Number of slices of the dielectric layer (calculates radius, increase for accuracy)
D_z = N_m*(2*D_d+D_m); % Thickness of moth-eye structure in 'z'
N_z = N_m*(2*N_zd+1); % Number of layers in one carved elementary cell
h_zd = D_d/N_zd; % Resolution in 'z' for carved dielectric layers
h_z = repmat([repmat(h_zd,1,N_zd),D_m,repmat(h_zd,1,N_zd)],1,N_m); % Height of each carved layer
z_l = [0,cumsum(h_z(1:N_z-1))]; % 'z' coordinate of each carved layer
r_z = interp1([0,D_z-h_zd],[r_1,r_2],z_l); % Radius of pillars for each carved layer
N_l = N_z+3*NN_z+2; % Total layers
% Fourier modes to consider in calculations (ODD)
N_x = 9; %=1; for uncarved % Number of Fourier modes in x direction
N_y = 9; %=1; for uncarved % Number of Fourier modes in y direction
Nt = N_x*N_y;
% Preallocations
ky = 2*pi./y; % Wavenumbers
I = eye(Nt); II = eye(2*Nt); % Identity matrices
ZZ = zeros(2*Nt); zz = zeros(2*Nt,1);
EPS = zeros(Nt,Nt,N_l);
EPSyx = zeros(Nt,Nt,N_z);
EPSxy = zeros(Nt,Nt,N_z);
W = zeros(2*Nt,2*Nt,N_l);
V = zeros(2*Nt,2*Nt,N_l);
X = zeros(2*Nt,2*Nt,N_l);
qm = zeros(2*Nt,N_l);
% Foruier mode numbers and arrays
mx = (-(N_x-1)/2:(N_x-1)/2)';
my = (-(N_y-1)/2:(N_y-1)/2);
m0x = (N_x+1)/2;
m0y = (N_y+1)/2;
m0 = (Nt+1)/2;
%% Wavelength based permittivity
% Metal (with thin layer correction)
load Cu.txt; % Taken from refractiveindex.info
Me = interp1(Cu(:,1),Cu,y/1000); % Interpolation of experimental data
Metal = (Me(:,2) + 1i*Me(:,3)); Metal = Metal.';
kkk = D_m/39; % Thickness of W layer divided mean free path of W
kky = 1000*ky;
fi = (1./kkk+3/4*(1-1/12*kkk.^2).*expint(kkk)-3./(8*kkk.^2).*(1-exp(-kkk))-(5./(8*kkk)+1/16-kkk/16).*exp(-kkk)).^-1;
wp = 6.5; % Plasma frequency in micron^-1
gamma = 12*10^-2; gammakk = gamma*fi/kkk;
epsDrude = 1-(wp^2)./((kky/(2*pi)).^2+1i*(kky/(2*pi))*gamma); % Drude model for the Bulk
epsDrudekk = 1-(wp^2)./((kky/(2*pi)).^2+1i*(kky/(2*pi))*gammakk); % Drude Model corrected for the thin layer
epsMetal = Metal-epsDrude+epsDrudekk;
%% Index profile of moth eye structures
% Build the moth-eye structure as a binary 3D matrix
indexProfLog = false(D_x,D_y,N_z);
% Hexagon vertice positions
centre = [D_x,D_y]/2;
for nn=1:N_z
r = round(r_z(nn));
% Pillars
[rr,cc] = meshgrid(1:2*r+1);
circ = (rr-1-r).^2+(cc-1-r).^2<=r^2;
% Positioning the pillars
indexProfLog((centre(1)-r):(centre(1)+r),(centre(2)-r):(centre(2)+r),nn) = circ;
indexProfLog(1:r,1:r,nn) = indexProfLog(1:r,1:r,nn) | circ(r+2:2*r+1,r+2:2*r+1);
indexProfLog(1:r,D_y-r:D_y,nn) = indexProfLog(1:r,D_y-r:D_y,nn) | circ(r+2:2*r+1,1:r+1);
indexProfLog(D_x-r:D_x,D_y-r:D_y,nn) = indexProfLog(D_x-r:D_x,D_y-r:D_y,nn) | circ(1:r+1,1:r+1);
indexProfLog(D_x-r:D_x,1:r,nn) = indexProfLog(D_x-r:D_x,1:r,nn) | circ(1:r+1,r+2:2*r+1);
end
%% Structure Analysis
uy = y/1000; % Wavelength in micron
% Permittivity and refractive indices
n_1 = 1; % Refractive index of incidence medium
n_3 = sqrt(1 + 0.28604141 + 1.07044083/(1 - 1.00585997e-2/uy^2) +...
1.10202242/(1 - 100/uy^2)); % Refractive index for transmission medium (quartz)
e_m = epsMetal; % Permittivty of metal for current wavelength
e_d = 1 + 0.6961663/(1-(0.0684043/uy)^2) + ...
0.4079426/(1 - (0.1162414/uy)^2) + 0.8974794/(1 - (9.896161/uy)^2); % Permittivity of dielectric layer (SiO2)
e_f = 1; % Permittivity of material between pillars
e_p = repmat([repmat(conj(e_d),N_zd,1);conj(e_m);repmat(conj(e_d),N_zd,1)],N_m,1);
% Permittivity matrices for each layer
indexProf = ones(D_x,D_y);
for nn=1:N_z
indexProf(indexProfLog(:,:,nn)) = e_p(nn);
indexProf(~indexProfLog(:,:,nn)) = e_f;
eps_fft = fftshift(fft2(indexProf))/(D_x*D_y);
eps_mn = eps_fft(D_x/2+2-N_x:D_x/2+N_x,D_y/2+2-N_y:D_y/2+N_y);
for pp = 1:N_x
for qq = 1:N_y
EE = rot90(eps_mn(pp:pp+N_x-1,qq:qq+N_y-1),2);
EPS(pp+N_x*(qq-1),:,nn+1) = reshape(EE,1,[]);
end
end
i_iepsx_mj = zeros(N_x,D_y,N_x);
iepsx_fft = fftshift(fft(indexProf.^(-1),[],1),1)/D_x;
for qq=1:D_y
iepsx_m = iepsx_fft(D_x/2+2-N_x:D_x/2+N_x,qq);
iepsx_mj = toeplitz(iepsx_m(N_x:2*N_x-1),flip(iepsx_m(1:N_x)));
i_iepsx_mj(:,qq,:) = inv(iepsx_mj);
end
epsxy_fft = fftshift(fft(i_iepsx_mj,[],2),2)/D_y;
epsxy_mnj = epsxy_fft(:,D_y/2+2-N_y:D_y/2+N_y,:);
E4 = zeros(N_x,N_y,N_x,N_y);
for pp = 1:N_x
for qq = 1:N_x
E4(pp,:,qq,:) = toeplitz(epsxy_mnj(pp,N_y:2*N_y-1,qq),flip(epsxy_mnj(pp,1:N_y,qq)));
end
end
EPSxy(:,:,nn) = reshape(E4,[N_x*N_y,N_x*N_y]);
i_iepsy_nl = zeros(D_x,N_y,N_y);
iepsy_fft = fftshift(fft(indexProf.^(-1),[],2),2)/D_y;
for pp=1:D_x
iepsy_n = iepsy_fft(pp,D_y/2+2-N_y:D_y/2+N_y);
iepsy_nl = toeplitz(iepsy_n(N_y:2*N_y-1),flip(iepsy_n(1:N_y)));
i_iepsy_nl(pp,:,:) = inv(iepsy_nl);
end
epsyx_fft = fftshift(fft(i_iepsy_nl,[],1),1)/D_x;
epsyx_mnl = epsyx_fft(D_x/2+2-N_x:D_x/2+N_x,:,:);
E4 = zeros(N_x,N_y,N_x,N_y);
for rr = 1:N_y
for ss = 1:N_y
E4(:,rr,:,ss) = toeplitz(epsyx_mnl(N_x:2*N_x-1,rr,ss),flip(epsyx_mnl(1:N_x,rr,ss)));
end
end
EPSyx(:,:,nn) = reshape(E4,[N_x*N_y,N_x*N_y]);
end
%% Initializing variables
% Incident Field
u_x = cosd(psi)*cosd(phi)*cosd(th) - sind(psi)*sind(phi);
u_y = cosd(psi)*sind(phi)*cosd(th) + sind(psi)*cosd(phi);
inc = zz; inc(m0) = u_x; inc(Nt+m0) = u_y;
% Wavenumber Matrices
k_0 = ky;
k_xi = k_0*n_1*sind(th)*cosd(phi) + 2*pi*mx/D_x;
k_yi = k_0*n_1*sind(th)*sind(phi) + 2*pi*my/D_y;
k_x_mn = reshape(repmat(k_xi,1,N_y),[],1);
k_y_mn = reshape(repmat(k_yi,N_x,1),[],1);
k_1_zmn = conj(((n_1*k_0)^2*ones(Nt,1) - k_x_mn.^2 - k_y_mn.^2).^(1/2));
k_d_zmn = conj((e_d*(k_0)^2*ones(Nt,1) - k_x_mn.^2 - k_y_mn.^2).^(1/2));
k_m_zmn = conj((e_m*(k_0)^2*ones(Nt,1) - k_x_mn.^2 - k_y_mn.^2).^(1/2));
k_3_zmn = conj(((n_3*k_0)^2*ones(Nt,1) - k_x_mn.^2 - k_y_mn.^2).^(1/2));
K_x = 1i*diag(k_x_mn/k_0);
K_y = 1i*diag(k_y_mn/k_0);
%% Eigenmatrices
for nn=1:N_z
% Eigenvalue formulation for the grating region
F11 = -K_x*(EPS(:,:,nn+1)\K_y);
F12 = I + K_x*(EPS(:,:,nn+1)\K_x);
F21 = -I - K_y*(EPS(:,:,nn+1)\K_y);
F22 = K_y*(EPS(:,:,nn+1)\K_x);
F = [F11 F12;F21 F22];
G11 = -K_x*K_y;
G12 = EPSyx(:,:,nn) + K_x^2;
G21 = -EPSxy(:,:,nn) - K_y^2;
G22 = K_x*K_y;
G = [G11 G12;G21 G22];
[W(:,:,nn+1), Q] = eig(F*G);
Q = (sqrt(Q));
qm(:,nn+1) = diag(Q);
V(:,:,nn+1) = -G*W(:,:,nn+1)/Q;
X(:,:,nn+1) = diag(exp(-k_0*qm(:,nn+1)*h_z(nn)));
end
% Eigenmatrices for reflection and transmission layers
W(:,:,1) = II;
W(:,:,N_l) = II;
EPS(:,:,1) = n_1^2*I;
VIxx = (1i/k_0)*diag((k_x_mn.*k_y_mn)./k_1_zmn);
VIxy = (1i/k_0)*diag((k_y_mn.^2 + k_1_zmn.^2)./k_1_zmn);
VIyx = -(1i/k_0)*diag((k_x_mn.^2 + k_1_zmn.^2)./k_1_zmn);
VIyy = -VIxx;
V(:,:,1) = [VIxx VIxy;VIyx VIyy];
qm(:,1) = 1i*[k_1_zmn;k_1_zmn]/k_0;
X(:,:,1) = II;
EPS(:,:,N_l) = n_3^2*I;
VIIxx = (1i/k_0)*diag((k_x_mn.*k_y_mn)./k_3_zmn);
VIIxy = (1i/k_0)*diag((k_y_mn.^2 + k_3_zmn.^2)./k_3_zmn);
VIIyx = -(1i/k_0)*diag((k_x_mn.^2 + k_3_zmn.^2)./k_3_zmn);
VIIyy = -VIIxx;
V(:,:,N_l) = [VIIxx VIIxy;VIIyx VIIyy];
qm(:,N_l) = 1i*[k_3_zmn;k_3_zmn]/k_0;
X(:,:,N_l) = II;
% Eigenmatrices for the uncarved layers
for nn=1:3*NN_z
if (mod(nn+1,3))==0
k_n_zmn = k_m_zmn;
D_n = D_m;
EPS(:,:,1+N_z+nn) = conj(e_m)*I;
else
k_n_zmn = k_d_zmn;
D_n = D_d;
EPS(:,:,1+N_z+nn) = conj(e_d)*I;
end
W(:,:,1+N_z+nn) = II;
Vnxx = (1i/k_0)*diag((k_x_mn.*k_y_mn)./k_n_zmn);
Vnxy = (1i/k_0)*diag((k_y_mn.^2 + k_n_zmn.^2)./k_n_zmn);
Vnyx = -(1i/k_0)*diag((k_x_mn.^2 + k_n_zmn.^2)./k_n_zmn);
Vnyy = -Vnxx;
V(:,:,1+N_z+nn) = [Vnxx Vnxy;Vnyx Vnyy];
qm(:,1+N_z+nn) = 1i*[k_n_zmn;k_n_zmn]/k_0;
X(:,:,1+N_z+nn) = diag(exp(-1i*[k_n_zmn;k_n_zmn]*D_n));
end
%% Big matrix calculation
% Preallocation
B = repmat(ZZ,2*(N_l+1),2*(N_l+1));
% Building B using b
for ll=1:N_l
b = [W(:,:,ll) W(:,:,ll)*X(:,:,ll);
V(:,:,ll) -V(:,:,ll)*X(:,:,ll);
⛳️ 运行结果
🔗 参考文献
-
M. G. Moharam and T. K. Gaylord, “Rigorous coupled-wave analysis of planar-grating diffraction,” J. Opt. Soc. Am. 71, 811-818 (1995).
-
L. Li, “New formulation of the Fourier modal method for crossed surface-relief gratings,” J. Opt. Soc. Am. A 14, 2758-2767 (1997).
-
L. Li, “Formulation and comparison of two recursive matrix algorithms for modeling layered diffraction gratings,” J. Opt. Soc. Am. A 13(5), 1024-1035 (1996).
🎈 部分理论引用网络文献,若有侵权联系博主删除
🎁 关注我领取海量matlab电子书和数学建模资料
👇 私信完整代码和数据获取及论文数模仿真定制
1 各类智能优化算法改进及应用
生产调度、经济调度、装配线调度、充电优化、车间调度、发车优化、水库调度、三维装箱、物流选址、货位优化、公交排班优化、充电桩布局优化、车间布局优化、集装箱船配载优化、水泵组合优化、解医疗资源分配优化、设施布局优化、可视域基站和无人机选址优化、背包问题、 风电场布局、时隙分配优化、 最佳分布式发电单元分配、多阶段管道维修、 工厂-中心-需求点三级选址问题、 应急生活物质配送中心选址、 基站选址、 道路灯柱布置、 枢纽节点部署、 输电线路台风监测装置、 集装箱船配载优化、 机组优化、 投资优化组合、云服务器组合优化、 天线线性阵列分布优化、CVRP问题、VRPPD问题、多中心VRP问题、多层网络的VRP问题、多中心多车型的VRP问题、 动态VRP问题、双层车辆路径规划(2E-VRP)、充电车辆路径规划(EVRP)、油电混合车辆路径规划、混合流水车间问题、 订单拆分调度问题、 公交车的调度排班优化问题、航班摆渡车辆调度问题、选址路径规划问题
2 机器学习和深度学习方面
2.1 bp时序、回归预测和分类
2.2 ENS声神经网络时序、回归预测和分类
2.3 SVM/CNN-SVM/LSSVM/RVM支持向量机系列时序、回归预测和分类
2.4 CNN/TCN卷积神经网络系列时序、回归预测和分类
2.5 ELM/KELM/RELM/DELM极限学习机系列时序、回归预测和分类
2.6 GRU/Bi-GRU/CNN-GRU/CNN-BiGRU门控神经网络时序、回归预测和分类
2.7 ELMAN递归神经网络时序、回归\预测和分类
2.8 LSTM/BiLSTM/CNN-LSTM/CNN-BiLSTM/长短记忆神经网络系列时序、回归预测和分类
2.9 RBF径向基神经网络时序、回归预测和分类