RCS-FDTD方柱-球的RCS程序

%FDTD_2D_RCSclear;clc;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*******************************初始化***********************************%%tic;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%v=3*10^8;%波速f=15*10^(9);%频率lamda=v/f;%波长k=2*pi/lamda;%波数epsz=1/(4*pi*9*10^9);%真空介电常数mu=4*pi*10.^(-7);%真空磁导率Z=sqrt(mu/epsz);%真空波阻抗epsilon=1;%相对介电常数sigma=0.0;%电导率%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%N=120;%网格数量L=1000;%迭代次数ddx=lamda/20;%网格尺寸dt=ddx/(2*v);%时间间隔lstart=1;%空间起始行坐标lend=N;%空间终止行坐标rstart=1;%空间起始列坐标rend=N;%空间终止列坐标ia=N/4;%总场区域x左ib=3*N/4;%总场区域x右ja=ia;%总场区域x下jb=ib;%总场区域x上pa=ia-5;%外推区域x左pb=ib+5;%外推区域x右qa=pa;%外推区域x下qb=pb;%外推区域x上length=pb-pa+1;%每个边上的总长度npml=N/8;%PML点数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%方柱参数side=2*lamda;%方柱边长halfgrid=round(side/(2*ddx));%方柱占据网格大小%%媒质参数fori=lstart:lend;%控制媒质分部区域forj=rstart:rend;ga(i,j)=1/(epsilon+sigma*dt/epsz);%求和参量gb(i,j)=sigma*dt/epsz;%求和参量end;end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%源参数spread=8;%脉冲宽度t0=25;%脉冲高度is=N/2;%源的X位置js=N/2;%源的Y位置%%输入平面波ez_inc=zeros(1,N);hx_inc=zeros(1,N);%%平面波吸收条件变量初始化ez_v1=0;ez_v2=0;ez_v3=0;ez_v4=0;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%迭待电磁场参量dz=zeros(N,N);%z方向电荷密度ez=zeros(N,N);%z方向电场iz=zeros(N,N);%z方向电场求和参量hx=zeros(N,N);%x方向磁场hy=zeros(N,N);%y方向磁场ihx=zeros(N,N);%x方向磁场参量ihy=zeros(N,N);%y方向磁场参量%%相位和幅度提取(傅里叶变换)ine=0;inh=0;ez_out=zeros(4,length+1);hx_out=zeros(4,length);hy_out=zeros(4,length);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*********************************PML设置********************************%%%%PML初始化fori=1:N;gi2(i)=1;gi3(i)=1;fi1(i)=0;fi2(i)=1;fi3(i)=1;end;forj=1:N;gj2(j)=1;gj3(j)=1;fj1(j)=0;fj2(j)=1;fj3(j)=1;end;%%阻抗渐变设置fori=1:npml+1;xnum=npml-i+1;xxn=xnum/npml;xn=0.33*(xxn^3);gi2(i)=1/(1+xn);gi2(N-i+1)=1/(1+xn);gi3(i)=(1-xn)/(1+xn);gi3(N-i+1)=(1-xn)/(1+xn);xxn=(xnum-0.5)/npml;xn=0.25*(xxn^3);fi1(i)=xn;fi1(N-i)=xn;fi2(i)=1/(1+xn);fi2(N-i)=1/(1+xn);fi3(i)=(1-xn)/(1+xn);fi3(N-i)=(1-xn)/(1+xn);end;forj=1:npml+1;xnum=npml-j+1;xxn=xnum/npml;xn=0.33*(xxn^3);gj2(j)=1/(1+xn);gj2(N-j+1)=1/(1+xn);gj3(j)=(1-xn)/(1+xn);gj3(N-j+1)=(1-xn)/(1+xn);xxn=(xnum-0.5)/npml;xn=0.25*(xxn^3);fj1(j)=xn;fj1(N-j)=xn;fj2(j)=1/(1+xn);fj2(N-j)=1/(1+xn);fj3(j)=(1-xn)/(1+xn);fj3(N-j)=(1-xn)/(1+xn);end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%***************************迭代求解电场和磁场****************************%%forT=1:L;disp(T);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*************************电场由磁场更新*********************************%%%%TM波Y方向传播forj=2:N-1;ez_inc(j)=ez_inc(j)+0.5*(hx_inc(j-1)-hx_inc(j));end;%%平面波吸收条件ez_inc(1)=ez_v2;ez_v2=ez_v1;ez_v1=ez_inc(2);ez_inc(N)=ez_v3;ez_v3=ez_v4;ez_v4=ez_inc(N-1);%%入射电场傅里叶变换(必须用边界加入)%注意：当计算频谱用的是其他的位置的入射波时候，会产生很大的误差；ine=ine+exp(-sqrt(-1)*2*pi*f*T*dt)*ez_inc(ja-1);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%电荷密度Dzfori=2:N-1;%为了使每个电场周围都有磁场进行数组下标处理forj=2:N;dz(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));end;end;%%脉冲或正弦源的加入source=exp((-0.5)*((t0-T)/spread).^2);%脉冲源ez_inc(5)=source;%TM平面波源%%电荷密度Dz加入射场fori=ia:ib;dz(i,ja)=dz(i,ja)+0.5*hx_inc(ja-1);dz(i,jb)=dz(i,jb)-0.5*hx_inc(jb);end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%电场Ezfori=1:N;forj=1:N;ez(i,j)=ga(i,j)*(dz(i,j)-iz(i,j));iz(i,j)=iz(i,j)+gb(i,j)*ez(i,j);end;end;%%边界电场置零,作为辅助PMLforj=1:N;ez(1,j)=0;ez(N,j)=0;end;fori=1:N;ez(i,1)=0;ez(i,N)=0;end;%%金属边界条件fori=N/2-halfgrid:N/2+halfgrid-1;forj=N/2-halfgrid:N/2+halfgrid-1;ez(i,j)=0;end;end;%%傅里叶变换Ez;s=0;fori=pa:pb+1;s=s+1;ez_out(1,s)=ez_out(1,s)+exp(-sqrt(-1)*2*pi*f*T*dt)*ez(i,qa);%下ez_out(2,s)=ez_out(2,s)+exp(-sqrt(-1)*2*pi*f*T*dt)*ez(i,qb);%上end;s=0;forj=qa:qb+1;s=s+1;ez_out(3,s)=ez_out(3,s)+exp(-sqrt(-1)*2*pi*f*T*dt)*ez(pa,j);%左ez_out(4,s)=ez_out(4,s)+exp(-sqrt(-1)*2*pi*f*T*dt)*ez(pb,j);%右end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*************************磁场由电场更新*********************************%%%%更新平面波磁场forj=1:N-1;hx_inc(j)=hx_inc(j)+0.5*(ez_inc(j)-ez_inc(j+1));end;%%入射磁场傅里叶变换inh=inh+exp(-sqrt(-1)*2*pi*f*(T+1/2)*dt)*hx_inc(ja-1);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%磁场Hxfori=1:N;forj=1:N-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));end;end;%%磁场Hx加入射场fori=ia:ib;hx(i,ja-1)=hx(i,ja-1)+0.5*ez_inc(ja);hx(i,jb)=hx(i,jb)-0.5*ez_inc(jb);end;%%傅里叶变换Hx;s=0;%注意：电场和磁场相差半个时间步fori=pa:pb;s=s+1;hx_out(1,s)=hx_out(1,s)+exp(-sqrt(-1)*2*pi*f*(T+1/2)*dt)*hx(i,qa);%下hx_out(2,s)=hx_out(2,s)+exp(-sqrt(-1)*2*pi*f*(T+1/2)*dt)*hx(i,qb);%上end;s=0;forj=qa:qb;s=s+1;hx_out