fdtd42.m

Matlab code FDTD 2D ( dennis sullivan book )

  end
end

for i=ia:ib
  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);
        dx(i,j,k)=gj3(j)*gk3(k)*dx(i,j,k)+gj2(j)*gk2(k)*0.5*curl_h;  %review
     end
  end
end

for i=ib+1:ie
  ixh=i-ib;    %review
  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);
        idxh(ixh,j,k)=idxh(ixh,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)*idxh(ixh,j,k));
     end
  end
end

%***********************************************************************
%     Update Dy field
%***********************************************************************
for i=2:ie
  for j=2:ja-1
     for k=2:ke
        curl_h=hx(i,j,k)-hx(i,j,k-1)-hz(i,j,k)+hz(i-1,j,k);
        idyl(i,j,k)=idyl(i,j,k)+curl_h;
        dy(i,j,k)=gi3(i)*gk3(k)*dy(i,j,k)+gi2(i)*gk2(k)*0.5*(curl_h+gj1(j)*idyl(i,j,k));
     end
  end
end

for i=2:ie
  for j=ja:jb
     for k=2:ke
        curl_h=hx(i,j,k)-hx(i,j,k-1)-hz(i,j,k)+hz(i-1,j,k);
        dy(i,j,k)=gi3(i)*gk3(k)*dy(i,j,k)+gi2(i)*gk2(k)*0.5*curl_h;
     end
  end
end

for i=2:ie
   for j=jb+1:je
      jyh=j-jb;
      for k=2:ke
        curl_h=hx(i,j,k)-hx(i,j,k-1)-hz(i,j,k)+hz(i-1,j,k);
        idyh(i,jyh,k)=idyh(i,jyh,k)+curl_h;
        dy(i,j,k)=gi3(i)*gk3(k)*dy(i,j,k)+gi2(i)*gk2(k)*0.5*(curl_h+gj1(j)*idyh(i,jyh,k));
     end
  end
end

%Incident Dy
%***********
for i=ia:ib
   for j=ja:jb-1
      dy(i,j,ka)=dy(i,j,ka)-0.5*hx_inc(j);
      dy(i,j,kb+1)=dy(i,j,kb+1)+0.5*hx_inc(j);
   end
end

%***********************************************************************
%     Update DZ field
%***********************************************************************
for i=2:ie
  for j=2:je
     for k=1:ka-1
        curl_h=hy(i,j,k)-hy(i-1,j,k)-hx(i,j,k)+hx(i,j-1,k);
        idzl(i,j,k)=idzl(i,j,k)+curl_h;
        dz(i,j,k)=gi3(i)*gj3(j)*dz(i,j,k)+gi2(i)*gj2(j)*0.5*(curl_h+gk1(k)*idzl(i,j,k));
     end
  end
end

for i=2:ie
  for j=2:je
     for k=ka:kb
        curl_h=hy(i,j,k)-hy(i-1,j,k)-hx(i,j,k)+hx(i,j-1,k);
        dz(i,j,k)=gi3(i)*gj3(j)*dz(i,j,k)+gi2(i)*gj2(j)*0.5*curl_h;
     end
  end
end

for i=2:ie
   for j=2:je
      for k=kb+1:ke
         kzh=k-kb;
      
         curl_h=hy(i,j,k)-hy(i-1,j,k)-hx(i,j,k)+hx(i,j-1,k);
         idzh(i,j,kzh)=idzh(i,j,kzh)+curl_h;
         dz(i,j,k)=gi3(i)*gj3(j)*dz(i,j,k)+gi2(i)*gj2(j)*0.5*(curl_h+gk1(k)*idzh(i,j,kzh));
     end
  end
end

%Incident DZ
%***********
for i=ia:ib
   for k=ka:kb
      dz(i,ja,k)=dz(i,ja,k)+0.5*hx_inc(ja-1);
      dz(i,jb,k)=dz(i,jb,k)-0.5*hx_inc(jb);
   end
end

 % The source
  
  % pulse=sin(2*pi*400*1e6*dt*T);
  %for k=kc-6:kc+6
  %   dz(ic,jc,k)=0;
  %end
  %pulse=exp(-0.5*(power((to-T)/spread,2)));
  %dz(ic,jc,kc)=pulse;
     
