带挡板的二维平行板传输线FDTD仿真，使用高斯激励

clc;clear; c=3.0e8; muz=4.0*pi*1.0e-7; epsz=1.0/(c*c*muz); T=0.5e-9; t0=3*T; freq=1.0e+9; lambda=c/freq; L=6; d=0.18; dz=lambda/20; dx=lambda/20; dt=dz/(2*c); Nz=L/dz; Nx=d/dx; nmax=1000; z=linspace(0,Nz*dz,Nz+1); x=linspace(0,Nx*dx,Nx+1); t=linspace(0,nmax*dt,nmax); ex=zeros(Nx,Nz+1); ez=zeros(Nx+1,Nz); hy=zeros(Nx,Nz); ex_air=zeros(Nx,Nz+1); ez_air=zeros(Nx+1,Nz); hy_air=zeros(Nx,Nz); etotal=zeros(Nx,nmax); etr=zeros(Nx,nmax); ein=zeros(Nx,nmax); ex1=zeros(Nx,1);ex2=zeros(Nx,1); ex3=zeros(Nx,1);ex4=zeros(Nx,1); for n=1:nmax if n<120 ex(1:Nx,1)=exp(-((n-1)*dt-t0)^2/(T^2)); ex_air(1:Nx,1)=exp(-((n-1)*dt-t0)^2/(T^2)); end ex1_2=ex1;ex2_2=ex2; ex3_2=ex3;ex4_2=ex4; ex1=ex(1:Nx,1);ex2=ex(1:Nx,2); ex3=ex(1:Nx,Nz+1);ex4=ex(1:Nx,Nz); ex1_air=ex_air(1:Nx,1);ex2_air=ex_air(1:Nx,2); ex3_air=ex_air(1:Nx,Nz+1);ex4_air=ex_air(1:Nx,Nz); hy(1:Nx,1:Nz)=hy(1:Nx,1:Nz)+(dt./muz).*((ez(2:Nx+1,1:Nz)-ez(1:Nx,1:Nz))./dx+(ex(1:Nx,1:Nz)-ex(1:Nx,2:Nz+1))./dz); ex(1:Nx,2:Nz)=ex(1:Nx,2:Nz)+(dt./epsz).*((hy(1:Nx,1:Nz-1)-hy(1:Nx,2:Nz))./dz); ez(2:Nx,1:Nz)=ez(2:Nx,1:Nz)+(dt./epsz).*((hy(2:Nx,1:Nz)-hy(1:Nx-1,1:Nz))./dx); hy_air(1:Nx,1:Nz)=hy_air(1:Nx,1:Nz)+(dt./muz).*((ez_air(2:Nx+1,1:Nz)-ez_air(1:Nx,1:Nz))./dx+(ex_air(1:Nx,1:Nz)-ex_air(1:Nx,2:Nz+1))./dz); ex_air(1:Nx,2:Nz)=ex_air(1:Nx,2:Nz)+(dt./epsz).*((hy_air(1:Nx,1:Nz-1)-hy_air(1:Nx,2:Nz))./dz); ez_air(2:Nx,1:Nz)=ez_air(2:Nx,1:Nz)+(dt./epsz).*((hy_air(2:Nx,1:Nz)-hy_air(1:Nx-1,1:Nz))./dx); ex(1:fix(Nx/3),Nz/2+1)=0; ex(fix(2*Nx/3):Nx,Nz/2+1)=0; ex(2:Nx-1,1)=-ex2_2(2:Nx-1)+((c*dt-dz)/(c*dt+dz)).*(ex(2:Nx-1,2)+ex1_2(2:Nx-1))+(2*dz/(c*dt+dz)).*(ex2(2:Nx-1)+ex1(2:Nx-1))+... ((c*dt)^2*dz/(2*dx^2*(c*dt+dz))).*(ex1(3:Nx)-2*ex1(2:Nx-1)+ex1(1:Nx-2)+ex2(3:Nx)-2*ex2(2:Nx-1)+ex2(1:Nx-2)); ex(2:Nx-1,Nz+1)=-ex4_2(2:Nx-1)+((c*dt-dz)/(c*dt+dz))*(ex(2:Nx-1,Nz)+ex3_2(2:Nx-1))+(2*dz/(c*dt+dz)).*(ex4(2:Nx-1)+ex3(2:Nx-1))+... ((c*dt)^2*dz/(2*dx^2*(c*dt+dz))).*(ex3(3:Nx)-2*ex3(2:Nx-1)+ex3(1:Nx-2)+ex4(3:Nx)-2*ex4(2:Nx-1)+ex4(1:Nx-2)); ex(1,1)=ex2(1)+((c*dt-dz)/(c*dt+dz)).*(ex(1,2)-ex1(1)); ex(Nx,1)=ex2(Nx)+((c*dt-dz)/(c*dt+dz)).*(ex(Nx,2)-ex1(Nx)); ex(1,Nz+1)=ex4(1)+((c*dt-dz)/(c*dt+dz)).*(ex(1,Nz)-ex3(1)); ex(Nx,Nz+1)=ex4(Nx)+((c*dt-dz)/(c*dt+dz)).*(ex(Nx,Nz)-ex3(Nx)); ez(1,1:Nz)=0; ez(Nx+1,1:Nz)=0; % ex(1:Nx,1)=ex2+((c*dt-dz)/(c*dt+dz)).*(ex(1:Nx,2)-ex1); % ex(1:Nx,Nz+1)=ex4+((c*dt-dz)/(c*dt+dz)).*(ex(1:Nx,Nz)-ex3); % ez(1,1:Nz)=0; % ez(Nx+1,1:Nz)=0; etotal(1:Nx,n)=ex(1:Nx,10); etr(1:Nx,n)=ex(1:Nx,390); ein(1:Nx,n)=ex_air(1:Nx,10); if mod(n,10)==0 rtime=num2str(round(n*dt/1.0e-9)); figure(1); subplot(2,1,1); mesh(ex); axis([0 400 0 12 -1 1]); caxis([-1 1]); title(['time = ',rtime,' ns']); subplot(2,1,2) imagesc(ex); caxis([-1 1]); pause(0.01) end end ere=etotal-ein; vin=sum(dx.*ein); vre=sum(dx.*ere); vtr=sum(dx.*etr); figure(2);plot(t,vin);axis([0 nmax*dt -0.2 0.2]); title('入射电压');xlabel('t (s)');ylabel('V'); figure(3);plot(t,vre);axis([0 nmax*dt -0.2 0.2]); title('反射电压');xlabel('t (s)');ylabel('V'); figure(4);plot(t,vtr);axis([0 nmax*dt -0.2 0.2]); title('透射电压');xlabel('t (s)');ylabel('V'); M=0:nmax-1; f=M/(nmax*dt); Vin=fft(vin); Vre=fft(vre); Vtr=fft(vtr); figure(5);plot(f,abs(Vin));axis([0 1.0e9 0 10]); title('入射电压');xlabel('f (Hz)');ylabel('|Vre|'); figure(6);plot(f,abs(Vre));axis([0 1.0e9 0 10]); title('反射电压');xlabel('f (Hz)');ylabel('|Vtr|'); figure(7);plot(f,abs(Vtr));axis([0 1.0e9 0 10]); title('透射电压');xlabel('f (Hz)');ylabel('|Vin|'); s11=Vre./Vin; s21=Vtr./Vin; figure(8); subplot(3,1,1),plot(f,abs(s11));axis([0 1.0e9 0 2]);xlabel('f (Hz)');ylabel('S11'); subplot(3,1,2),plot(f,abs(s21));axis([0 1.0e9 0 2]);xlabel('f (Hz)');ylabel('S21'); figure(9); plot(f,abs(s11).^2+abs(s21).^2);axis([0 1.0e9 0 2]);xlabel('f (Hz)');