FDTD仿真二维TE波方金属柱和圆金属柱的双站RCS

% FDTD仿真二维TE波目标双站RCS；

% 编写时间：2012年4月28日；

% 编写者：tingo;

% Email:guoqianghe22@https://www.360docs.net/doc/552731276.html,.

%%

clear;

clf;

clc;

close all;

%% 计算中用到的一些电磁常量；

mu0=4*pi*(1e-7); %磁导率；

epsilon0=(1/36/pi)*(1e-9); %介电常数；

Vvacuum=1/sqrt(mu0*epsilon0); % 真空中的波传播速度；Zc=sqrt(mu0/epsilon0);

%% 求解区域及网格划分

namda=1e-2; % 1GHz对应的波长；

dx=namda/40;

dy=dx;

nx=150;

ny=nx;

Lx=nx*dx;

Ly=Lx;

a=2*namda;

dt=dx/Vvacuum/2; % 时间步长；

itmax=1000; % 仿真时间步；

%% PML层电磁参数设置；

Rmax=1e-3; % 电磁波垂直入射的反射系数；iebc=4; % 左右两边PML层网格数目；

jebc=4; % 上下两边PML层网格数目；orderbc=4; % pML层中电导率分布阶数；iebcf=nx+2*iebc;

iebcb=nx+2*iebc;

%% 为储存电磁数据的矩阵变量分配内存；%========FDTD 区域=========%

Ex=zeros(nx,ny+1);

Ey=zeros(nx+1,ny);

Hz=zeros(nx,ny);

caex=zeros(nx,ny+1);

cbex=zeros(nx,ny+1);

caey=zeros(nx+1,ny);

cbey=zeros(nx+1,ny);

cahz=zeros(nx,ny);

cbhz=zeros(nx,ny);

%=========PML层==========%

% front;

Exbcf=zeros(iebcf,jebc);

Eybcf=zeros(iebcf+1,jebc);

Hzxbcf=zeros(iebcf,jebc);

Hzybcf=zeros(iebcf,jebc);

caexbcf=zeros(iebcf,jebc); %Exbcf计算系数1, cbexbcf=zeros(iebcf,jebc); %Exbcf计算系数2；caeybcf=zeros(iebcf+1,jebc);

cbeybcf=zeros(iebcf+1,jebc);

cahzxbcf=zeros(iebcf,jebc);

cbhzybcf=zeros(iebcf,jebc);

% back;

Exbcb=zeros(iebcb,jebc); % jebc+1??????????????????????????? Eybcb=zeros(iebcb+1,jebc);

Hzxbcb=zeros(iebcb,jebc);

Hzybcb=zeros(iebcb,jebc);

caexbcb=zeros(iebcb,jebc);

cbexbcb=zeros(iebcb,jebc);

caeybcb=zeros(iebcb+1,jebc);

cbeybcb=zeros(iebcb+1,jebc);

cahzxbcb=zeros(iebcb,jebc);

cbhzxbcb=zeros(iebcb,jebc);

cahzybcb=zeros(iebcb,jebc);

cbhzybcb=zeros(iebcb,jebc);

% left

Exbcl=zeros(iebc,ny+1);

Eybcl=zeros(iebc,ny);

Hzxbcl=zeros(iebc,ny);

Hzybcl=zeros(iebc,ny);

caexbcl=zeros(iebc,ny+1);

caeybcl=zeros(iebc,ny); cbeybcl=zeros(iebc,ny); cahzxbcl=zeros(iebc,ny); cbhzxbcl=zeros(iebc,ny); cahzybcl=zeros(iebc,ny); cbhzybcl=zeros(iebc,ny); % right;

Exbcr=zeros(iebc,ny+1); Eybcr=zeros(iebc,ny); Hzxbcr=zeros(iebc,ny); Hzybcr=zeros(iebc,ny);

caexbcr=zeros(iebc,ny+1); cbexbcr=zeros(iebc,ny+1); caeybcr=zeros(iebc,ny); cbeybcr=zeros(iebc,ny); cahzxbcr=zeros(iebc,ny); cbhzxbcr=zeros(iebc,ny); cahzybcr=zeros(iebc,ny);

%% 仿真区域媒质参数设定；nmedium=2; % 媒质数目；epsilonr=[1,1]; % 相对介电常数；mur=[1,1]; % 相对磁导率

