untitled.m

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

%***********************************************************************
%     2-D FDTD code with CPML absorbing boundary conditions
%***********************************************************************  
%***********************************************************************
% fdtd2D_CPML 
%  ..................................
 clear;
 clear all;
 clc;
 %function [TP,eps0,dt,dxyz,freq,lambda,EzPMLpower]=fdtd2D_cpml_e01_newxx(ax)
 
%  Input Fundamental Constants (MKS units) 
      c0 = 2.99792458E8;
      mu0 = 4.0 * pi * 1.0E-7; 
      eps0 = 1.0/(c0*c0*mu0);

%  ..................................
%  Specify Material Relative Permittivity and Conductivity 
      epsR = 1.0;
      sigM1 = 0.0 ; % free space

%  ..................................
%  Specify Number of Time Steps and Grid Size Parameters 
      MaxTime  = 500;  % total number of time steps
%  ..................................
ie=300;          % Size of main grid
je=300;
ke=300;
ih=ie+1;
jh=je+1;   
kh=ke+1;   
%  ..................................
%  Specify the CPML Thickness in Each Direction (Value of Zero 
%  Corresponds to No PML, and the Grid is Terminated with a PEC) 
      % PML thickness in each direction 
      iBC = 21;    
      jBC = iBC;  
      kBC = iBC; 
%
ieT=ie+2*iBC;          % Size of total computational domain
jeT=je+2*jBC;        
keT=ke+2*kBC;        
ihT=ieT+1;
jhT=jeT+1;          
khT=keT+1;          
                                        % number of Ez field components %      
      Ex =zeros(ieT, jhT) ;
      Ey =zeros(ihT, jeT);
      Ez =zeros(ihT, jhT);
      
      EzPMLpower=zeros(1,ie);
      EzPMLpowert=zeros(1,ie);
%    
      Hx =zeros(ihT, jeT) ;
      Hy =zeros(ieT, jhT);
      Hz =zeros(ieT, jeT); 
      
      CA  =zeros(ihT, jhT);
      CB  =zeros(ihT, jhT);
      sig =zeros(ihT, jhT);
      epsi=zeros(ihT, jhT); 
%  ..................................
%  Specify Grid Cell Size in Each Direction and Calculate the 
%  Resulting Courant-Stable Time Step
      
      lambda=5.5e-7;             %wavelength bandwidth of source excitation
      
      freq=c0/lambda;                %frequency bandwidth of source excitation
      
      omega=2.0*pi*freq; 
      
      dxyz = lambda/20;    % cell size in each direction
      
      dt = 0.99 / (c0*(1.0/dxyz^2+1.0/dxyz^2 )^0.5);
      
      t_periodic=1/(freq*dt);                                           
      
      ha=1/freq;
      
      ta=zeros(MaxTime,1);               % time step increment
      EzPMLpower=zeros(MaxTime,1);
      for n=1:MaxTime
         ta(n)=(n-1)*dt;
      end
  
%  ..................................
%  Specify the Impulsive Source (See Equation 7.134) 
      tw = 53.0E-12; t0 = 4.0*tw  ;
      

rtau=50.0e-12;
tau=rtau/dt;
delay=tau;
srcconst=-dt*3.0e+11;

   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  ..................................
%  Specify the Time Step at which the Grid is Recorded for Visualization 
      record_grid = 300 ;
      NoFig=99;
      Nplt_jmp=10;

%  ..................................
%  Source point.  
   isrc=round(ihT/2);         % Location of     source
   jsrc=round(jhT/2); 
%      
      iview= floor(ieT/2);
      jview= floor(jeT/2); 
