fdtdnodiff.m

A SIMPLE BUT USFUL 2D FDTD MATLAB CODE2

% 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)*((Ey(:,:,2:Nz+1)-Ey(:,:,1:Nz))*Cz ...
                    - (Ez(:,2:Ny+1,:)-Ez(:,1:Ny,:))*Cy);
  Hy = Hy + (Dt/mu0)*((Ez(2:Nx+1,:,:)-Ez(1:Nx,:,:))*Cx ...
                    - (Ex(:,:,2:Nz+1)-Ex(:,:,1:Nz))*Cz);
  Hz = Hz + (Dt/mu0)*((Ex(:,2:Ny+1,:)-Ex(:,1:Ny,:))*Cy ...
                    - (Ey(2:Nx+1,:,:)-Ey(1:Nx,:,:))*Cx);

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

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