✅作者简介：热爱科研的Matlab仿真开发者，修心和技术同步精进，代码获取、论文复现及科研仿真合作可私信。
🍎个人主页：Matlab科研工作室
🍊个人信条：格物致知。
更多Matlab完整代码及仿真定制内容点击👇
🔥 内容介绍
超材料是一种具有独特电磁特性的新型人工材料，其结构尺寸远小于入射光的波长。三维超材料由于其结构的复杂性，可以实现对光波的更精细控制，从而在光学领域具有广泛的应用前景。本文研究了三维超材料对光的衍射特性，分析了衍射引起的反射、透射和场分布的变化规律。
引言
光波在通过周期性结构时会发生衍射，产生多个衍射波。三维超材料的周期性结构可以有效地控制衍射波的传播方向和强度，从而实现对光波的调控。本文通过数值模拟的方法，研究了三维超材料对光的衍射特性，分析了衍射引起的反射、透射和场分布的变化规律。
数值模拟方法
本文采用有限元方法对三维超材料的衍射特性进行数值模拟。模拟模型为一个周期性的三维超材料结构，入射光为平面波。通过求解麦克斯韦方程组，可以得到超材料结构对入射光的衍射场分布。
结果与讨论
反射率和透射率
图1给出了三维超材料结构的反射率和透射率随入射角的变化曲线。可以看出，当入射角较小时，超材料结构的反射率较低，透射率较高。随着入射角的增加，反射率逐渐增加，透射率逐渐减小。在某些特定的入射角下，超材料结构可以实现对入射光的全反射或全透射。
场分布
图2给出了三维超材料结构在不同入射角下的场分布。可以看出，当入射角较小时，衍射波主要集中在超材料结构的表面附近。随着入射角的增加，衍射波逐渐向超材料结构内部传播。在某些特定的入射角下，衍射波可以在超材料结构内部形成驻波。
结论
本文通过数值模拟的方法，研究了三维超材料对光的衍射特性，分析了衍射引起的反射、透射和场分布的变化规律。研究结果表明，三维超材料可以有效地控制衍射波的传播方向和强度，从而实现对光波的调控。这些特性为三维超材料在光学器件、光通信和光计算等领域提供了广阔的应用前景。波导：**三维超材料可以作为光波导，引导和控制光在特定路径中传播，实现光信号的传输和处理。
-
**隐形斗篷：**通过调控超材料的衍射特性，可以实现对物体的光学隐形，使物体在特定波段下对光波透明。
结论
三维超材料光的衍射引起的反射、透射和场分布特性是影响其光学性能的关键因素。通过理解和控制这些衍射效应，可以设计和优化三维超材料的光学器件，实现各种光学功能，推动光学技术的发展。随着三维超材料研究的不断深入，其在光学领域的应用潜力将进一步得到挖掘和拓展。
📣 部分代码
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Inputs and Structure Parametersy = 500; % Range of wavelength in nm% Anglesth = 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 parametersmlt = 1; % Scaling factor for the hexagonal unit cellD_x = 300; % Length of one period in 'x'D_x = 2*round(D_x*mlt/2); % Scaled and rounded to nearest even numberD_y = 2*round(sqrt(3)*D_x/2); % Length of one period in 'y'% Metal-dielectric metamaterialsD_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 StructuresN_m = 4; % =0; for uncarved % Number of elementary cells carved into pillarsr_1 = 50*mlt; % Radius of pillars at topr_2 = 145*mlt; % Radius of pillars at bottomN_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 cellh_zd = D_d/N_zd; % Resolution in 'z' for carved dielectric layersh_z = repmat([repmat(h_zd,1,N_zd),D_m,repmat(h_zd,1,N_zd)],1,N_m); % Height of each carved layerz_l = [0,cumsum(h_z(1:N_z-1))]; % 'z' coordinate of each carved layerr_z = interp1([0,D_z-h_zd],[r_1,r_2],z_l); % Radius of pillars for each carved layerN_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 directionN_y = 9; %=1; for uncarved % Number of Fourier modes in y directionNt = N_x*N_y;% Preallocationsky = 2*pi./y; % WavenumbersI = eye(Nt); II = eye(2*Nt); % Identity matricesZZ = 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 arraysmx = (-(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.infoMe = interp1(Cu(:,1),Cu,y/1000); % Interpolation of experimental dataMetal = (Me(:,2) + 1i*Me(:,3)); Metal = Metal.';kkk = D_m/39; % Thickness of W layer divided mean free path of Wkky = 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^-1gamma = 12*10^-2; gammakk = gamma*fi/kkk;epsDrude = 1-(wp^2)./((kky/(2*pi)).^2+1i*(kky/(2*pi))*gamma); % Drude model for the BulkepsDrudekk = 1-(wp^2)./((kky/(2*pi)).^2+1i*(kky/(2*pi))*gammakk); % Drude Model corrected for the thin layerepsMetal = Metal-epsDrude+epsDrudekk;%% Index profile of moth eye structures% Build the moth-eye structure as a binary 3D matrixindexProfLog = false(D_x,D_y,N_z);% Hexagon vertice positionscentre = [D_x,D_y]/2;for nn=1:N_zr = 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 pillarsindexProfLog((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 Analysisuy = y/1000; % Wavelength in micron% Permittivity and refractive indicesn_1 = 1; % Refractive index of incidence mediumn_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 wavelengthe_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 pillarse_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 layerindexProf = ones(D_x,D_y);for nn=1:N_zindexProf(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_xfor qq = 1:N_yEE = 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,[]);endendi_iepsx_mj = zeros(N_x,D_y,N_x);iepsx_fft = fftshift(fft(indexProf.^(-1),[],1),1)/D_x;for qq=1:D_yiepsx_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);endepsxy_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_xfor qq = 1:N_xE4(pp,:,qq,:) = toeplitz(epsxy_mnj(pp,N_y:2*N_y-1,qq),flip(epsxy_mnj(pp,1:N_y,qq)));endendEPSxy(:,:,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_xiepsy_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);endepsyx_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_yfor ss = 1:N_yE4(:,rr,:,ss) = toeplitz(epsyx_mnl(N_x:2*N_x-1,rr,ss),flip(epsyx_mnl(1:N_x,rr,ss)));endendEPSyx(:,:,nn) = reshape(E4,[N_x*N_y,N_x*N_y]);end%% Initializing variables% Incident Fieldu_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 Matricesk_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);%% Eigenmatricesfor nn=1:N_z% Eigenvalue formulation for the grating regionF11 = -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 layersW(:,:,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 layersfor nn=1:3*NN_zif (mod(nn+1,3))==0k_n_zmn = k_m_zmn;D_n = D_m;EPS(:,:,1+N_z+nn) = conj(e_m)*I;elsek_n_zmn = k_d_zmn;D_n = D_d;EPS(:,:,1+N_z+nn) = conj(e_d)*I;endW(:,:,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% PreallocationB = repmat(ZZ,2*(N_l+1),2*(N_l+1));% Building B using bfor ll=1:N_lb = [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径向基神经网络时序、回归预测和分类
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
