fdtd42.m

Matlab code FDTD 2D ( dennis sullivan book )

clc;clear;
%**********************************************************************
%                      fdtd4-2c program
%**********************************************************************
%Three dimensional FDTD program of a plane wave impinging on a 
%dielectric sphere 
%The source is in position at (ic,jc,kc) cell of the problem space with
%PML absorbing boundary conditions.
%**********************************************************************

%***********************************************************************
%     Fundamental constants
%***********************************************************************
pi   = 3.14159;
ddx  = 0.01 ;               % cell size
epsz = 8.85e-12;            %permittivity of free space
muz  = 4.0*pi*1.0e-7;          %permeability of free space

%***********************************************************************
%     Grid parameters
%***********************************************************************

ie   = 40;                     % number of grid cells in X direction.
je   = 40;                     % number of grid cells in Y-direction. 
ke   = 40;                     % number of grid cells in Z-direction. 
nmax=65;
ddx  = 0.01;
dt   = ddx/(2*3e8);		 %The center of the incident pulse

%Total/Scattered field boundaries.
%********************************
ia=7;
ja=7;
ka=7;

ib=ie-ia-1;
jb=je-ja-1;
kb=ke-ka-1;

npml = 9;                %The number of PML cells

%***********************************************************************
%     Initialize Field arrays
%***********************************************************************
ex(1:ie,1:je,1:ke) = 0.0;
ey(1:ie,1:je,1:ke) = 0.0;
ez(1:ie,1:je,1:ke) = 0.0;

dx(1:ie,1:je,1:ke) = 0.0;
dy(1:ie,1:je,1:ke) = 0.0;
dz(1:ie,1:je,1:ke) = 0.0;

hx(1:ie,1:je,1:ke) = 0.0;
hy(1:ie,1:je,1:ke) = 0.0;
hz(1:ie,1:je,1:ke) = 0.0;

gax(1:ie,1:je,1:ke) = 1.0;
gay(1:ie,1:je,1:ke) = 1.0;
gaz(1:ie,1:je,1:ke) = 1.0;

gbx(1:ie,1:je,1:ke) = 0.0;
gby(1:ie,1:je,1:ke) = 0.0;
gbz(1:ie,1:je,1:ke) = 0.0;

ix(1:ie,1:je,1:ke) = 0.0;
iy(1:ie,1:je,1:ke) = 0.0;
iz(1:ie,1:je,1:ke) = 0.0;

ez_inc(1:je+1)   =0.0;
hx_inc(1:je+1)   =0.0;
ez_low_m2  =0.0;
ez_high_m2 =0.0;
ez_low_m1  =0.0;
ez_high_m1 =0.0;

%************************************************************************
%Fourier analysis parameters
%************************************************************************

NFREQS=3;
real_in(1:NFREQS)=0.0;
imag_in(1:NFREQS)=0.0;

real_pt(1:NFREQS,1:ie,1:je)=0.0;
imag_pt(1:NFREQS,1:ie,1:je)=0.0;

amp(1:ie,1:je)=0.0;
phase(1:ie,1:je)=0.0;

amp_in(1:NFREQS)=0.0;
phase_in(1:NFREQS)=0.0;

freq(1)=10e6;
freq(2)=100e6;
freq(3)=433e6;

for n=1:NFREQS
   arg(n)=2*pi*freq(n)*dt;
end

idxl(1:ia+1,1:je,1:ke)=0.0; idxh(1:ia+1,1:je,1:ke)=0.0;  %review
ihxl(1:ia+1,1:je,1:ke)=0.0; ihxh(1:ia+1,1:je,1:ke)=0.0; 
idyl(1:ie,1:ja+1,1:ke)=0.0; idyh(1:ie,1:ja+1,1:ke)=0.0; 
ihyl(1:ie,1:ja+1,1:ke)=0.0; ihyh(1:ie,1:ja+1,1:ke)=0.0; 
idzl(1:ie,1:je,1:ka+1)=0.0; idzh(1:ie,1:je,1:ka+1)=0.0; 
ihzl(1:ie,1:je,1:ka+1)=0.0; ihzh(1:ie,1:je,1:ka+1)=0.0; 


