电磁波 基于matlab 三维FDTD方法电磁波传播仿真【含Matlab源码 3870期】.zip

%微带线的FDTD仿真 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%最基本设定 TimeStop=2000; %计算时间步 f0=6.7e9; %设定中心频点 w0=2*pi*f0; %中心角频率 delta_t=3.3083e-013; W0=w0*delta_t; %数字角频率 f_max=10e9; f_min=2e9; Gaussian_T=0.5/f_max; Gaussian_t0=3*Gaussian_T; Patch_W=15e-3; %物理尺寸的设置(天线尺寸) Patch_H=0.794e-3; Dielectric_er=2.32; %基板相对介质常数 c=3*10^8; Mu0=pi*4e-7; Eps0=1/c^2/Mu0; Eps1=Eps0*Dielectric_er; Eps10=(Eps0+Eps1)/2; Total_L=400; %设定计算领域 Total_W=50; Total_H=9; Antenna_W=40; %天线的格子数 Antenna_H=4; delta_y=Patch_W/Antenna_W; %网格长度 delta_x=delta_y; delta_z=Patch_H/Antenna_H; Pml_L=8; %吸收边界参数设定 R0=0.00001; n=2; Ipmin=1;Imin=Ipmin+Pml_L;Imax=Imin+Total_L;Ipmax=Imax+Pml_L; %吸收边界 Jpmin=1;Jmin=Jpmin+Pml_L;Jmax=Total_W+Jmin;Jpmax=Jmax+Pml_L; Kmin=1;Kmax=Total_H+Kmin;Kpmax=Kmax+Pml_L; Iamin=Imin+2; Jamin=Jmin+round((Total_W-Antenna_W)/2);Jamax=Antenna_W+Jamin; Kamax=1+Antenna_H; J_ref0=Jamin+1; I_ref0=Iamin+3; I_ref1=Imax-3; I_ref2=Imax+3; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%在计算微带线激励时要用到的场 Ez(Ipmin:Ipmax,Jpmin:Jpmax,Kmin:(Kpmax-1))=0; %FDTD采样点Ez，Ez(i,j,k+1/2,t=n) Ex(Ipmin:(Ipmax-1),Jpmin:Jpmax,Kmin:Kpmax)=0; %FDTD采样点Ex，Ex(i+1/2,j,k,t=n) Ey(Ipmin:Ipmax,Jpmin:(Jpmax-1),Kmin:Kpmax)=0; %FDTD采样点Ey，Ey(i,j+1/2,k,t=n) Hz(Ipmin:(Ipmax-1),Jpmin:(Jpmax-1),Kmin:Kpmax)=0; %FDTD采样点Hz，Hz(i+1/2,j+1/2,k,t=n+1/2) Hx(Ipmin:Ipmax,Jpmin:(Jpmax-1),Kmin:(Kpmax-1))=0; %FDTD采样点Hx，Ex(i,j+1/2,k+1/2,t=n+1/2) Hy(Ipmin:(Ipmax-1),Jpmin:Jpmax,Kmin:(Kpmax-1))=0; %FDTD采样点Hy，Hy(i+1/2,j,k+1/2,t=n+1/2) Ezx(Ipmin:Ipmax,Jpmin:Jpmax,Kmin:(Kpmax-1))=0; %FDTD采样点Ez，Ez(i,j,k+1/2,t=n) Ezy(Ipmin:Ipmax,Jpmin:Jpmax,Kmin:(Kpmax-1))=0; %FDTD采样点Ez，Ez(i,j,k+1/2,t=n) Exy(Ipmin:(Ipmax-1),Jpmin:Jpmax,Kmin:Kpmax)=0; %FDTD采样点Ex，Ex(i+1/2,j,k,t=n) Exz(Ipmin:(Ipmax-1),Jpmin:Jpmax,Kmin:Kpmax)=0; %FDTD采样点Ex，Ex(i+1/2,j,k,t=n) Eyx(Ipmin:Ipmax,Jpmin:(Jpmax-1),Kmin:Kpmax)=0; %FDTD采样点Ey，Ey(i,j+1/2,k,t=n) Eyz(Ipmin:Ipmax,Jpmin:(Jpmax-1),Kmin:Kpmax)=0; %FDTD采样点Ey，Ey(i,j+1/2,k,t=n) Hzx(Ipmin:(Ipmax-1),Jpmin:(Jpmax-1),Kmin:Kpmax)=0; %FDTD采样点Hz，Hz(i+1/2,j+1/2,k,t=n+1/2) Hzy(Ipmin:(Ipmax-1),Jpmin:(Jpmax-1),Kmin:Kpmax)=0; %FDTD采样点Hz，Hz(i+1/2,j+1/2,k,t=n+1/2) Hxy(Ipmin:Ipmax,Jpmin:(Jpmax-1),Kmin:(Kpmax-1))=0; %FDTD采样点Hx，Ex(i,j+1/2,k+1/2,t=n+1/2) Hxz(Ipmin:Ipmax,Jpmin:(Jpmax-1),Kmin:(Kpmax-1))=0; %FDTD采样点Hx，Ex(i,j+1/2,k+1/2,t=n+1/2) Hyx(Ipmin:(Ipmax-1),Jpmin:Jpmax,Kmin:(Kpmax-1))=0; %FDTD采样点Hy，Hy(i+1/2,j,k+1/2,t=n+1/2) Hyz(Ipmin:(Ipmax-1),Jpmin:Jpmax,Kmin:(Kpmax-1))=0; %FDTD采样点Hy，Hy(i+1/2,j,k+1/2,t=n+1/2) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%在非PML区FDTD迭代要用到的系数 MediaNumber=2; %表示只有两种不同的介质：空气，Duroid %epsilon_r Media(1,1)=1; %epsilon_r Media(2,1)=1; %mu_r Media(3,1)=0; %sigma Media(4,1)=0; %sigma_m Media(1,2)=Dielectric_er; Media(2,2)=1; Media(3,2)=0; Media(4,2)=0; Media(1,1:MediaNumber)=Media(1,1:MediaNumber)*Eps0; Media(2,1:MediaNumber)=Media(2,1:MediaNumber)*Mu0; for I=1:MediaNumber for J=1:MediaNumber CA1(I,J)=((Media(1,I)+Media(1,J))-0.5*(Media(3,I)+Media(3,J))*delta_t)/((Media(1,I)+Media(1,J))+0.5*(Media(3,I)+Media(3,J))*delta_t); %电场计算公式第一项系数 CB1(I,J)=2*delta_t/((Media(1,I)+Media(1,J))+0.5*(Media(3,I)+Media(3,J))*delta_t); %电场计算公式第二项系数 CP1(I,J)=((Media(2,I)+Media(2,J))-0.5*(Media(4,I)+Media(4,J))*delta_t)/((Media(2,I)+Media(2,J))+0.5*(Media(4,I)+Media(4,J))*delta_t); %磁场场计算公式第一项系数 CQ1(I,J)=2*delta_t/((Media(2,I)+Media(2,J))+0.5*(Media(4,I)+Media(4,J))*delta_t); %磁场场计算公式第二项系数 end end CA(Kmin:(Kamax-1))=CA1(2,2);CB(Kmin:(Kamax-1))=CB1(2,2);CP(Kmin:(Kamax-1))=CP1(2,2);CQ(Kmin:(Kamax-1))=CQ1(2,2); CA((Kamax+1):Kmax)=CA1(1,1);CB((Kamax+1):Kmax)=CB1(1,1);CP((Kamax+1):Kmax)=CP1(1,1);CQ((Kamax+1):Kmax)=CQ1(1,1); CA(Kamax)=CA1(1,2);CB(Kamax)=CB1(1,2);CP(Kamax)=CP1(1,2);CP(Kamax)=CP1(1,2);CQ(Kamax)=CQ1(1,2); CA_half(Kmin:(Kamax-1))=CA1(2,2);CB_half(Kmin:(Kamax-1))=CB1(2,2);CP_half(Kmin:(Kamax-1))=CP1(2,2);CQ_half(Kmin:(Kamax-1))=CQ1(2,2); CA_half(Kamax:(Kmax-1))=CA1(1,1);CB_half(Kamax:(Kmax-1))=CB1(1,1);CP_half(Kamax:(Kmax-1))=CP1(1,1);CQ_half(Kamax:(Kmax-1))=CQ1(1,1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%在进行PML计算时要用到的电导率和磁导率的设定 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%sigma_x L=0; sigma_x0(1+L*2)=Eps0*c*(-log(R0))/(2^(n+2)*delta_x*Pml_L^(n+1)); %sigma_x(1,2,3,,...)=sigma_x(0,1/2,1,...) for L=1/2:1/2:Pml_L sigma_x0(1+L*2)=sigma_x0(1)*((2*L+1)^(n+1)-(2*L-1)^(n+1)); %可从Perfectly matched layer for the FDTD solution of wave-structure interaction problems.pdf查到相应公式 end sigma_m_x0=sigma_x0/Eps0*Mu0; L=0; sigma_x1(1+L*2)=Eps1*c*(-log(R0))/(2^(n+2)*delta_x*Pml_L^(n+1)); for L=1/2:1/2:Pml_L sigma_x1(1+L*2)=sigma_x1(1)*((2*L+1)^(n+1)-(2*L-1)^(n+1)); end sigma_m_x1=sigma_x1/Eps1*Mu0; L=0; sigma_x10(1+L*2)=Eps10*c*(-log(R0))/(2^(n+2)*delta_x*Pml_L^(n+1)); for L=1/2:1/2:Pml_L sigma_x10(1+L*2)=sigma_x10(1)*((2*L+1)^(n+1)-(2*L-1)^(n+1)); end sigma_m_x10=sigma_x10/Eps10*Mu0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%sigma_y L=0; sigma_y0(1+L*2)=Eps0*c*(-log(R0))/(2^(n+2)*delta_y*Pml_L^(n+1)); %sigma_y(1,2,3,,...)=sigma_y(0,1/2,1,...) for L=1/2:1/2:Pml_L sigma_y0(1+L*2)=sigma_y0(1)*((2*L+1)^(n+1)-(2*L-1)^(n+1)); end sigma_m_y0=sigma_y0/Eps0*Mu0; L=0; sigma_y1(1+L*2)=Eps1*c*(-log(R0))/(2^(n+2)*delta_y*Pml_L^(n+1)); for L=1/2:1/2:Pml_L sigma_y1(1+L*2)=sigma_y1(1)*((2*L+1)^(n+1)-(2*L-1)^(n+1)); end sigma_m_y1=sigma_y1/Eps1*Mu0; L=0; sigma_y10(1+L*2)=Eps10*c*(-log(R0))/(2^(n+2)*delta_y*Pml_L^(n+1)); for L=1/2:1/2:Pml_L sigma_y10(1+L*2)=sigma_y10(1)*((2*L+1)^(n+1)-(2*