总场散射场FDTD

% Total/scattered Field Formulation % plane wave in free space% The Perfectly Matched Layer% FDTD in Two Dimensions 3_2_3dKE=60;IE=KE;JE=KE;kc=KE/2;T=0;NSTEP=200;% Initialize arrays % dz=zeros(KE,KE);ez=zeros(KE,KE);hx=zeros(KE,KE);hy=zeros(KE,KE);ihx=zeros(KE,KE);ihy=zeros(KE,KE); ga=ones(KE,KE);ez_inc=zeros(1,JE);hx_inc=zeros(1,JE);% the PML parameters % gi2=ones(1,IE);gi3=gi2;fi1=zeros(1,IE);fi2=ones(1,IE);fi3=fi2;gj2=ones(1,JE);gj3=gj2;fj1=zeros(1,JE); fj2=ones(1,JE);fj3=fj2;npml=8;for i=1:8 xnum=npml-i; xd=npml; xxn=xnum/xd;xn = 0.33*xxn3;gi2(i)=1.0/(1.0+xn);gi2(IE-1-i)=1.0/(1.0+xn);gi3(i)=(1.0-xn)/(1.0+xn);gi3(IE-1-i)=(1.0-xn)/(1.0+xn);xxn=(xnum-0.5)/xd;xn = 0.25*xxn3;fi1(i)=xn;fi1(IE-2-i)=xn; fi2(i)=1.0/(1.0+xn);fi2(IE-2-i)=1.0/(1.0+xn);fi3(i)=(1.0-xn)/(1.0+xn);fi3(IE-2-i)=(1.0-xn)/(1.0+xn);endfor j=1:8xnum=npml-j;xd=npml;xxn=xnum/xd;xn = 0.33*xxn3;gj2(j)=1.0/(1.0+xn);gj2(JE-1-j)=1.0/(1.0+xn);gj3(j)=(1.0-xn)/(1.0+xn);gj3(JE-1-j)=(1.0-xn)/(1.0+xn);xxn=(xnum-0.5)/xd;xn = 0.25*xxn3;fj1(j)=xn;fj1(JE-2-j)=xn;fj2(j)=1.0/(1.0+xn);fj2(JE-2-j)=1.0/(1.0+xn);fj3(j)=(1-xn)/(1.0+xn);fj3(JE-2-j)=(1.0-xn)/(1.0+xn);end% total/scattered field boundaries % ia=7;ib=IE-ia-1;ja=7;jb=JE-ja-1;ez_inc_low_m1=0;ez_inc_low_m2=0;ez_inc_high_m1=0;ez_inc_high_m2=0;% sigal source %t0=40;spread=12;% grid %ddx=0.01; % set the cell size to 1cm % dt=ddx/(2*3e8); % calculate the time step % epsz=8.851*1e-12;% material parameter % epsilon=4; % relative dielectric constant % kstart=100; % position of the dielectric medium % epsz=8.85419e-12; % absolute dielectric constant % sigma=0.04; % conductivity %for n=1:NSTEPa=1;T=T+1;for j=2:JEez_inc(j)= ez_inc(j)+0.5*( hx_inc(j-1)- hx_inc(j);end% ABC for the incident buffer % ez_inc(1) =ez_inc_low_m2;ez_inc_low_m2=ez_inc_low_m1; ez_inc_low_m1=ez_inc(2);ez_inc(JE-1) =ez_inc_high_m2; ez_inc_high_m2=ez_inc_high_m1;ez_inc_high_m1=ez_inc(JE-2);% calculate dx field %for i=2:KEfor j=2:KEdz(i,j)=gi3(i)*gj3(j)*dz(i,j)+gi2(i)*gj2(j)*0.5*(hy(i,j)-hy(i-1,j)-hx(i,j)+hx(i,j-1);endend% insert source %pulse=exp(-0.5*(t0-T)人2/sp read2);%sin(2*pi*1500*1e6*dt*T);ez_inc(3)=pulse;% Incident DZ values %for i=ia:ibdz(i,ja)= dz(i,ja)+0.5*hx_inc(ja-1); dz(i,jb)= dz(i,jb)-0.5*hx_inc(jb);end% calculate E from dx %for i=1:KEfor j=1:KEez(i,j)=ga(i,j)*dz(i,j);endend% set the EZ edges to 0, as part of the PML % for j=1:JE-1ez(1,j)=0.0; ez(IE,j)=0.0;endfor i=1:IE-1ez(i,1)=0.0; ez(i,JE)=0.0;endfor j=1:JE-1hx_inc(j)=hx_inc(j)+0.5*( ez_inc(j)-ez_inc(j+1);end% calculate hx field %for i=1:KEfor j=1:KE-1 curl_e=ez(i,j)-ez(i,j+1);ihx(i,j)=ihx(i,j)+fi1(i)*curl_e; hx(i,j)=fj3(j)*hx(i,j)+fj2(j)*0.5* (curl_e+ihx(i,j);endend% Incident Hx vzlues %for i=ia:ibhx(i,ja-1)= hx(i,ja-1)+0.5*ez_inc(ja); hx(i,jb)= hx(i,jb)-0.5*ez_inc(jb);end% calculate hy field %for i=1:KE-1for j=1:KEcurl_e=ez(i+1,j)-ez(i,j); ihy(i,j)=ihy(i,j)+fj1(j)*curl_e;hy(i,j)=fi3(i)*hy(i,j)+fi2(i)*0.5* (curl_e+ihy(i,j);endend% Incident Hy vzlues %for j=ja:jbhy(ia-1,j)= hy(ia-1,j)-0.5*ez_inc(j); hy(ib,j)= hy(ib,j)+0.5*ez_inc(j);end% diaplay in the picture %m=moviein(500);x=1:KE;y=1:KE;X,Y=meshgrid(x,y);surf(X,Y,ez);axis(0 60 0 60 -1 1 );set(gcf,Color, white, Number, off, Name, sprintf(Simulation FDTD 2D, Iteration = %i, n);title(sprintf(t = %d,n),Color,1 0 0,FontSize,14); % sprintf(t = %.3f nsec,n*dt*1e9) %xlabel(x);ylabel(y);zlabel(Ez);m(a)=getframe(gcf);a=a+1;%pause(0.005)end