%  ..................................
%  Specify the CPML Order and Other Parameters 
      m = 3; ma = 1 ;
      sigXmax =  (0.8*(m+1)/(dxyz*(mu0/eps0*epsR)^0.5));
      sigYmax =  (0.8*(m+1)/(dxyz*(mu0/eps0*epsR)^0.5));
      sigZmax =  (0.8*(m+1)/(dxyz*(mu0/eps0*epsR)^0.5));
      aXmax = 0.0003;
      aYmax = aXmax; 
      aZmax = aXmax;
      kXmax = 15.0;
      kYmax = kXmax; 
      kZmax = kXmax;   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  CPML (7.100)
      psi_Exy1=zeros(ieT,jBC)  ;     
      psi_Exy2=zeros(ieT,jBC);
      
      psi_Exz1=zeros(ieT,jhT); 
      psi_Exz2=zeros(ieT,jhT) ;
 
      psi_Eyx1=zeros(iBC,jeT);
      psi_Eyx2=zeros(iBC,jeT);
      
      psi_Eyz1=zeros(ihT,jeT); 
      psi_Eyz2=zeros(ihT,jeT) ;
      
      psi_Ezx1=zeros(iBC,jhT); 
      psi_Ezx2=zeros(iBC,jhT) ; 
 
      psi_Ezy1=zeros(ihT,jBC);    
      psi_Ezy2=zeros(ihT,jBC);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
      psi_Hxy1=zeros(ihT,jBC-1) ;  
      psi_Hxy2=zeros(ihT,jBC-1) ;
 
      psi_Hxz1=zeros(ihT,jeT); 
      psi_Hxz2=zeros(ihT,jeT) ;
 
      psi_Hyz1=zeros(ieT,jhT) ; 
      psi_Hyz2=zeros(ieT,jhT);
 
      psi_Hyx1=zeros(iBC-1,jhT); 
      psi_Hyx2=zeros(iBC-1,jhT) ; 
 
      psi_Hzx1=zeros(iBC-1,jeT);  
      psi_Hzx2=zeros(iBC-1,jeT);  
      
      psi_Hzy1=zeros(ieT,jBC-1) ;  
      psi_Hzy2=zeros(ieT,jBC-1);         
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
      bEx=zeros(iBC,1); 
      cEx=zeros(iBC,1); 
      A_ExBC=zeros(iBC,1); 
      sigExBC=zeros(iBC,1); 
      K_ExBC=zeros(iBC,1);   
 
      bEy=zeros(jBC,1);
      cEy=zeros(jBC,1);
      A_EyBC=zeros(jBC,1);
      sigEyBC=zeros(jBC,1);
      K_EyBC=zeros(jBC,1); 
      
      bEz=zeros(kBC,1) ;
      cEz=zeros(kBC,1) ;
      A_EzBC=zeros(kBC,1) ;
      sigEzBC=zeros(kBC,1) ;
      K_EzBC=zeros(kBC,1) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
% (7.105)
      bHx=zeros(iBC-1,1);
      cHx=zeros(iBC-1,1);
      A_HxBC=zeros(iBC-1,1);
      sigHx_BC=zeros(iBC-1,1);
      K_HxBC=zeros(iBC-1,1); 
      
      bHy=zeros(jBC-1,1) ;
      cHy=zeros(jBC-1,1) ;
      A_HyBC=zeros(jBC-1,1) ;
      sigHy_BC=zeros(jBC-1,1) ;
      K_HyBC=zeros(jBC-1,1) ; 
      
      bHz=zeros(kBC-1,1) ;
      cHz=zeros(kBC-1,1) ;
      A_HzBC=zeros(kBC-1,1) ;
      sigHz_BC=zeros(kBC-1) ;
      K_HzBC=zeros(kBC-1,1) ;

%     denominators for the update equations 
%   (7.109)
      F_ex=zeros(ieT,1);
      F_hx=zeros(ieT,1);
 
      F_ey=zeros(jeT,1);
      F_hy=zeros(jeT,1);
 
      F_ez=zeros(keT,1);
      F_hz=zeros(keT,1);
 
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  INITIALIZE VARIABLES
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
 
   sig(:,:) = sigM1;
   sig(1:isrc-18,jBC+90:jBC+120)=2e+12;
   sig(isrc+18:ihT,jBC+90:jBC+120)=2e+12;
   sig(1:isrc-18,jBC+160:jBC+190)=2e+12;
   sig(isrc+18:ihT,jBC+160:jBC+190)=2e+12;
   epsi(:,:)=epsR*eps0;
   
   %{
   neff=1.5;
   b_periodic=floor(lambda/(2*neff*dxyz));
   delta_n=neff*(0.0001);
   for ia=1:ihT
       for ja=1:jhT
           if ia>=40 && ia<=80 && ja>=30 && ja<=90
                 epsi(ia,ja)=eps0*(neff+delta_n*(1+cos(2*pi*ia/b_periodic)))^2;
           elseif ia>=26 && ia<40 && ja>=30 && ja<=90
                 epsi(ia,ja)=(neff*(1-0.003))^2*eps0;
           elseif ia>80 && ia<=94 && ja>=30 && ja<=90
                 epsi(ia,ja)=(neff*(1-0.003))^2*eps0;
           else
                 epsi(ia,ja)=epsR*eps0;
           end
       end
   end
   %}
 
