untitled.m

FDTD CPML 計算介質端平均電場功率

   for i = 1:ihT
      for j = 1:jhT 
            CA(i,j) = (1.0 - sig(i,j)*dt / (2.0*epsi(i,j))) / (1.0 + sig(i,j) * dt / (2.0*epsi(i,j)));
            CB(i,j) = (dt/(epsi(i,j))) / (1.0 + sig(i,j)*dt / (2.0*epsi(i,j)));
      end
    end
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  FILL IN DENOMINATORS FOR FIELD UPDATES
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  CPML, Psi_w,v
% (7.109)
   ii = iBC-1;
   for i = 1:ieT
      if i <= iBC-1
         F_hx(i) = 1.0/(K_HxBC(i)*dxyz);
      elseif i >= ihT+1-iBC
         F_hx(i) = 1.0/(K_HxBC(ii)*dxyz);
         ii = ii-1;
      else
         F_hx(i) = 1.0/dxyz;
      end
   end
   jj = jBC-1;
   for j = 1:jeT;
      if j <= jBC-1 
         F_hy(j) = 1.0/(K_HyBC(j)*dxyz);
      elseif j >= jhT+1-jBC
         F_hy(j) = 1.0/(K_HyBC(jj)*dxyz);
         jj = jj-1;
      else
         F_hy(j) = 1.0/dxyz;
      end
   end 
   ii = iBC;
   for i = 1:ieT
      if i <= iBC 
         F_ex(i) = 1.0/(K_ExBC(i)*dxyz);
      elseif i >= ihT+1-iBC
         F_ex(i) = 1.0/(K_ExBC(ii)*dxyz);
         ii = ii-1;
      else
         F_ex(i) = 1.0/dxyz;
      end
   end
   jj = jBC;
   for j = 1:jeT
      if j <= jBC
         F_ey(j) = 1.0/(K_EyBC(j)*dxyz);
      elseif j >= jhT+1-jBC
         F_ey(j) = 1.0/(K_EyBC(jj)*dxyz);
         jj = jj-1;
      else
         F_ey(j) = 1.0/dxyz;
      end
   end 

%.:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:.
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  BEGIN TIME STEP
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
  for nTimeStep = 1:MaxTime 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  UPDATE Ex
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
      for i = 1:ieT
	    for j = 2:jeT 
	         Ex(i,j) = CA(i,j) * Ex(i,j) + CB(i,j) * ( (Hz(i,j) - Hz(i,j-1))*F_ey(j) ); % The first four terms in (7.106)
	    end
      end 
      for i = 1:ieT
%.....................................................................
%  CPML for bottom Ex, j-direction
%.....................................................................
         for j = 2:jBC
  	         psi_Exy1(i,j) = bEy(j)*psi_Exy1(i,j) + cEy(j) *(Hz(i,j) - Hz(i,j-1))/dxyz;  %(7.105)
             Ex(i,j) = Ex(i,j) + CB(i,j)*psi_Exy1(i,j); %The last two terms in (7.106)
         end 
      end  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  UPDATE Ey
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
      for i = 2:ieT
	     for j = 1:jeT 
             Ey(i,j) = CA(i,j) * Ey(i,j) + CB(i,j) * ( (Hz(i-1,j) - Hz(i,j))*F_ex(i) ); 
         end 
       end
      for j = 1:jeT
%.....................................................................
%  CPML for bottom Ey, i-direction
%.....................................................................

         for i = 2:iBC
	      psi_Eyx1(i,j) = bEx(i)*psi_Eyx1(i,j) + cEx(i)*(Hz(i-1,j) - Hz(i,j))/dxyz;
	      Ey(i,j) = Ey(i,j) + CB(i,j)*psi_Eyx1(i,j);
         end
