fdtd.m

A SIMPLE BUT USFUL 2D FDTD MATLAB CODE3

% Physical constants
eps0 = 8.8541878e-12;         % Permittivity of vacuum
mu0  = 4e-7 * pi;             % Permeability of vacuum
c0   = 299792458;             % Speed of light in vacuum

% Parameter initiation
Lx = .05; Ly = .04; Lz = .03; % Cavity dimensions in meters
Nx =  25; Ny =  20; Nz =  15; % Number of cells in each direction
Cx = Nx / Lx;                 % Inverse cell dimensions
Cy = Ny / Ly;
Cz = Nz / Lz;
Nt = 1024;                    % Number of time steps
Dt = 1/(c0*norm([Cx Cy Cz])); % Time step

% Allocate field matrices
Ex = zeros(Nx  , Ny+1, Nz+1);
Ey = zeros(Nx+1, Ny  , Nz+1);
Ez = zeros(Nx+1, Ny+1, Nz  );
Hx = zeros(Nx+1, Ny  , Nz  );
Hy = zeros(Nx  , Ny+1, Nz  );
Hz = zeros(Nx  , Ny  , Nz+1);

% Allocate time signals
Et = zeros(Nt,3);

% Initiate fields (near but not on the boundary)
Ex(1,2,2) = 1;
Ey(2,1,2) = 2;
Ez(2,2,1) = 3;

% Time stepping
for n = 1:Nt;
 
  % Update H everywhere
  Hx = Hx + (Dt/mu0)*(diff(Ey,1,3)*Cz - diff(Ez,1,2)*Cy);
  Hy = Hy + (Dt/mu0)*(diff(Ez,1,1)*Cx - diff(Ex,1,3)*Cz);
  Hz = Hz + (Dt/mu0)*(diff(Ex,1,2)*Cy - diff(Ey,1,1)*Cx);

  % Update E everywhere except on boundary
  Ex(:,2:Ny,2:Nz) = Ex(:,2:Ny,2:Nz) + (Dt /eps0) * ...
      (diff(Hz(:,:,2:Nz),1,2)*Cy - diff(Hy(:,2:Ny,:),1,3)*Cz);
  Ey(2:Nx,:,2:Nz) = Ey(2:Nx,:,2:Nz) + (Dt /eps0) * ...
      (diff(Hx(2:Nx,:,:),1,3)*Cz - diff(Hz(:,:,2:Nz),1,1)*Cx);
  Ez(2:Nx,2:Ny,:) = Ez(2:Nx,2:Ny,:) + (Dt /eps0) * ...
      (diff(Hy(:,2:Ny,:),1,1)*Cx - diff(Hx(2:Nx,:,:),1,2)*Cy);

  % Sample the electric field at chosen points
%   Eto(n,:) = [Ex(4,4,4) Ey(4,4,4) Ez(4,4,4)];
   Eto(n,:) = [Ex(10,10,10) Ey(10,10,10) Ez(10,10,10)];  
end

%===============================================================
Eto1(1:1024)=sqrt(Eto(1:1024,1).^2+Eto(1:1024,2).^2+Eto(1:1024,3).^2);
% plot(Eto(:,1),':');
% hold on
% plot(Eto(:,2),'-.');
% plot(Eto(:,3),'--');
% plot(Eto1,'o');
plot(Eto(:,1));
hold on
plot(Eto(:,2),'m');
plot(Eto(:,3),'g');
plot(Eto1,'r');