Tot_Src_ball_Mur.zip_FDTD 二维 matlab_Oblique incidence_PML_二维FDTD

clear; IsFigure=1; %画图=1 WantToSeeEp=1; %看电磁常数分布=1 %定义 基本常数 mu0=4*pi*1e-7; e0=8.85*1e-12; c=1/sqrt(mu0*e0); factor_te=mu0/e0; factor_tm=e0/mu0; Zac=sqrt(factor_te); epr=1; %基底介电常数 mur=1; ep=epr*e0; mu=mur*mu0; a=1e-6; %晶格常数。 Zmax=0.6; %The maximum value for z axis when plotting figures. Colormax=0.6; %The maximum value for colormap when plotting figures. %定义 频率 W1=1; %归一化频率——能带图中的纵坐标，单位为（a/lamda） W=W1*(2*pi*c/a); %波源cos(Wt),W为角速度，频率为W/(2*pi) v=c/sqrt(epr); %%%%%%%%%%定义时间循环次数 （总时长为： NTimeSteps * dt） NTimeSteps=300; %循环总数 Meach=2; %每隔Meach个时间步显示一次 %%%%%%%%%% %定义 空间网格数量 MLatx=11; % x方向晶胞个数 MLaty=11; % y方向 NMlat=21; %一个晶胞中包含的基本网格数 %必须为奇数,目的是为了将晶胞中心点位于网格点上 if mod(NMlat,2)==0 NMlat=NMlat+1; end %强迫NMlat为奇数！ NTx=MLatx*NMlat; %总基本网格数! NTy=MLaty*NMlat; Lc=40; %散射场宽度 cita=pi/6; %定义 网格离散长度 dx=a/NMlat; %定义FDTD基本网格大小——决定精确度 dy=dx; dt=1/sqrt(1/(dx*dx)+1/(dy*dy))/c; %时间长度——决定稳定性 %定义 几个用于化简计算量的中间参数 c1=(v*dt-dx)/(v*dt+dx); c2=v*dt+dx; c3=(v*dt-sqrt(2)*dx)/(v*dt+sqrt(2)*dx); %定义 含/不含有波导结构时 电磁参数-空间分布 Ep=ones(NTx,NTy)*ep; Mu=ones(NTx,NTy)*mu; %是否显示空间介电常数分布 if WantToSeeEp==1 x=0:dx:(NTx-1)*dx; x=x-NTx*dx/2; y=0:dy:(NTy-1)*dy; y=y-NTy*dy/2; [X,Y]=meshgrid(x,y); X=X'; Y=Y'; surf(X,Y,Ep(1:NTx,1:NTy)/e0); shading interp; view(0,90); axis([min(x), max(x),min(y), max(y)]) axis equal; axis off; disp('Press any key to continue...'); pause end %各个电磁场量——初始赋值 Ez=zeros(NTx,NTy); %此处认为的矩阵（i,j）与坐标（x,y）相反,从而在空间的形式上一致,方便想象 Hx=zeros(NTx,NTy-1); Hy=zeros(NTx-1,NTy); %开始离散Maxwell方程的时间大循环： tic for m=1:NTimeSteps %中心区——X、Y分量的离散方程 Hx=Hx-2*dt./(Mu(:,1:NTy-1)+Mu(:,2:NTy))./dy.*(Ez(:,2:NTy)-Ez(:,1:NTy-1)); Hy=Hy+2*dt./(Mu(1:NTx-1,:)+Mu(2:NTx,:))./dx.*(Ez(2:NTx,:)-Ez(1:NTx-1,:)); %定义Mur贰阶边界条件:（部分dx、dy混用了，均匀FDTD网格时暂不予考虑） %为Mur吸收边界保留处前一时刻的场值（m-1时刻的 Ez、Hz） 所以插在XY与Z分量中间 Ezl=Ez(2,2:NTy-1); Ezr=Ez(NTx-1,2:NTy-1); Ezu=Ez(2:NTx-1,NTy-1); Ezd=Ez(2:NTx-1,2); %%%%%%%%%%%%% 去激励源E 总场E-散射场H： %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %算左侧Hy→E for j=Lc+1:NTy-Lc i=Lc; T0=((i+0.5)*cos(cita)*dx+j*sin(cita)*dy)*epr^0.5/c; if T0>=m*dt ts=0; else ts=m*dt-T0; end Ez(i,j)=Ez(i,j)-dt/dx/Ep(i,j)/Zac*cos(cita)*sin(ts*W); end % %算右侧Hy→E for j=Lc+1:NTy-Lc i=NTx-Lc+1; T0=((i+0.5)*cos(cita)*dx+j*sin(cita)*dy)*epr^0.5/c; if T0>=m*dt ts=0; else ts=m*dt-T0; end Ez(i,j)=Ez(i,j)+dt/dx/Ep(i,j)/Zac*cos(cita)*sin(ts*W); end % %算下侧Hx→E for i=Lc+1:NTx-Lc j=Lc; T0=(i*cos(cita)*dx+(j+0.5)*sin(cita)*dy)*epr^0.5/c; if T0>=m*dt ts=0; else ts=m*dt-T0; end Ez(i,j)=Ez(i,j)-dt/dx/Ep(i,j)/Zac*sin(cita)*sin(ts*W); end %算上侧Hx→E for i=Lc+1:NTx-Lc j=NTy-Lc+1; T0=(i*cos(cita)*dx+(j+0.5)*sin(cita)*dy)*epr^0.5/c; if T0>=m*dt ts=0; else ts=m*dt-T0; end Ez(i,j)=Ez(i,j)+dt/dx/Ep(i,j)/Zac*sin(cita)*sin(ts*W); end %%%%%四顶角 T1l=((i+0.5)*cos(cita)*dx+j*sin(cita)*dy)*epr^0.5/c; T1d=(i*cos(cita)*dx+(j+0.5)*sin(cita)*dy)*epr^0.5/c; if T1l>=m*dt t1l=0; else t1l=m*dt-T1l; end if T1d>=m*dt t1d=0; else t1d=m*dt-T1d; end Ez(Lc,Lc)=Ez(Lc,Lc)-dt/dx/Ep(Lc,Lc)/Zac*cos(cita)*sin(t1l*W)-dt/dx/Ep(Lc,Lc)/Zac*sin(cita)*sin(t1d*W); %%%%%%%%% T2l=((i+0.5)*cos(cita)*dx+j*sin(cita)*dy)*epr^0.5/c; T2u=(i*cos(cita)*dx+(j+0.5)*sin(cita)*dy)*epr^0.5/c; if T2l>=m*dt t2l=0; else t2l=m*dt-T2l; end if T2u>=m*dt t2u=0; else t2u=m*dt-T2u; end Ez(Lc,NTy-Lc+1)=Ez(Lc,NTy-Lc+1)-dt/dx/Ep(Lc,NTy-Lc+1)/Zac*cos(cita)*sin(t2l*W)-dt/dx/Ep(Lc,NTy-Lc+1)/Zac*sin(cita)*sin(t2u*W); %%%%%%%%% T3r=((i+0.5)*cos(cita)*dx+j*sin(cita)*dy)*epr^0.5/c; T3d=(i*cos(cita)*dx+(j+0.5)*sin(cita)*dy)*epr^0.5/c; if T3r>=m*dt t3r=0; else t3r=m*dt-T3r; end if T3d>=m*dt t3d=0; else t3d=m*dt-T3d; end Ez(NTx-Lc+1,Lc)=Ez(NTx-Lc+1,Lc)+dt/dx/Ep(NTx-Lc+1,Lc)/Zac*cos(cita)*sin(t3r*W)+dt/dx/Ep(NTx-Lc+1,Lc)/Zac*sin(cita)*sin(t3d*W); %%%%%%%%% T4r=((i+0.5)*cos(cita)*dx+j*sin(cita)*dy)*epr^0.5/c; T4u=(i*cos(cita)*dx+(j+0.5)*sin(cita)*dy)*epr^0.5/c; if T4r>=m*dt t4r=0; else t4r=m*dt-T4r; end if T4u>=m*dt t4u=0; else t4u=m*dt-T4u; end Ez(NTx-Lc+1,NTy-Lc+1)=Ez(NTx-Lc+1,NTy-Lc+1)+dt/dx/Ep(NTx-Lc+1,NTy-Lc+1)/Zac*cos(cita)*sin(t4r*W)+dt/dx/Ep(NTx-Lc+1,NTy-Lc+1)/Zac*sin(cita)*sin(t4u*W); %%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%% %中心区——Z分量的离散方程：（注意！！→计算时选取的矩阵在边界以内一格，边界处场值由后面程序去处理） Ez(2:NTx-1,2:NTy-1)=Ez(2:NTx-1,2:NTy-1)+dt./Ep(2:NTx-1,2:NTy-1).*(... (Hy(2:NTx-1,2:NTy-1)-Hy(1:NTx-2,2:NTy-1))./dx-(Hx(2:NTx-1,2:NTy-1)-Hx(2:NTx-1,1:NTy-2))./dy); %Mur吸收边界 计算处理 %左边界： Ez(1,2:NTy-1)=Ezl+c1.*(Ez(2,2:NTy-1)-Ez(1,2:NTy-1))... -Mu(1,2:NTy-1).*v^2*dt*dx/(2*dy*c2).*(Hx(1,2:NTy-1)-Hx(1,1:NTy-2)+Hx(2,2:NTy-1)-Hx(2,1:NTy-2)); %右边界： Ez(NTx,2:NTy-1)=Ezr+c1.*(Ez(NTx-1,2:NTy-1)-Ez(NTx,2:NTy-1))... -Mu(NTx,2:NTy-1).*v^2*dt*dx/(2*dy*c2).*(Hx(NTx,2:NTy-1)-Hx(NTx,1:NTy-2)+Hx(NTx-1,2:NTy-1)-Hx(NTx-1,1:NTy-2)); %下边界： Ez(2:NTx-1,1)=Ezd+c1.*(Ez(2:NTx-1,2)-Ez(2:NTx-1,1))... +Mu(2:NTx-1,1).*v^2*dt*dx/(2*dy*c2).*(Hy(2:NTx-1,1)-Hy(1:NTx-2,1)+Hy(2:NTx-1,2)-Hy(1:NTx-2,2)); %上边界： Ez(2:NTx-1,NTy)=Ezu+c1.*(Ez(2:NTx-1,NTy-1)-Ez(2:NTx-1,NTy))... +Mu(2:NTx-1,NTy).*v^2*dt*dx/(2*dy*c2).*(Hy(2:NTx-1,NTy)-Hy(1:NTx-2,NTy)+Hy(2:NTx-1,NTy-1)-Hy(1:NTx-2,NTy-1)); %左下角： Ez(1,1)=Ezl(1)+c3.*(Ez(2,2)-Ez(1,1)); %右上角： Ez(NTx,NTy)=Ezr(NTy-2)+c3.*(Ez(NTx-1,NTy-1)-Ez(NTx,NTy)); %左上角： Ez(1,NTy)=Ezu(1)+c3.*(Ez(2,NTy-1)-Ez(1,NTy)); %右下角： Ez(NTx,1)=Ezr(1)+c3.*(Ez(NTx-1,2)-Ez(NTx,1)); %%%%%%%%%%%%% 激励源E 总场E-散射场H： %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %算左侧Hy→E for j=Lc+1:NTy-Lc i=Lc; T0=(i*cos(cita)*dx+j*sin(cita)*dy)*epr^0.5/c; if T0>=m*dt ts=0; else ts=m*dt-T0; end Hy(i-1,j)=Hy(i-1,j)+dt/dx/Mu(i,j)*sin(ts*W); end % %算右侧Hy→E for j=Lc+1:NTy-Lc i=NTx-Lc+1; T0=(i*cos(cita)*dx+j*sin(cita)*dy)*epr^0.5/c; if T0>=m*dt ts=0; else ts=m*dt-T0; end Hy(i,j)=Hy(i,j)-dt/dx/Mu(i,j)*sin(ts*W); end % %算下侧Hx→E for i=Lc+1:NTx-Lc j=Lc-1; T0=(i*cos(cita)*dx+j*sin(cita)*dy)*epr^0.5/c; if T0>=m*dt ts=0; else ts=m*dt-T0; end Hx(i-1,j)=Hx(i-1,j)-dt/dx/Mu(i,j)*sin(ts*W); end %算