sigma=[0,6e7]; % 电导率；

rou=[0,0]; % 磁阻率；

ca=zeros(1,nmedium);

cb=zeros(1,nmedium);

cp=zeros(1,nmedium);

cq=zeros(1,nmedium);

for i=1:nmedium

epsilon=epsilon0*epsilonr(i);

mu=mu0*mur(i);

eaf=sigma(i)*dt/2/epsilon;

haf=rou(i)*dt/2/mu;

ca(i)=(1-eaf)/(1+eaf);

cb(i)=dt/epsilon/(1+eaf);

cp(i)=(1-haf)/(1+haf);

cq(i)=dt/mu/(1+haf);

end

%% 自由空间FDTD迭代计算系数；

% Ex;

caex(:,:)=ca(1);

cbex(:,:)=cb(1);

% Ey;

caey(:,:)=ca(1);

cbey(:,:)=cb(1);

% Hz;

cahz(:,:)=cp(1);

cbhz(:,:)=cq(1);

%% 散射体模型；

ObjType=input('请输入散射体模型（1. 方金属柱；2. 圆金属柱）：'); switch ObjType

case 1

% 金属方柱；

nx1=round( (Lx-a)/2/dx ); % 金属方柱x方向开始网格编码；nx2=round( (Lx+a)/2/dx ); % 金属方柱x方向截止处网格编码；ny1=nx1; % 金属方柱y方向开始网格编码；

ny2=nx2; % 金属方柱y方向截止处网格编码；

% 金属方柱区域FDTD计算系数；

caex( (nx1+1):(nx2 ) , ( ny1+1):(ny2+1) )=ca(2);

cbex( (nx1+1):(nx2 ) , ( ny1+1):(ny2+1) )=cb(2);

caey( (nx1+1):(nx2+1) , (ny1+1):(ny2 ) )=ca(2);

cbey( (nx1+1):(nx2+1) , (ny1+1):(ny2 ) )=cb(2);

cphz( nx1+1 : nx2 , (ny1+1) : ny2 )=cp(2);

cqhz( nx1+1 : nx2 , (ny1+1) : ny2 )=cq(2);

case 2

% 圆金属柱；

r=namda; % 圆介质柱半径；

xcenter=Lx/2; % 圆柱的中心x坐标；

ycenter=Ly/2; % 圆柱的中心y坐标；

% 圆柱区域FDTD计算系数；

% Ex;

for i=1:nx

for j=2:ny

x=(i-0.5)*dx;

y=(j-1)*dy;

dist=sqrt( (x-xcenter)^2+(y-ycenter)^2 ); if(dist<=r)

caex(i,j)=ca(2);

cbex(i,j)=cb(2);

end

end

end

% Ey;

for i=2:nx

for j=1:ny

x=(i-1)*dx;

y=(j-0.5)*dy;

dist=sqrt( (x-xcenter)^2+(y-xcenter)^2 ); if(dist<=r)

caey(i,j)=ca(2);

cbey(i,j)=cb(2);

end

end

end

otherwise

disp('没有电磁目标，请重新运行程序，正确输入！');

end

%% PML层区域煤质参数；

%==============================================%

% Front

%==============================================%

delbc=jebc*dy;

sigmabcm=-log(Rmax/100.0)*(orderbc+1)*epsilon0*Vvacuum/2/delbc; bcfactor=sigmabcm/dy/(delbc^orderbc)/(orderbc+1);

caexbcf(1:iebcf,1)=1;

cbexbcf(1:iebcf,1)=0;

for j=2:jebc

y1=(jebc-j+1.5)*dy;

y2=(jebc-j+0.5)*dy;

sigmay=bcfactor*(y1^(orderbc+1)-y2^(orderbc+1) );

ca1=exp(-sigmay*dt/epsilon0);

cb1=(1-ca1)/sigmay/dy ;

caexbcf(1:iebcf,j)=ca1;

cbexbcf(1:iebcf,j)=cb1;

end

sigmay=bcfactor*(0.5*dy)^(orderbc+1);

ca1=exp(-sigmay*dt/epsilon0);

cb1=(1-ca1)/sigmay/dy ;

caex(1:nx,1)=ca1;

cbex(1:nx,1)=cb1;

%左、右边PML区域与前PML区域的交界面处，j=1; % exbcl 和exbcr 在j=1;

