FDTD_work_fdtd_matlab_电磁波仿真_

% FDTD_2D_RCS % 本程序实现FDTD仿真2维TM波金属柱散射 % TM波沿Y方向传播，X方向极化 % 采用UPML吸收，脉冲源激励 clear;clc; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%*******************************初始化***********************************%% tic; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% v=3*10^8; % 波速 f=0.3e9; % 频率 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; % N=input('Please input the number of grid cells N:'); % 迭代次数 L=500; % L=input('Please input the number of iterative steps L:'); ddx=lamda/20; % 网格尺寸 dt=ddx/(2*v); % 时间间隔 lstart=1; % 空间起始行坐标 lend=N; % 空间终止行坐标 rstart=1; % 空间起始列坐标 rend=N; % 空间终止列坐标 %********************************总场区域**********************************% ia=N/5; % 总场区域x左 ib=4*N/5; % 总场区域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; % 每个边上的总长度 % PML区域范围 npml=10; % npml=input('Please input the number of PML cells npml:'); %*******************************方柱参数***********************************% side=1.5*lamda; % 柱边长 halfgrid=round(side/(2*ddx)); % 柱占据网格大小,这里只是一半的尺寸 %********************************媒质参数***********************************% for i=lstart:lend; % 控制媒质分部区域 for j=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方向磁场参量,PML中间变量 ihy=zeros(N,N); % y方向磁场参量,PML中间变量 %% 相位和幅度提取(傅里叶变换) ine=0; inh=0; ez_out=zeros(4,length+1);%length为外推边界的长度 hx_out=zeros(4,length); hy_out=zeros(4,length); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%*********************************PML设置********************************%% %% PML初始化 %% 若fj1 和fi1 为0，而其他参数为1，则为自由空间中的情况 for i=1:N; gi2(i)=1; gi3(i)=1; fi1(i)=0; fi2(i)=1; fi3(i)=1; end; for j=1:N; gj2(j)=1; gj3(j)=1; fj1(j)=0; fj2(j)=1; fj3(j)=1; end; %% 阻抗渐变设置 %% 匹配层介质数据也具有对称性 for i=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; for j=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; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%***************************迭代求解电场和磁场****************************%% %%*******************************总循环开始********************************%% for T=1:L; disp(T); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%*************************电场由磁场更新*********************************%% %% TM波Y方向传播 for j=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(-j*2*pi*f*T*dt)*ez_inc(ja-1);%连接边界处的频率电场值 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% 电场通量Dz for i=2:N-1; % 为了使每个电场周围都有磁场进行数组下标处理，没有计算最边上的四组数据。 for j=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); % 脉冲源 % source=sin(2*pi*f*T*dt); % 正弦波源 ez_inc(1)=source; % TM平面波源 % dz(is,js)=source; % 源 % dz(is,js)=source; %% 电场通量Dz加入射场 %% 连接边界分离 for i=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; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% 电场Ez for i=1:N; for j=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; %% 边界电场置零,作为辅助PML %% 场已经耗散贻尽，几乎不影响正常场分布 for j=1:N; ez(1,j)=0; ez(N,j)=0; end; for i=1:N; ez(i,1)=0; ez(i,N)=0; end; r=round(halfgrid);%圆半径 lcx=N/2 lcy=N/2-halfgrid;%左半圆圆心 rcx=N/2; rcy=N/2+halfgrid-1;%右半圆圆心 %% 金属边界条件 for i=qa:qb; for j=pa:pb; if (sqrt((i-lcx)^2+(j-lcy)^2)<=r)||(sqrt((i-rcx)^2+(j-rcy)^2)<=r)||((i>=(N/2-halfgrid))&&(i<=(N/2+halfgrid-1)))&&((j>=(N/2-halfgrid))&&(j<=(N/2+halfgrid-1))) ez(i,j)=0; else ez(i,j)=ez(i,j)*1; end; end; end; % 傅里叶变换Ez; s=0; for i=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; for j=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; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%