%   write(*,*)"ihT: ", ihT
%   write(*,*)"jhT: ", jhT
%   write(*,*)"khT: ", khT
%   write(*,*)"dt: ", dt
%   write(*,*)"MaxTime : ", MaxTime 
%   write(*,*)"max time: ", MaxTime *dt
%   write(*,*)"record grid after ", record_grid, "dt"
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  SET CPML PARAMETERS IN EACH DIRECTION
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  x-dir
%
  for i = 1:iBC
      sigExBC(i) = sigXmax * ( (iBC - i) / (iBC - 1.0) )^m ;
      K_ExBC(i) = 1.0+(kXmax-1.0)*((iBC - i) / (iBC - 1.0))^m ;
      A_ExBC(i) = aXmax*((i-1.0)/(iBC-1.0))^ma ;%(7.79)
      bEx(i) = exp(-(sigExBC(i) / K_ExBC(i) + A_ExBC(i))*dt/eps0) ;
      if ( (sigExBC(i) == 0.0) && (A_ExBC(i) == 0.0) && (i == iBC) )  
         cEx(i) = 0.0;
      else
         cEx(i) = sigExBC(i)*(bEx(i)-1.0)/(sigExBC(i)+K_ExBC(i)*A_ExBC(i)) / K_ExBC(i);
      end
  end
  for i = 1:iBC-1
      sigHx_BC(i) = sigXmax * ( (iBC - i - 0.5)/(iBC-1.0))^m;
      K_HxBC(i) = 1.0+(kXmax-1.0)*((iBC - i - 0.5) / (iBC - 1.0))^m ;
      A_HxBC(i) = aXmax*((i-0.5)/(iBC-1.0))^ma;
      bHx(i) = exp(-(sigHx_BC(i) / K_HxBC(i) + A_HxBC(i))*dt/eps0) ;
      cHx(i) = sigHx_BC(i)*(bHx(i)-1.0)/ (sigHx_BC(i)+K_HxBC(i)*A_HxBC(i)) / K_HxBC(i);
  end 
%  y-dir
%
   for j = 1:jBC
      sigEyBC(j) = sigYmax * ( (jBC - j ) / (jBC - 1.0) )^m;
      K_EyBC(j) = 1.0+(kYmax-1.0)* ((jBC - j) / (jBC - 1.0))^m;
      A_EyBC(j) = aYmax*((j-1)/(jBC-1.0))^ma;
      bEy(j) = exp(-(sigEyBC(j) / K_EyBC(j) + A_EyBC(j))*dt/eps0);
      if ( (sigEyBC(j) == 0.0) && (A_EyBC(j) == 0.0) && (j == jBC)) 
         cEy(j) = 0.0;
      else
         cEy(j) = sigEyBC(j)*(bEy(j)-1.0)/(sigEyBC(j)+K_EyBC(j)*A_EyBC(j)) / K_EyBC(j);
      end
   end
   for j = 1:jBC-1
      sigHy_BC(j) = sigYmax * ( (jBC - j - 0.5)/(jBC-1.0))^m;
      K_HyBC(j) = 1.0+(kYmax-1.0)*((jBC - j - 0.5) / (jBC - 1.0))^m;
      A_HyBC(j) = aYmax*((j-0.5)/(jBC-1.0))^ma;
      bHy(j) = exp(-(sigHy_BC(j) / K_HyBC(j) + A_HyBC(j))*dt/eps0);
      cHy(j) = sigHy_BC(j)*(bHy(j)-1.0)/ (sigHy_BC(j)+K_HyBC(j)*A_HyBC(j)) / K_HyBC(j);
   end 
 
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  FILL IN UPDATING COEFFICIENTS
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   DA = 1.0; %(7-108)
   DB = dt/mu0; 
% (7.107)