fdtdfor.m

A SIMPLE BUT USFUL 2D FDTD MATLAB CODE

% 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 = 32;                      % 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
 
  for i = 1:Nx+1
     for j = 1:Ny
        for k = 1:Nz
           Hx(i,j,k) = Hx(i,j,k) + (Dt/mu0)* ...
           ((Ey(i,j,k+1)-Ey(i,j,k))*Cz - (Ez(i,j+1,k)-Ez(i,j,k))*Cy);
        end
     end
  end
  
  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

  for i = 1:Nx
     for j = 2:Ny
        for k = 2:Nz
           Ex(i,j,k) = Ex(i,j,k) + (Dt /eps0) * ...
              ((Hz(i,j,k)-Hz(i,j-1,k))*Cy-(Hy(i,j,k)-Hy(i,j,k-1))*Cz);
        end
     end
  end

  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
  Etss(n,:) = [Ex(10,10,10) Ey(10,10,10) Ez(10,10,10)];
end

%  ===============================================================
Etss1(1:32)=sqrt(Etss(1:32,1).^2+Etss(1:32,2).^2+Etss(1:32,3).^2);
% plot(Etss(:,1),':');
% hold on
% plot(Etss(:,2),'-.');
% plot(Etss(:,3),'--');
% plot(Etss1,'o');
plot(Etss(:,1));
hold on
plot(Etss(:,2),'m');
plot(Etss(:,3),'g');
plot(Etss1,'r');