%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% fdtd
% Carlos J. Cela, March 2012
%
% One-dimensional, free space FDTD simulation
% using a Gaussian source (D-H formulation).
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all

% constants
c0 = 3e8;
e0 =8.8541e-12;
u0 =4*pi*1e-7;

% Number of cells in Z direction
numCells = 200;

% Initialize fields to zero
Ex(1:numCells) = 0;
Hy(1:numCells) = 0;
Ix(1:numCells) = 0;
Dx(1:numCells) = 0;

% Parameters
dz = 0.01;
dt = dz/(2*c0);

% Sigmoid envelope pulse parameters
t0 = 50;    
spread = 20; 

% Boundary conditions
Eleft(3) = 0;
Eright(2) = 0;

% Material properties
sigma(1:numCells) = 0;
er(1:numCells) =1;

er(100:120) = 5;
sigma(100:120) = 0.1;

% Auxiliary calculations
for k=1:numCells
  gbx(k) = sigma(k)*dt/e0;
  gax(k) = 1/(er(k)+gbx(k));
end

% Main loop (Loop C)
for t = 1:1500

  % Electric field
  for k = 2:numCells
    Dx(k) = Dx(k) + 0.5*(Hy(k-1)-Hy(k));
    Ex(k) = gax(k)*(Dx(k)-Ix(k));
    Ix(k) = Ix(k)+gbx(k)*Ex(k);
  end

  % Add E field source now
  Dx(2) = Dx(2)+exp(-0.5*((t0-t)/spread)^2);

   % Boundary conditions
  Eleft(3) = Eleft(2);
  Eleft(2) = Eleft(1);
  Eleft(1) = Ex(2);

  Eright(3) = Eright(2);
  Eright(2) = Eright(1);
  Eright(1) = Ex(199);

  Ex(1) = Eleft(3);
  Ex(200) = Eright(3);

  % Magnetic field
  for k = 1:numCells-1
    Hy(k) = Hy(k) + 0.5*(Ex(k)-Ex(k+1));
  end

  % Plot  
  plot(1:numCells,(er-1)/3,'linewidth',4);
  hold on
  plot(1:numCells, Ex,'r','linewidth',2);
  hold off
  axis([1 200 -2 2]);
  grid on;
  title(['Cell size:' num2str(dz) 'm - Time step:' num2str(dt) 's']);
  pause(0.00001);
end






