%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% fdtd 
% Carlos J. Cela, March 2012
%
% One-dimensional, free space, ABC, FDTD simulation.
% Uses D-H formulation, and frequency-dependent 
% (dispersive) material model (Debye); calculates response at 
% given frequencies at source and arbitrary position 
% in model.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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;
Sx(1:numCells) = 0;

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

% Sigmoid envelope pulse parameters
delay = 10;    
spread = 5; 

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

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

sigma(100:numCells-2) = 0.001;
er(100:numCells-2) =2;
chi1(100:numCells-2) = 2;
t0 = 0.001e-7;

% Frequencies of interest
freq(1) = 50e6;
freq(2) = 200e6;
freq(3) = 500e6;

% Auxiliary calculations
delExp = exp(-dt/t0);
for k=1:numCells
  gbx(k) = sigma(k)*dt/e0;
  gax(k) = 1/(er(k)+gbx(k)+(chi1(k)*dt/t0));
  gbc(k) = chi1(k)*dt/t0;
end

% Fourier components of medium
realPart(1:length(freq),1:numCells-1) = 0;
imagPart(1:length(freq),1:numCells-1) = 0;
amplitude(1:length(freq),1:numCells-1) = 0;
phase(1:length(freq),1:numCells-1) = 0;

% Fourier components of source
realSource(1:length(freq)) = 0;
imagSource(1:length(freq)) = 0;

% Main loop 
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)- delExp*Sx(k));
    Ix(k) = Ix(k)+gbx(k)*Ex(k);
    Sx(k) = delExp*Sx(k) + gbc(k)*Ex(k);
  end

  % Calculate Fourier transform of Ex
  arg = 2*pi*t*dt;
  for k = 1:numCells-1
    for m = 1:length(freq)
      realPart(m, k) = realPart(m, k) + cos(arg*freq(m))* Ex(k);
      imagPart(m, k) = imagPart(m, k) - sin(arg*freq(m))* Ex(k);
    end
  end

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

  % Calculate Fourier transform of input
  if (t<100)
    sourceE = sourceD;%/(e0*er(sourcePos)); 
    for m = 1:length(freq)
      realSource(m) = realSource(m) + cos(arg*freq(m))* sourceE;
      imagSource(m) = imagSource(m) - sin(arg*freq(m))* sourceE;
    end
  end

  % 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(numCells-1);

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

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

  % Plot routines
  if(mod(t,10)==0)
    % First, convert to polar form
    for m = 1:length(freq)
      % Source
      amplitudeSource(m) = sqrt(imagSource(m)^2 + realSource(m)^2);
      phaseSource(m) = atan2(imagSource(m), realSource(m));

      % Everywhere in the model (normalized to source)
      for k = 1:numCells-1
        amplitude(m,k) = 1/amplitudeSource(m)*sqrt(imagPart(m,k)^2 + realPart(m,k)^2);
        phase(m,k) = atan2(imagPart(m,k), realPart(m,k))-phaseSource(m);
      end
    end

    % Now, plot
    subplot (2, 1, 1)  
    plot(1:numCells,(er-1),'k','linewidth',4);
    hold on
    plot(1:numCells, Ex,'r','linewidth',2);
    hold off
    title(['Cell size:' num2str(dz) 'm - Time step:' num2str(dt) 's']);
    grid on
    axis([1 numCells -0.1 1.1]);
    subplot (2, 1, 2)
    plot(amplitude(1,:),'r','linewidth',2);
    hold on
    plot(amplitude(2,:),'b','linewidth',2);
    plot(amplitude(3,:),'k','linewidth',2);
    legend('50MHz', '200MHz', '500MHz');
    %axis([1 numCells 0 1.5]);
    hold off
    grid on;
    
    pause(0.4);
  end
end






