fdtd.m

基于Matlab的FDTD算法

% MATLAB FDTD Program 
% 
% 2D FDTD Code with Mur's ABCs
% 22 July, 1999.
% 

clear all;
tic
% --------- Define Constants -----------
c= 2.99795637e+8;
eps0 = 8.854e-12;
mu0 = 1.2566306e-6;

% --------- Define some parameters ----------
dx = 1e-3;      % Space increment
dy = 1e-3;      % Space increment

dt = 1/((sqrt((1/dx^2)+(1/dy^2))*c));     % Courants stablity criterion
a = 0;
nx = 221 + a;       % number of cells in X direction
ny = 101;           % number of cells in Y direction

nmax = 500;     % total number of time steps

w_pulse = 15e-12;   % Width of source gaussian pulse.
t0 = 4.0 * w_pulse; % pulse truncation value

% ------- Initialise field vectors and update constants --------
ey = zeros(nx+1,ny+1);
hz = zeros(nx+1,ny+1);
ex = zeros(nx+1,ny+1);

cbx = ones(nx+1,ny+1);
dbx = ones(nx+1,ny+1);
cby = ones(nx+1,ny+1);
dby = ones(nx+1,ny+1);

erx = ones(nx+1,ny+1);
ery = ones(nx+1,ny+1);
%---
erx(1+a:151+a,1:41)=inf;
erx(151+a:221+a,1:16)=inf;

erx(1+a:151+a,60:101)=inf; %altered
erx(151+a:221+a,75:101)=inf;

erx(166+a:221+a,30:61)=inf;
%---
ery(1+a:150+a,1:41)=inf; %altered
ery(151+a:221+a,1:16)=inf;

ery(1+a:150+a,61:101)=inf;%altered
ery(151+a:221+a,76:101)=inf;

ery(166+a:221+a,31:61)=inf;
%---
%---
%exc_rec = zeros(nmax,1);
%time_rec = zeros(nmax,1);

exbufl1 = zeros(nx+1,2);
exbufl2 = zeros(nx+1,2);
exbufr1 = zeros(nx+1,2);
exbufr2 = zeros(nx+1,2);

eybufl1 = zeros(2,ny+1);
eybufl2 = zeros(2,ny+1);
eybufr1 = zeros(2,ny+1);
eybufr2 = zeros(2,ny+1);



%for i=1:nx+1
%   for j=1:ny+1
cbx(1:nx+1,1:ny+1) = (dt ./ (eps0.*erx(1:nx+1,1:ny+1)*dx));
dbx(1:nx+1,1:ny+1) = (dt / (mu0*dx));
%   end
%end

%for i=1:nx+1
%   for j=1:ny+1
cby(1:nx+1,1:ny+1) = (dt ./ (eps0.*ery(1:nx+1,1:ny+1)*dy));
dby(1:nx+1,1:ny+1) = (dt / (mu0 * dy));
%   end
%end



kx =  (c*dt - dx) / (c*dt + dx);  % MUR's ABC constant
ky =  (c*dt - dy) / (c*dt + dy);  % MUR's ABC constant


%=================================
% ----- Main FDTD Loop --------
%=================================
x=1;
for n=1:nmax,

   % Calculate Ey
%   for i=2:nx
%      for j=1:ny
ey(2:nx,1:ny) = ey(2:nx,1:ny) - cbx(2:nx,1:ny).*(hz(2:nx,1:ny) - hz(1:nx-1,1:ny));
%     end
%  end
  
  % Calculate Ex
%  for i=1:nx
%     for j=2:ny
ex(1:nx,2:ny) = ex(1:nx,2:ny) + cby(1:nx,2:ny).*(hz(1:nx,2:ny) - hz(1:nx,1:ny-1));
%     end
%  end


  % Apply excitation
  t = n * dt;
  exc = exp((-1 * ((t - t0)^2)) / (w_pulse * w_pulse));
  %if(n<=8)
  %   exc = sin(2*pi*21.2e9*n*dt);
  %else
  %   exc=1.052*exp(-0.116389*n);
%	end
  
  %exc = 0.01*sin(2*pi*1e9*n*dt);
  
  %exc_rec0(n) = exc;            % store a record of the excitation
  
  hz(20+a,41:60) = hz(20+a,41:60) + 0.1*exc;
  
  % Apply MUR ABC ------ ex field ---------------------
	exbufl1(:,2) = exbufl1(:,1);
   exbufl1(:,1) = ex(:,2);

   exbufr1(:,2) = exbufr1(:,1);
   exbufr1(:,1) = ex(:,ny);

   ex(:,1)    = exbufl1(:,2) + kx*(exbufl1(:,1)-exbufl1(:,2));
   ex(:,ny+1) = exbufr1(:,2) + kx*(exbufr1(:,1)-exbufr1(:,2));
   
   % Apply MUR ABC ------ ey field ---------------------
  	eybufl1(2,:) = eybufl1(1,:);
   eybufl1(1,:) = ey(2,:);

   eybufr1(2,:) = eybufr1(1,:);
   eybufr1(1,:) = ey(nx,:);

   ey(1,:)    = eybufl1(2,:) + ky*(eybufl1(1,:)-eybufl1(2,:));
   ey(nx+1,:) = eybufr1(2,:) + ky*(eybufr1(1,:)-eybufr1(2,:));
   
%  time_rec(n) = ey(100);
%for i=1:nx+1
%   for j=1:ny+1

if(mod(n,10)==0)
   e_rec(:,:,x) = sqrt(ey(:,:).^2+ex(:,:).^2);
   %ey_rec(:,:,x)=ey(:,:);
   %ex_rec(:,:,x)=ex(:,:);
   
   %e_s11(x) = sqrt(ey(15+a,50).^2+ex(15+a,50).^2);
   %e_s21(x) = sqrt(ey(207+a,67).^2+ex(207+a,67).^2);
   %e_s31(x) = sqrt(ey(207+a,22).^2+ex(207+a,22).^2);
      
   x=x+1;
end
%e_rec(n)= sqrt(ey(30,50).^2+ex(30,50).^2);
%   end
%end

  % Calculate Hz
%  for i=1:nx
%     for j=1:ny
hz(1:nx,1:ny) = hz(1:nx,1:ny) - dbx(1:nx,1:ny).*(ey(2:nx+1,1:ny) - ey(1:nx,1:ny)) + dby(1:nx,1:ny).*(ex(1:nx,2:ny+1) - ex(1:nx,1:ny));
%    end
% end

% --- Define the Metal Area ---

%hz(70,60)=0;
%hz(70,40)=0;

%ex(70:90,40)=0;

end  % Main time stepping loop.

toc
% Post Processing.
%mt2    % animated wave propagation
w=size(e_rec);
n=w(3);

M=moviein(n);

for j=1:n
     pcolor(e_rec(:,:,j));
     shading interp;
     M(:,j)=getframe;
  end
movie(M)
%ssp
%time_rec_refl = time_rec;

%pulse_r = time_rec_in-time_rec_refl;
%pulse_i = time_rec_in;
%fft_i = fft(pulse_i);
%fft_r = fft(pulse_r);

%df = 1/dt;

%refl = abs(fft_r/fft_i); 


%for i=1:100,
%  freq(i) = i*df;
%end

%plot(freq(1:100),20*log10(refl(1:100)))