1 D FDTD with dielectric medium and absorbing boundary

%%%%%%%%%%%%%%%%%%%%%
close all;
clear all;

Nz = 200;                    %// Number of cells in z-direction
Nt = 500;                    %// Number of time steps

N = Nz;                      %// Global number of cells

Ex         = zeros(Nz,1);      %// Ex component
Hy         = zeros(Nz,1);      %// Hy component
eps        = ones(1,Nz);       %// Relative permittivity 
sig        = ones(1,Nz);
mu         = ones(1,Nz);
ca          = ones(1,Nz);
cb          = ones(1,Nz);


cc = 2.99792458e8;          %// Speed of light in free space
muz = 4.0*pi*1.0e-7;        %// Permeability of free space
epsz = 1.0/(cc*cc*muz);     %// Permittivty of free space

freq = 10e6;               %// Excitation frequency
omega = 2.0*pi*freq;       %// Excitation circular frequency
lambda = cc/freq;   %// Reference wavelength

eps(1:N-120)= 1.0;
eps(N-120:N)=4.0;

sig(1:N)=0.0;

Dx = lambda/10;            %// Reference cells width
Z = Nz*Dx;
% Dt = Dx/cc;                %// Reference time step
Dt=1*(Dx/(cc*1/sqrt(min(eps))));


rbc=(cc*Dt-Dx)/(cc*Dt-Dx);
scfact=Dt/muz/Dx;

Ex(1:N)=0.0;
Hy(1:N)=0.0;

for nt=1:Nt
     
    t = (Dt*(nt-1))/2;       %// Actual time
    
   Ex(10)= sin(omega*t);
   
      for n=2:N
      Hy(n) = Hy(n) -( Ex(n) - Ex(n-1) ); 
      end
    
                
    for n=2:N-1 
        eaf=(Dt*sig(n))/(2.0*epsz*eps(n));
        ca=(1.0-eaf)/(1.0+eaf);
        cb=scfact*(Dt/(epsz*eps(n)*Dx))/(1.0+eaf);
  
        if n==1
        Ex(n) = ca * Ex(n) - cb *( Hy(n+1) - Hy(n) );
        else
          if n==2           %MUR ABC for z=0, writes to E(1)
           Ex(n-1)=Ex(n)+ rbc*(Ex(n)+Ex(n-1));
           Ex(n-1)=Ex(n);     
          elseif n==N-1     %MUR ABC for z=zmax, writes to E(N)
           Ex(n+1)=Ex(n)+ rbc*(Ex(n)+Ex(n+1));
           Ex(n+1)=Ex(n);
          end
         end
          Ex(n) = ca * Ex(n) - cb *( Hy(n+1) - Hy(n) );    
    end
    

    %// Generate output figures  
    subplot(2,1,1)
    plot(Ex,'b')
    axis([1 Nz -2.5 2.5]);   
    xlabel('z in m');    
    ylabel('Amplitude of E_x');
    grid on 

    subplot(2,1,2)
    plot(Hy,'r')
    axis([1 Nz -2.5 2.5]);   
    xlabel('z in m');    
    ylabel('Amplitude of H_y');
    grid on   
    
    pause(0.0005)
    getframe;
   
end %// End of marching-on-in-time loop