%***********************************************************************
%     Specify the dipole Source 
%***********************************************************************
ic  = ie/2; jc=je/2; kc=ke/2;      %The position of the source

to      = 40;
spread  = 10;		         %The width of the incident pulse
T       = 0 ;

%*************************************************************************
% Specify the dielectric sphere
%*************************************************************************
radius=10;
epsilon=[1.0 30];
sigma =[0.0 0.3];

%Calculate the gax and gbx
%*************************

for i=ia:ib-1
   for j=ja:jb-1
      for k=ka:kb-1
         eps=epsilon(1);
         cond=sigma(1);
         xdist=(ic-i-0.5);
         ydist=(jc-j);
         zdist=(kc-k);
         dist=sqrt(power(xdist,2)+power(ydist,2)+power(zdist,2));
           if dist<=radius
               eps=epsilon(2);
               cond=sigma(2);
            end
         gax(i,j,k)=1.0/(eps+(cond*dt/epsz));
         gbx(i,j,k)=cond*dt/epsz;
      end
   end
end

%Calculate the gay and gby
%*************************

for i=ia:ib-1
   for j=ja:jb-1
      for k=ka:kb-1
         eps=epsilon(1);
         cond=sigma(1);
         xdist=(ic-i);
         ydist=(jc-j-0.5);
         zdist=(kc-k);
         dist=sqrt(power(xdist,2)+power(ydist,2)+power(zdist,2));
           if dist<=radius
               eps=epsilon(2);
               cond=sigma(2);
            end
         gay(i,j,k)=1.0/(eps+(cond*dt/epsz));
         gby(i,j,k)=cond*dt/epsz;
      end
   end
end

%Calculate the gaz and gbz
%*************************

for i=ia:ib-1
   for j=ja:jb-1
      for k=ka:kb-1
         eps=epsilon(1);
         cond=sigma(1);
         xdist=(ic-i);
         ydist=(jc-j);
         zdist=(kc-k-0.5);
         dist=sqrt(power(xdist,2)+power(ydist,2)+power(zdist,2));
           if dist<=radius
               eps=epsilon(2);
               cond=sigma(2);
            end
         gaz(i,j,k)=1.0/(eps+(cond*dt/epsz));
         gbz(i,j,k)=cond*dt/epsz;
      end
   end
end

%***********************************************************************
%     Calculate the PML parameters
%***********************************************************************
for i=1:ie
   gi1(i)=0.0;
   fi1(i)=0.0;
   gi2(i)=1.0;
   fi2(i)=1.0;
   gi3(i)=1.0;
   fi3(i)=1.0;
 end

for j=1:je
   gj1(j)=0.0;
   fj1(j)=0.0;
   gj2(j)=1.0;
   fj2(j)=1.0;
   gj3(j)=1.0;
   fj3(j)=1.0;
end

for k=1:ke
   gk1(k)=0.0;
   fk1(k)=0.0;
   gk2(k)=1.0;
   fk2(k)=1.0;
   gk3(k)=1.0;
   fk3(k)=1.0;
end

for i=1:npml
   xxn=(npml-i)/npml;
   xn=0.33*power(xxn,3);
   
   fi1(i)     =xn;
   fi1(ie-1-i)=xn;

   gi2(i)     =1.0/(1.0+xn);
   gi2(ie-1-i)=1.0/(1.0+xn);
   
   gi3(i)     =(1.0-xn)/(1.0+xn);
   gi3(ie-1-i)=(1.0-xn)/(1.0+xn);
   
   xxn=(npml-i-0.5)/npml;
   xn=0.25*power(xxn,3);
   
   gi1(i) =xn;
   gi1(ie-2-i)=xn;
   
   fi2(i)     =1.0/(1.0+xn);
   fi2(ie-2-i)=1.0/(1.0+xn);
   
   fi3(i)     =(1.0-xn)/(1.0+xn);
   fi3(ie-2-i)=(1.0-xn)/(1.0+xn);