%***********************************************************************
%     Caluculate E fields from D fields
%***********************************************************************

for i=2:ie-1
   for j=2:je-1
      for k=2:ke-1
         ex(i,j,k)=gax(i,j,k)*(dx(i,j,k)-ix(i,j,k));
         ix(i,j,k)=ix(i,j,k)+gbx(i,j,k)*ex(i,j,k);
         
         ey(i,j,k)=gay(i,j,k)*(dy(i,j,k)-iy(i,j,k));
         iy(i,j,k)=iy(i,j,k)+gby(i,j,k)*ey(i,j,k);

         ez(i,j,k)=gaz(i,j,k)*(dz(i,j,k)-iz(i,j,k));
         iz(i,j,k)=iz(i,j,k)+gbz(i,j,k)*ez(i,j,k);
      end
   end
end

%Calculate the fourier transform of Ex
%***************************************
for j=1:je
   for i=1:ie
      for m=1:NFREQS
         real_pt(m,i,j)=real_pt(m,i,j)+cos(arg(m)*T)*ez(i,j,kc);
         imag_pt(m,i,j)=imag_pt(m,i,j)+sin(arg(m)*T)*ez(i,j,kc);
      end
   end
end

% Calculate the incident field
%******************************

for j=1:je-1
   hx_inc(j)=hx_inc(j)+0.5*(ez_inc(j)-ez_inc(j+1));
end


%***********************************************************************
%     Update magnetic fields
%***********************************************************************
                      %1-  Culculate Hx
                      %*****************

for i=1:ia-1
   for j=1:je-1
       for k=1:ke-1
         curl_e=ey(i,j,k+1)-ey(i,j,k)-ez(i,j+1,k)+ez(i,j,k);
         ihxl(i,j,k)=ihxl(i,j,k)+curl_e;
         hx(i,j,k)=fj3(j)*fk3(k)*hx(i,j,k)+fj2(j)*fk2(k)*0.5*(curl_e+fi1(i)*ihxl(i,j,k));

       end
     end
  end
  
  for i=ia:ib
   for j=1:je-1
       for k=1:ke-1
         curl_e=ey(i,j,k+1)-ey(i,j,k)-ez(i,j+1,k)+ez(i,j,k);
         hx(i,j,k)=fj3(j)*fk3(k)*hx(i,j,k)+fj2(j)*fk2(k)*0.5*curl_e;
        end
     end
  end
  
for i=ib+1:ie
  ixh=i-ib;
  for j=1:je-1
     for k=1:ke-1
         curl_e=ey(i,j,k+1)-ey(i,j,k)-ez(i,j+1,k)+ez(i,j,k);
        ihxh(ixh,j,k)=ihxh(ixh,j,k)+curl_e;
        hx(i,j,k)=fj3(j)*fk3(k)*hx(i,j,k)+fj2(j)*fk2(k)*0.5*(curl_e+fi1(i)*ihxh(ixh,j,k));


     end
  end
end

     %Incident Hx value
     %&&&&&&&&&&&&&&&&&&
     for i=ia:ib
        for k=ka:kb
           hx(i,ja-1,k)= hx(i,ja-1,k)+0.5*ez_inc(ja);
           hx(i,jb,k)  = hx(i,jb,k)-0.5*ez_inc(jb);
        end
     end
     

%***********************************************************************
%     Update Hy field
%***********************************************************************
for i=1:ie-1
  for j=1:ja-1
     for k=1:ke-1
        curl_e=ez(i+1,j,k)-ez(i,j,k)-ex(i,j,k+1)+ex(i,j,k);
        ihyl(i,j,k)=ihyl(i,j,k)+curl_e;
        hy(i,j,k)=fi3(i)*fk3(k)*hy(i,j,k)+fi2(i)*fk3(k)*0.5*(curl_e+fj1(j)*ihyl(i,j,k));  
     end
  end
end