caexbcl(1:iebc,1)=ca1;

cbexbcl(1:iebc,1)=cb1;

caexbcr(1:iebc,1)=ca1;

cbexbcr(1:iebc,1)=cb1;

for j=1:jebc;

y1=(jebc-j+1)*dy;

y2=(jebc-j)*dy;

sigmay=bcfactor*(y1^(orderbc+1)-y2^(orderbc+1) );

rouy=sigmay*mu0/epsilon0;

cp1=exp(-rouy*dt/mu0);

cq1=(1-cp1)/rouy/dy;

cahzybcf(1:iebcf,j)=cp1;

cbhzybcf(1:iebcf,j)=cq1;

end

% PML层中Ey和Hzx迭代系数计算：sigmax=0;,roux=0;

for j=1:jebc

ca1=1;

cb1=dt/dx/epsilon0;

cp1=1;

cq1=dt/dx/mu0;

caeybcf(iebc+2:iebc+nx,j)=ca1; % 修正此处“iebc+1:iebc+nx+1”为"iebc+2:iebc+nx"，2012/4/27.

cbeybcf(iebc+2:iebc+nx,j)=cb1; % 修正此处“iebc+1:iebc+nx+1”为"iebc+2:iebc+nx"，2012/4/27.

cahzxbcf(iebc+1:iebc+nx,j)=cp1;

cbhzxbcf(iebc+1:iebc+nx,j)=cq1;

end

%==============================================% % Back

%==============================================%

caexbcb(1:iebcf,jebc)=1;

cbexbcb(1:iebcf,jebc)=0;

for j=1:jebc-1

y1=(j+0.5)*dy;

y2=(j-0.5)*dy;

sigmay=bcfactor*(y1^(orderbc+1)-y2^(orderbc+1) );

ca1=exp(-sigmay*dt/epsilon);

cb1=(1-ca1)/sigmay/dy ;

caexbcb(1:iebcf,j)=ca1;

cbexbcb(1:iebcf,j)=cb1;

end

sigmay=bcfactor*(0.5*dy)^(orderbc+1);

ca1=exp(-sigmay*dt/epsilon);

cb1=(1-ca1)/sigmay/dy;

caex(1:nx,ny+1)=ca1;

cbex(1:nx,ny+1)=cb1;

% 左、右边PML区域与后PML区域的交界面处，j=ny+1; % exbcl 和exbcr 在j=1;

caexbcl(1:iebc,ny+1)=ca1;

cbexbcl(1:iebc,ny+1)=cb1;

caexbcr(1:iebc,ny+1)=ca1;

cbexbcr(1:iebc,ny+1)=cb1;

for j=1:jebc;

y1=j*dy;

y2=(j-1)*dy;

sigmay=bcfactor*(y1^(orderbc+1)-y2^(orderbc+1) );

rouy=sigmay*mu0/epsilon0;

cp1=exp(-rouy*dt/mu0);

cq1=(1-cp1)/rouy/dy;

cahzybcb(1:iebcf,j)=cp1;

cbhzybcb(1:iebcf,j)=cq1;

end

% PML层中Ey和Hzx迭代系数计算：sigmax=0;,roux=0;

for j=1:jebc

ca1=1;

cb1=dt/dx/epsilon0;

cp1=1;

cq1=dt/dx/mu0;

caeybcb(iebc+2:iebc+nx,j)=ca1; % 修正此处“iebc+1:iebc+nx+1”为"iebc+2:iebc+nx"，2012/4/27.

cbeybcb(iebc+2:iebc+nx,j)=cb1; % 修正此处“iebc+1:iebc+nx+1”为"iebc+2:iebc+nx"，2012/4/27.

cahzxbcb(iebc+1:iebc+nx,j)=cp1;

cbhzxbcb(iebc+1:iebc+nx,j)=cq1;

end

%==============================================%

% Left

%==============================================%

delbc=iebc*dx;

sigmabcm=-log(Rmax/100)*(orderbc+1)*epsilon0*Vvacuum/2/delbc; bcfactor=sigmabcm/dx/(delbc^orderbc)/(orderbc+1);

caeybcl(1,1:ny)=1;

cbeybcl(1,1:ny)=0;

caeybcf(1,1:jebc)=1;

cbeybcf(1,1:jebc)=0;

caeybcb(1,1:jebc)=1;

