fdtd2.m

简单的fdtd算法的matlab实现

% 简化版，介质（非磁性）不能接触边界
N=6;                  %  运算步数
I=200;                    %  空间大小
J=400;                   %  空间大小
c=3e8;                   %  真空光速
L=632.8e-9;
w=2*pi*c/L;              %  角频率=2.9788e+015
e0=(1e-9)/(36*pi);       %  真空的介电常数=8.8419e-012
e=ones(I,J)*e0;
u0=pi*4e-7;                               %  真空的磁导率=1.2566e-006
ds=4e-9;                                 %  波长=2*pi*c/w 步长=波长/40
dt=ds/(2*c);                             %  时间步
g0=0;
g=zeros(I,J);                    %  真空的电导率
gm=0;                                    %  真空的导磁率

sourcePosX=50;
sourcePosY=40;
ArrayPosX=10;
ArrayPosY=100;
ArrayLength=100;
ArrayWidth=3;
period=6;
ArrayNumber=floor( (I-2*ArrayPosX)/period );
%  电磁场分量定义
Hz=zeros(I,J);
Ex=zeros(I,J);
Ey=zeros(I,J);
Hz1=zeros(I,J);
Ex1=zeros(I,J);
Ey1=zeros(I,J);

% Source -----------------------------------------------------------
Source=zeros(I,J);
%Source(sourcePosX,sourcePosY)=1;
x1=sourcePosX;
x2=I-sourcePosX*2;
 for x=1:1:x2-x1+1;
     Source(x+x1-1,sourcePosY)=exp(-4*pi*(x-(x2-x1)/2)^2/(x2-x1)^2);
 end
% --------------------------------------------------------------------
% Object --------------------------------------------------------------
nn=L*1e6;   % L is wave length
coon=(-6.21198035452752*nn.^5+ 36.9470450716619*nn.^4-83.1702487377279*nn.^3+...
   94.5079035255191*nn.^2-48.61825472792*nn+9.64235361750187)
eee=(-108.909148353719*nn.^5+ 671.044293885632*nn.^4-1547.0900731642*nn.^3+ ...
   1674.14670244375*nn.^2-773.494216929895*nn+130.284211701206)


g1=(coon+i*eee)*e0.*w                   % 电导率


%   Set Ag
% g(20:80,15:J)=g1;                      % 3.72e8; %0.91*(1e-9)/(36*pi)*w;

%g(1:100,LL:R)=g1; %(-16.4)*(-i)*e0.*w;                      % 3.72e8; %0.91*(1e-9)/(36*pi)*w;
%e(20:80,80:J)=e0; 

%   Set Air

for count=1:ArrayNumber
    g(ArrayPosX+(count-1)*period:ArrayPosX+(count-1)*period+ArrayWidth,ArrayPosY:ArrayLength)=g1;
end
%k=49;
%for ii=1:1   
   % g(k:k+2,20:J)=0;
   % k=k+5;
%end
%e(20:80,80:J)=e0;

%for k=1:30     
 %  g(i:i+2,140:180)=g1;
 % i=i+6;
    %e(i:i+1,15:60)=e0;
%end

% -----------------------------------------------------------------------

CA=(1-g*dt/2./e)./(1+g*dt/2./e);    % CA=1
CB=dt./e./(1+g*dt/2./e);            % CB=0.094248
CP=(1-gm*dt/2/u0)./(1+gm*dt/2/u0);    % CP=1
CQ=dt/u0/(1+gm*dt/2/u0);              % CQ=6.6315e-007

% 将运算移出循环外，加快程序运行。-----------(Used in Mur boundary condition)
Boundary1=(c*dt-ds)/(c*dt+ds);
Boundary2=c^2*e0*dt/(2*(c*dt+ds));
CORNER=(c*dt-sqrt(2)*ds)./(c*dt+sqrt(2)*ds);

%开始计算
disp('正在计算......');

for n=0:N;
% Switch ----------------------------------------------------------------
if n<pi/w/dt
    Switch=0.5*(1-cos(pi*n/(pi/w/dt)));
else 
    Switch=1;
end
% ------------------------------------------------------------------------
% time evolution ----------------------------------------------------------
x=1:I-1;
y=1:J-1;
Ex1(x,y)=CA(x,y).*Ex(x,y)+CB(x,y).*(Hz(x,y+1)-Hz(x,y))/ds;
Ey1(x,y)=CA(x,y).*Ey(x,y)-CB(x,y).*(Hz(x+1,y)-Hz(x,y))/ds;

x=2:I-1;
y=2:J-1;
Hz1(x,y)=CP*Hz(x,y)-CQ*(Ey1(x,y)-Ey1(x-1,y)-Ex1(x,y)+Ex1(x,y-1))/ds;
Hz1=exp(-i*dt*n*w)*Source*Switch/sqrt(u0/e0)+Hz1;
% -------------------------------------------------------------------------
%  计算截断边界------------------------------------------------------------
%左边界
Hz1(1,y)=Hz(2,y)+Boundary1*(Hz1(2,y)-Hz(1,y))+Boundary2*(Ex1(1,y)-Ex1(1,y-1)+Ex1(2,y)-Ex1(2,y-1));
%右边界
Hz1(I,y)=Hz(I-1,y)-Boundary1*(Hz(I,y)-Hz1(I-1,y))+Boundary2*(Ex1(I-1,y)-Ex1(I-1,y-1)+Ex1(I,y)-Ex1(I,y-1));
%下边界
Hz1(x,1)=Hz(x,2)+Boundary1*(Hz1(x,2)-Hz(x,1))-Boundary2*(Ey1(x,1)-Ey1(x-1,1)+Ey1(x,2)-Ey1(x-1,2));
%上边界
Hz1(x,J)=Hz(x,J-1)-Boundary1*(Hz(x,J)-Hz1(x,J-1))-Boundary2*(Ey1(x,J-1)-Ey1(x-1,J-1)+Ey1(x,J)-Ey1(x-1,J));
%左下角
Hz1(1,1)=Hz(2,2)+CORNER*(Hz1(2,2)-Hz(1,1));
%左上角
Hz1(1,J)=Hz(2,J-1)+CORNER*(Hz1(2,J-1)-Hz(1,J));
%右上角
Hz1(I,J)=Hz(I-1,J-1)-CORNER*(Hz(I,J)-Hz1(I-1,J-1));
%右下角
Hz1(I,1)=Hz(I-1,2)-CORNER*(Hz(I,1)-Hz1(I-1,2));
% -------------------------------------------------------------------------
Hz=Hz1;
Ex=Ex1;
Ey=Ey1;

if mod(n,100)==0
fprintf('the %d step is finished\n',n);
pcolor(0:ds:ds*(J-1),0:ds:ds*(I-1),abs(Hz));
%E=abs(Ex(1:I,1:J)).^2+abs(Ey(1:I,1:J)).^2;
%E=abs(Ex).^2+abs(Ey).^2;
%pcolor(E);
axis equal
shading interp
%shading flat
colorbar('horiz')
xlabel('Y')
ylabel('X')
%pause(0.01)
end

end