for i=1:ie-1
  for j=ja:jb
     for k=1:ke-1
        curl_e=ez(i+1,j,k)-ez(i,j,k)-ex(i,j,k+1)+ex(i,j,k);
        hy(i,j,k)=fi3(i)*fk3(k)*hy(i,j,k)+fi2(i)*fk3(k)*0.5*curl_e;
     end
  end
end

for i=1:ie-1
   for j=jb+1:je
      jyh=j-jb;
      for k=1:ke-1
        curl_e=ez(i+1,j,k)-ez(i,j,k)-ex(i,j,k+1)+ex(i,j,k);
        ihyh(i,jyh,k)=ihyh(i,jyh,k)+curl_e;
        hy(i,j,k)=fi3(i)*fk3(k)*hy(i,j,k)+fi2(i)*fk3(k)*0.5*(curl_e+fj1(j)*ihyh(i,jyh,k));
     end
  end
end

%Incident Hy
%***********
for j=ja:jb
   for k=ka:kb
      hy(ia-1,j,k)=hy(ia-1,j,k)-0.5*ez_inc(j);
      hy(ib,j,k)=hy(ib,j,k)+0.5*ez_inc(j);
   end
end

%***********************************************************************
%     Update HZ field
%***********************************************************************
for i=1:ie-1
  for j=1:je-1
     for k=1:ka-1
        curl_e=ex(i,j+1,k)-ex(i,j,k)-ey(i+1,j,k)+ey(i,j,k);
        ihzl(i,j,k)=ihzl(i,j,k)+curl_e;
        hz(i,j,k)=fi3(i)*fj3(j)*hz(i,j,k)+fi2(i)*fj2(j)*0.5*(curl_e+fk1(k)*ihzl(i,j,k));
     end
  end
end

for i=1:ie-1
  for j=1:je-1
     for k=ka:kb
        curl_e=ex(i,j+1,k)-ex(i,j,k)-ey(i+1,j,k)+ey(i,j,k);
        hz(i,j,k)=fi3(i)*fj3(j)*hz(i,j,k)+fi2(i)*fj2(j)*0.5*curl_e;
     end
  end
end

for i=1:ie-1
   for j=1:je-1
      for k=kb+1:ke
         kzh=k-kb;
      
         curl_e=ex(i,j+1,k)-ex(i,j,k)-ey(i+1,j,k)+ey(i,j,k);

         ihzh(i,j,kzh)=ihzh(i,j,kzh)+curl_e;
         hz(i,j,k)=fi3(i)*fj3(j)*hz(i,j,k)+fi2(i)*fj2(j)*0.5*(curl_e+fk1(k)*ihzh(i,j,kzh));

     end
  end
end

%saving ez at different time
      if T==30
         ez30=ez;
      elseif T==40
          ez40=ez;
      elseif T==50
          ez50=ez;
      elseif T==60
          ez60=ez;
      end
%end saving

%*************************************************************************     
% End of the main loop
%*************************************************************************     
end

%Calculate the Fourier amplitudee and phase of the incident pulse
%*****************************************************************
for m=1:NFREQS
   amp_in(m)=sqrt(power(real_in(m),2)+power(imag_in(m),2));
   phase_in(m)=atan2(imag_in(m),real_in(m)); 
  end
  
%Calculate the Fourier amplitudee and phase of the total field
%**************************************************************
 for m=1:NFREQS
    for j=ja:jb
       amp(ic,j)=(1./amp_in(m))* sqrt(power(real_pt(m,ic,j),2)+power(imag_pt(m,ic,j),2));
    end
 end


%Plotting the electric field at various time 

     figure
     subplot(4,2,1),surfc(ez30(:,:,30));
     title(['T = 30']);
     xlabel('cm');
     ylabel('cm');
     zlabel('EZ');

    subplot(4,2,2),surfc(ez40(:,:,30));
    title(['T = 40']);
    ylabel('cm');
    zlabel('EZ');

     
    subplot(4,2,3),surfc(ez50(:,:,30));
    title(['T = 50']);
    ylabel('EZ');
    xlabel('cm');
    zlabel('EZ');

    
    subplot(4,2,4),surfc(ez60(:,:,30));
    title(['T = 60']);
    ylabel('cm');
    xlabel('cm');
    zlabel('EZ');


%*****************************************************************************