% Title: FDTD1D
% One dimensional FDTD code
% TEM mode with dielectric transition at x=BB
% PMC left boundary x=0
% PEC right boundary x=xN; 
% David Thiel - 18 September 2000; revised May 2003
clear; hold off;
xN=100; 		% Set the solution space size in x-direction
H=zeros(xN,1);  % Set initial electric field components = zero	
E=zeros(xN,1);	% Set initial magnetic field components = zero
eta=ones(1,xN);	% Set material properties (impedance) arrays
BB=80;					% position of dielectric boundary
nmax=150;				% maximum number of time steps
delX=1e-3;				% Cell size in x-direction
ST=40;					% Source location (cell number in x-direction)
%Set material parameters of air and dielectric
e0=8.854e-12;	% Permittivity of air (free space)
m0=pi*4e-7;	    % Permeability of air (free space)
e1=4.9*e0;	    % Permittivity of dielectric (e.g., PCB)
m1=m0;			% Permeability of dielectric (e.g., PCB) 
c=1/sqrt(e0*m0);	     % Velocity of propagation in dielectric
delT=delX/(c*sqrt(1)); 	% Courant condition for minimum time step
eta0=sqrt(m0/e0);	    % Instrinsic impedance of air (free space)
eta1=sqrt(m1/e1);	    % Instrinsic impedance of dielectric (e.g., PCB)
eta(1:BB)=eta0*eta(1:BB);		% Set impedance of air in 
                        % material property array
eta(BB+1:xN)=eta1*eta(BB+1:xN);	% Set impedance of dielectric in 
                        % material property array
% Set source parameters for Gaussian pulse
ndecay=3;		
n0=4*ndecay;

% The time loop
for nt=1:nmax-1
   T(nt)=nt;
      Eactive=exp(-((nt-n0)/ndecay)^2); % Gaussian pulse source
      E(ST,nt)= Eactive+E(ST,nt);	% Insert pulse at x=20
% Update equation for E field 
   for xi=2:xN
   E(xi,nt+1)=(H(xi,nt)-H(xi-1,nt))*eta(xi)+E(xi,nt);
	end
% Reset the PEC boundary conditions
      E(xN,nt+1)=0;		
% Update equation for H field
   for xi=1:xN-1
	H(xi,nt+1)=(E(xi+1,nt+1)-E(xi,nt+1))/eta(xi)+H(xi,nt);
	end
% Reset the PMC boundary condition
    H(1,nt+1)=0;
end		% End of time loop

X=delX*[1:xN];	% This is the X axis array
% Now set up the movie
hl=delX*[BB BB];	% Vertical boundary line on plot
vl=[-1.2 1.2];		% Vertical boundary line on plot
B='Frame number';
for nt=1:nmax
   plot(hl,vl,X,E(:,nt)); % plot E field and boundary line
	axis([0 xN*delX -1.2 1.2]); % fix the axes
	xlabel('Distance (m)');
	ylabel('E field (V/m)');
   A=num2str(nt);
   Fno=strcat(B,A); % this is the frame number
   title(Fno);
   M(nt)=getframe;
end
movie(M);	% Now run the movie