photonic_crystal_resonant_cavities.asv

FDTD CPML 計算方型共振腔電場頻率

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Photonic_crystal_resonant_cavities
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
clear;
clear all;
clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Total number of time steps
MaxTime  = 1024;
  
epsR = 1.0;
sigM1 = 0.0 ; 
  
c0 = 2.99792458E8;
mu0 = 4.0 * pi * 1.0E-7; 
eps0 = 1.0/(c0*c0*mu0);
  
a=0.58652e-6;
freq=1.9341e+14;
lambda=c0/freq;         
omega=2.0*pi*freq;      wa=omega*a/2*pi*c0; 

dxyz = lambda/8;
dt = 0.99 / (c0*(1.0/dxyz^2+1.0/dxyz^2 )^0.5);

% Specify the Impulsive Source
tw = 53.0E-12;      tO = 4.0*tw;
t0=70e-15;         T0=30e-15;

% Size of main grid 
ie=160;          
je=160;
ke=160;
ih=ie+1;
jh=je+1;   
kh=ke+1;
% PML thickness in each direction
iBC = 40;    
jBC = iBC;  
kBC = iBC;
% Size of total computational domain      
ieT=ie+2*iBC;          
jeT=je+2*jBC;        
keT=ke+2*kBC;        
ihT=ieT+1;
jhT=jeT+1;          
khT=keT+1;

Ex =zeros(ieT, jhT);Ey =zeros(ihT, jeT);Ez =zeros(ihT, jhT);
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); 
ta=zeros(MaxTime,1);
isrc=round(ihT/2);
jsrc=round(jhT/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.24;
aYmax = aXmax; 
aZmax = aXmax;
kXmax = 15.0;
kYmax = kXmax; 
kZmax = kXmax;
      
EzPMLpower=zeros(MaxTime,1);
for n=1:MaxTime
ta(n)=n+1;
end

record_grid = 1024 ;
NoFig=99;
Nplt_jmp=20;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  CPML 
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;
epsi(:,:) = epsR*eps0;  
%     Add metal cylinder
diam=39*0.24;          % diameter of cylinder
rad=diam/2.0;          % radius of cylinder

% 膀非翴蛾
icenter01=ieT/2;icenter001=70.2*ieT/234;          
jcenter01=70.2*jeT/234;jcenter001=jeT/2;     
%---------------------------------------
icenter02=ieT/2;icenter002=93.6*ieT/234;
jcenter02=93.6*jeT/234;jcenter002=jeT/2;
%---------------------------------------
icenter03=ieT/2;icenter003=140.4*ieT/234;
jcenter03=140.4*jeT/234;jcenter003=jeT/2;
%---------------------------------------
icenter04=ieT/2;icenter004=163.8*ieT/234;
jcenter04=163.8*jeT/234;jcenter004=jeT/2;

%