cbeybcb(1,1:jebc)=0;

for i=2:iebc

x1=(iebc-i+1.5)*dx;

x2=(iebc-i+0.5)*dx;

sigmax=bcfactor*(x1^(orderbc+1)-x2^(orderbc+1));

ca1=exp(-sigmax*dt/epsilon0);

cb1=(1-ca1)/(sigmax*dx);

caeybcl(i,1:ny)=ca1;

cbeybcl(i,1:ny)=cb1;

% front 左片区;

caeybcf(i,1:jebc)=ca1;

cbeybcf(i,1:jebc)=cb1;

%back 左片区

caeybcb(i,1:jebc)=ca1;

cbeybcb(i,1:jebc)=cb1;

end

sigmax= bcfactor*(0.5*dx)^(orderbc+1);

ca1=exp(-sigmax*dt/epsilon0);

cb1=(1-ca1)/(sigmax*dx);

caey(1,1:ny)=ca1;

cbey(1,1:ny)=cb1;

caeybcf(iebc+1,1:jebc)=ca1;

cbeybcf(iebc+1,1:jebc)=cb1;

caeybcb(iebc+1,1:jebc)=ca1;

cbeybcb(iebc+1,1:jebc)=cb1;

for i=1:iebc

x1=(iebc-i+1)*dx;

x2=(iebc-i)*dx;

sigmax=bcfactor*(x1^(orderbc+1)-x2^(orderbc+1)); roux=sigmax*mu0/epsilon0;

cp1=exp(-roux*dt/mu0);

cq1=(1-cp1)/roux/dx;

cahzxbcl(i,1:ny)=cp1;

cbhzxbcl(i,1:ny)=cq1;

% front 左片区；

cahzxbcf(i,1:jebc)=cp1;

cbhzxbcf(i,1:jebc)=cq1;

% back 左片区

cahzxbcb(i,1:jebc)=cp1;

cbhzxbcb(i,1:jebc)=cq1;

end

% PML层中Ex和Hzy迭代系数计算：sigmay=0;,rouy=0; for i=1:iebc

ca1=1;

cb1=dt/dy/epsilon0;

cp1=1;

cq1=dt/dy/mu0;

caexbcl(i,2:ny)=ca1;

cbexbcl(i,2:ny)=cb1;

cahzybcl(i,1:ny)=cp1;

cbhzybcl(i,1:ny)=cq1;

end

%==============================================% % Right

%==============================================% caeybcr(iebc,1:ny)=1;

cbeybcr(iebc,1:ny)=0;

caeybcf(iebcf+1,1:jebc)=1;

cbeybcf(iebcf+1,1:jebc)=0;

caeybcb(iebcb+1,1:jebc)=1;

cbeybcb(iebcb+1,1:jebc)=0;

for i=1:iebc-1

x1=(i+0.5)*dx;

x2=(i-0.5)*dx;

sigmax=bcfactor*(x1^(orderbc+1)-x2^(orderbc+1));

ca1=exp(-sigmax*dt/epsilon0);

cb1=(1-ca1)/(sigmax*dx);

caeybcr(i,1:ny)=ca1;

cbeybcr(i,1:ny)=cb1;

% front 右片区;

caeybcf(iebc+nx+1+i,1:jebc)=ca1;

cbeybcf(iebc+nx+1+i,1:jebc)=cb1;

%back 右片区

caeybcb(iebc+nx+1+i,1:jebc)=ca1;

cbeybcb(iebc+nx+1+i,1:jebc)=cb1;

end

% Ey在FDTD右边界与PML分界线上的计算系数；sigmax= bcfactor*(0.5*dx)^(orderbc+1);

ca1=exp(-sigmax*dt/epsilon0);

cb1=(1-ca1)/(sigmax*dx);

caey(nx+1,1:ny)=ca1;

cbey(nx+1,1:ny)=cb1;

caeybcf(iebc+nx+1,1:jebc)=ca1;

cbeybcf(iebc+nx+1,1:jebc)=cb1;

caeybcb(iebc+nx+1,1:jebc)=ca1;

cbeybcb(iebc+nx+1,1:jebc)=cb1;

for i=1:iebc

x1=i*dx;

x2=(i-1)*dx;

sigmax=bcfactor*(x1^(orderbc+1)-x2^(orderbc+1)); roux=sigmax*mu0/epsilon0;