fdtd_1d.m

一维FDTD 一维FDTD 一维FDTD 一维FDTD

% FDTD 1D problem

%parameters of the incident wave
f=3*10^9;
c=3.0*10^8;
%parameters of the mesh 
namta=c/f;
omiga=2*pi*f;
dx=namta/25;

dt=dx/c/2;
%parameters of the material 
mu=4*pi*10^(-7);
epsr=11.9;
eps=epsr*1/(c*c*mu);
% stimulus parameters
x=0.5;
NX=round(x/dx);
NT=200;



%initialization
ex=zeros(1,NX+1);
hy=zeros(1,NX);

xx=linspace(dx,x,NX);
subplot(2,1,1),plot(x,ex(1:NX),'r');
axis([0 x -1 1]);
xlabel('x'),ylabel('ex');
subplot(2,1,2),plot(x,hy(1:NX),'b');
axis([0 x -0.01 0.01]);
xlabel('x'),ylabel('hy');

rect=get(gcf,'Position');
rect(1:2)=[0 0];

M=moviein(NT/2,gcf,rect);

% main loop
for n=1:NT
    
    %if n>=1 & n<=10
        ex(1)=ex(1)+exp(-4*pi*(n*dt)^2);   %source
        %end
    % simple boundary condition
    temp=ex(NX); 
      for k=2:NX
        ex(k)=ex(k)-dt/(eps*dx)*(hy(k)-hy(k-1));
    end
    ex(NX+1)=temp;
   
    for k=1:NX
        hy(k)=hy(k)-dt/(mu*dx)*(ex(k+1)-ex(k));
    end


if mod(n,2)==0;

rtime=num2str(n);

subplot(2,1,1),plot(xx,ex(1:NX),'r');
axis([0 x -1 1]);
title(['time step = ',rtime]);
xlabel('x ');ylabel('EX');

subplot(2,1,2),plot(xx,hy,'b');
axis([0 x -0.01 0.01]);
title(['time step= ',rtime]);
xlabel('x ');ylabel('HY');

M(:,n/2)=getframe(gcf,rect);

end

end

movie(gcf,M,0,10,rect);