%.....................................................................
%  CPML for top Ey, i-direction
%.....................................................................
         ii = iBC;
         for i = ihT+1-iBC:ieT
	        psi_Eyx2(ii,j) = bEx(ii)*psi_Eyx2(ii,j) + cEx(ii)*(Hz(i-1,j) - Hz(i,j))/dxyz;
	        Ey(i,j) = Ey(i,j) + CB(i,j)*psi_Eyx2(ii,j);
            ii = ii-1;
	     end
      end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  UPDATE Ez
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
      for i = 2:ieT
         for j = 2:jeT
            Ez(i,j) = CA(i,j) * Ez(i,j) + CB(i,j) * ((Hy(i,j) - Hy(i-1,j))*F_ex(i) + (Hx(i,j-1) - Hx(i,j))*F_ey(j));
	     end
      end
      for j = 2:jeT
%.....................................................................
%  CPML for bottom Ez, x-direction
%.....................................................................
         for i = 2:iBC
   	        psi_Ezx1(i,j) = bEx(i)*psi_Ezx1(i,j) + cEx(i) *(Hy(i,j) - Hy(i-1,j))/dxyz;
	        Ez(i,j) = Ez(i,j) + CB(i,j)*psi_Ezx1(i,j);
         end
%.....................................................................
%  CPML for top Ez, x-direction
%.....................................................................
         ii = iBC;
         for i = ihT+1-iBC:ieT
   	         psi_Ezx2(ii,j) = bEx(ii)*psi_Ezx2(ii,j) + cEx(ii) *(Hy(i,j) - Hy(i-1,j))/dxyz;
	         Ez(i,j) = Ez(i,j) + CB(i,j)*psi_Ezx2(ii,j);
             ii = ii-1;
	     end
      end
      for i = 2:ieT

%.....................................................................
%  CPML for bottom Ez, y-direction
%.....................................................................
         for j = 2:jBC
            psi_Ezy1(i,j) = bEy(j)*psi_Ezy1(i,j) + cEy(j)*(Hx(i,j-1) - Hx(i,j))/dxyz;
	        Ez(i,j) = Ez(i,j) + CB(i,j)*psi_Ezy1(i,j);
         end
%.....................................................................
%  CPML for top Ez, y-direction
%.....................................................................
         jj = jBC;
         for j = jhT+1-jBC:jeT
             psi_Ezy2(i,jj) = bEy(jj)*psi_Ezy2(i,jj) + cEy(jj)*(Hx(i,j-1) -  Hx(i,j))/dxyz;
	         Ez(i,j) = Ez(i,j) + CB(i,j)*psi_Ezy2(i,jj);
            jj = jj-1;
	     end
      end 

%-----------------------------------------------------------------------
%   SOURCE
%-----------------------------------------------------------------------
%   i = isource;
%   j = jsource;
%   k = ksource;
  % source = -2.0*((nTimeStep*dt-tO)/tw) *
  % exp(-((nTimeStep*dt-tO)/tw)^2.0); 
    %source =exp(-((nTimeStep*dt-t0)/tw)^2);
    %if nTimeStep<=floor(t_periodic)
     %  source = sin(omega*dt*nTimeStep); 
    %else
     %   source = 0;
    %end
    source =sin(omega*nTimeStep*dt);
    
    Ez(isrc,30) = Ez(isrc,30)+source;

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  UPDATE Hx
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      for i = 2:ieT
	   for j = 1:jeT
	      Hx(i,j) = DA * Hx(i,j) + DB * ( (Ez(i,j) - Ez(i,j+1))*F_hy(j) ) ;
	   end
      end
      for i = 2:ieT
%.....................................................................
%  CPML for bottom Hx, j-direction
%.....................................................................
         for j = 1:jBC-1
  	         psi_Hxy1(i,j) = bHy(j)*psi_Hxy1(i,j) + cHy(j) *(Ez(i,j) - Ez(i,j+1))/dxyz;
             Hx(i,j) = Hx(i,j) + DB*psi_Hxy1(i,j);
         end
%.....................................................................
%  CPML for top Hx, j-direction
%.....................................................................
         jj = jBC-1;
         for j = jhT+1-jBC:jeT
  	         psi_Hxy2(i,jj) = bHy(jj)*psi_Hxy2(i,jj) + cHy(jj) *(Ez(i,j) - Ez(i,j+1))/dxyz ;
             Hx(i,j) = Hx(i,j) + DB*psi_Hxy2(i,jj) ;
             jj = jj-1 ;
         end
      end 