end

for j=1:npml
   xnum=npml-j;
   xd=npml;
   xxn=xnum/xd;
   xn=0.33*power(xxn,3);
   
   fj1(j)     =xn;
   fj1(je-1-j)=xn;

   gj2(j)     =1.0/(1.0+xn);
   gj2(je-1-j)=1.0/(1.0+xn);
   
   gj3(j)     =(1.0-xn)/(1.0+xn);
   gj3(je-1-j)=(1.0-xn)/(1.0+xn);
   
   xxn=(xnum-0.5)/xd;
   xn=0.33*power(xxn,3);
   
   gj1(j)     =xn;
   gj1(je-2-j)=xn;

   fj2(j)     =1.0/(1.0+xn);
   fj2(je-2-j)=1.0/(1.0+xn);
   
   fj3(j)     =(1.0-xn)/(1.0+xn);
   fj3(je-2-j)=(1.0-xn)/(1.0+xn);
end

for k=1:npml
   xnum=npml-k;
   xd=npml;
   xxn=xnum/xd;
   xn=0.33*power(xxn,3);
   
   fk1(k)     =xn;
   fk1(ke-1-k)=xn;

   gk2(k)     =1.0/(1.0+xn);
   gk2(ke-1-k)=1.0/(1.0+xn);
   
   gk3(k)     =(1.0-xn)/(1.0+xn);
   gk3(ke-1-k)=(1.0-xn)/(1.0+xn);
   
   xxn=(xnum-0.5)/xd;
   xn=0.33*power(xxn,3);
   
   gk1(k)     =xn;
   gk1(ke-2-k)=xn;

   fk2(k)     =1.0/(1.0+xn);
   fk2(ke-2-k)=1.0/(1.0+xn);
   
   fk3(k)     =(1.0-xn)/(1.0+xn);
   fk3(ke-2-k)=(1.0-xn)/(1.0+xn);
end

%***********************************************************************
%     BEGIN TIME-STEPPING LOOP
%***********************************************************************

for n=1:nmax
   
   T=T+1
   
   %Calculate the incident buffer
   %*****************************
   
   for j=2:je
      ez_inc(j)=ez_inc(j)+0.5*(hx_inc(j-1)-hx_inc(j));
   end
   
   %Fourier transform of the incident field
   %****************************************
   for m=1:NFREQS
      real_in(m)=real_in(m)+cos(arg(m)*T)*ez_inc(ja-1);
      imag_in(m)=imag_in(m)-sin(arg(m)*T)*ez_inc(ja-1);
   end
   
  %Puts the sinusoidal source in the middle
  %  pulse=sin(2*pi*400*1e6*dt*T);
  
  pulse=exp(-0.5*(power((to-T)/spread,2)));
  ez_inc(3)=pulse;

   %ABC for the incident buffer
   %&&&&&&&&&&&&&&&&&&&&&&&&&&&&
   ez_inc(1)=ez_low_m2;
   ez_low_m2=ez_low_m1;
   ez_low_m1=ez_inc(2);
   
   ez_inc(je)=ez_high_m2;
   ez_high_m2=ez_high_m1;
   ez_high_m1=ez_inc(je-1);
   
%***********************************************************************
%     Update Dx field
%***********************************************************************
for i=2:ia-1
  for j=2:je
     for k=2:ke
        curl_h=hz(i,j,k)-hz(i,j-1,k)-hy(i,j,k)+hy(i,j,k-1);
        idxl(i,j,k)=idxl(i,j,k)+curl_h;
        dx(i,j,k)=gj3(j)*gk3(k)*dx(i,j,k)+gj2(j)*gk2(k)*0.5*(curl_h+gi1(i)*idxl(i,j,k));
     end