fdtd_1d_ill.m

fdtd一维到三维算法

%***********************************************************************
%     1-D FDTD code with initial condition E(z,0) = E0
%   Illustrates the staggered grid i FDTD
%***********************************************************************
%
%     Program author: Mats Gustafsson
%                     Department of Electroscience
%                     Lund Institute of Technology
%                     Sweden
%                     
%     Date of this version:  October 2002
%
%
%  
%    PEC                                            PEC
%     ----------------------------------------------->
%    z=0                                            z=Lz
%    r=1                                            r=Nz+1
%
%***********************************************************************
% Physical constants
mu0  = 4e-7 * pi;             % Permeability of vacuum
c0   = 299792458;             % Speed of light in vacuum
eps0 = 1/c0^2/mu0;            % Permittivity of vacuum
eta0 = sqrt(mu0/eps0);        % Wave impedance in vacuum  
GHz = 10^9;
mm = 10^(-3);

%***********************************************************************
% Parameter initiation
%***********************************************************************
Lz = 1;                     % Length in meters
Tmax = Lz/c0;                 % Time interval [0,Tmax]
Dz = 20*mm;                    % spatial grid size

R = 1;                  % stabillity number (grid cells per time step)
freq = 1*GHz;                  % center frequency of source excitation

n_pict = 2;                    % plot each n_pict step 
Nz = round(Lz/Dz);              % Number of cells in the z-direction
Dt = Dz*R/c0;                  % Time step
Nt = round(Tmax/Dt);           % Number of time steps
z = linspace(0,Lz,Nz+1);         % z-coordinates
t = Dt*[0:Nt-1]';               % time
lambda = c0/freq;               % center wavelength of source excitation
ppw = lambda/Dz                 % sample points per wave length at the center frequency
omega = 3.0*pi*freq;            % angular frequency

% colormap(ho)
hott=hot(32);
hottt=hott(size(hott,1):-1:1,:);
ho=[1-hottt; hottt];
%***********************************************************************
% Allocate field vectors
%***********************************************************************
Ex = zeros(Nz+1  ,1);          % Electric field
Hy = zeros(Nz,1);          % Magnetic field

res = zeros(2*Nz+1,2*Nt);
%***********************************************************************
%     Wave excitation, initial condition
%***********************************************************************

E0 = sin(omega*(z-Lz/2)/c0).*exp(-omega^2*((z-Lz/2).^2/c0^2));
Ex = E0';
%***********************************************************************
% Time stepping
%***********************************************************************
for n = 1:Nt;
  res1 = zeros(2*Nz+1,2*Nt);
  
  % Update H from n*Dt-Dt/2 to n*Dt+Dt/2
  Hy = Hy - (Dt/mu0/Dz) * diff(Ex,1);                     % main grid
  
  res(2:2:2*Nz,2*n-1) = eta0*Hy;
  res1(2:2:2*Nz,2*n-1) = eta0*Hy;
  % Update E from n*Dt to (n+1)*Dt
  Ex(2:Nz) = Ex(2:Nz) - (Dt/eps0/Dz) * diff(Hy,1);    % main grid except boundary
  
  res(1:2:2*Nz+1,2*n) = Ex;
  res1(1:2:2*Nz+1,2*n) = Ex;
  % Sample the electric field at chosen points
  if mod(n,n_pict) == 0
      figure(2);
      subplot(2,1,1);
      ma = max(max(res));
      imagesc(res',[-ma ma])
      colormap(ho)
      colorbar
      axis xy
      title(['time is ',num2str(n*Dt*10^9),'ns'])
      axis tight;
      xlabel('2*space/Dz');
      ylabel('2*time/Dt');
      
      subplot(2,1,2);
      ma = max(max(res1));
      imagesc(res1',[-ma ma])
      colormap(ho)
      colorbar
      axis xy
      title(['time is ',num2str(n*Dt*10^9),'ns'])
      axis tight;
      xlabel('2*space/Dz');
      ylabel('2*time/Dt');
      
      figure(1)
      subplot(2,1,1);
      plot(z,Ex(1:Nz+1),'LineWidth',1);
      title(['time is ',num2str(n*Dt*10^9),'ns'])
      axis tight;
      ylabel('E-field');
      xlabel('z/m');
      
      subplot(2,1,2);
      plot(z(1:Nz),Hy(1:Nz),'LineWidth',1);
      title(['time is ',num2str(n*Dt*10^9),'ns'])
      axis tight;
      ylabel('H-field');
      xlabel('z/m');
      pause(0.05);
  end
end