%  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  UPDATE Hy
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
      for i = 1:ieT
	   for j = 2:jeT
            Hy(i,j) = DA * Hy(i,j) + DB * ( (Ez(i+1,j) - Ez(i,j))*F_hx(i) );
	   end 
      end
      for j = 2:jeT
%.....................................................................
%  CPML for bottom Hy, i-direction
%.....................................................................
         for i = 1:iBC-1
	        psi_Hyx1(i,j) = bHx(i)*psi_Hyx1(i,j) + cHx(i)*(Ez(i+1,j) - Ez(i,j))/dxyz;
	        Hy(i,j) = Hy(i,j) + DB*psi_Hyx1(i,j);
         end
%.....................................................................
%  CPML for top Hy, i-direction
%.....................................................................
         ii = iBC-1;
         for i = ihT+1-iBC:ieT
	         psi_Hyx2(ii,j) = bHx(ii)*psi_Hyx2(ii,j) + cHx(ii)*(Ez(i+1,j) - Ez(i,j))/dxyz;
	         Hy(i,j) = Hy(i,j) + DB*psi_Hyx2(ii,j);
             ii = ii-1;
         end
     end
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  UPDATE Hz
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      for i = 1:ieT
         for j = 1:jeT
            Hz(i,j) = DA * Hz(i,j) + DB * ((Ey(i,j) - Ey(i+1,j))*F_hx(i) + (Ex(i,j+1) - Ex(i,j))*F_hy(j));
	     end
      end
      for j = 1:jeT
%.....................................................................
%  CPML for bottom Hz, x-direction
%.....................................................................
         for i = 1:iBC-1
   	        psi_Hzx1(i,j) = bHx(i)*psi_Hzx1(i,j) + cHx(i) *(Ey(i,j) - Ey(i+1,j))/dxyz;
	        Hz(i,j) = Hz(i,j) + DB*psi_Hzx1(i,j);
         end
%.....................................................................
%  CPML for top Hz, x-direction
%.....................................................................
         ii = iBC-1;
         for i = ihT+1-iBC:ieT
   	         psi_Hzx2(ii,j) = bHx(ii)*psi_Hzx2(ii,j) + cHx(ii) *(Ey(i,j) - Ey(i+1,j))/dxyz;
	         Hz(i,j) = Hz(i,j) + DB*psi_Hzx2(ii,j);
             ii = ii-1;
	     end
      end
      for i = 1:ieT
%.....................................................................
%  CPML for bottom Hz, y-direction
%.....................................................................
         for j = 1:jBC-1
            psi_Hzy1(i,j) = bHy(j)*psi_Hzy1(i,j) + cHy(j)*(Ex(i,j+1) - Ex(i,j))/dxyz;
	        Hz(i,j) = Hz(i,j) + DB*psi_Hzy1(i,j);
         end
%.....................................................................
%  CPML for top Hz, y-direction
%.....................................................................
         jj = jBC-1;
         for j = jhT+1-jBC:jeT
            psi_Hzy2(i,jj) = bHy(jj)*psi_Hzy2(i,jj) + cHy(jj)*(Ex(i,j+1) - Ex(i,j))/dxyz;
	        Hz(i,j) = Hz(i,j) + DB*psi_Hzy2(i,jj);
            jj = jj-1;
         end
      end
      
%----------------------------------------------- 
%-----------------------------------------------  
%-----------------------------------------------  
      n=nTimeStep;
      EzPMLpower(n)=0;
      %{  
      for ja=jBC+1:jeT-jBC
          for ia=iBC+1:2*iBC
              EzPMLpower(n)= EzPMLpower(n) +  abs(Ez(ia,ja)) + abs(Ez(ie-iBC+ia,ja)) ;       %right,left power
          end  
      end
      %}  
% キА