photonic_crystal_two_core_power.m

FDTD CPML 同時計算方型and六角形共振腔

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% 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);
  
freq=1.9341e+14;
lambda=c0/freq;         
omega=2.0*pi*freq; 

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=320;        
je=160; 
ke=160;  
ih=ie+1; 
jh=je+1;   
kh=ke+1; 
% PML thickness in each direction
iBC = 40;    
jBC = iBC;  
kBC = iBC;
iBCm = 40;

ie1=160;
ie2=160;
% Size of total computational domain      
ieT=ie+3*iBC;
jeT=je+2*jBC;        
keT=ke+2*kBC;        
ihT=ieT+1;
jhT=jeT+1;          
khT=keT+1;

iBCM=ih+2*iBC;

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);

isrc01=round(3*ihT/11);
jsrc01=round(3*jhT/6);
isrc02=round(8*ihT/11);
jsrc02=round(3*jhT/6);
%  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;

EzPMLpower01=zeros(MaxTime,1);
EzPMLpower02=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);   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
psi_Eyx1M=zeros(iBCM,jeT);psi_Eyx2M=zeros(iBCM,jeT);
psi_Ezx1M=zeros(iBCM,jhT);psi_Ezx2M=zeros(iBCM,jhT);
psi_Hyx1M=zeros(iBCM-1,jhT);psi_Hyx2M=zeros(iBCM-1,jhT);
psi_Hzx1M=zeros(iBCM-1,jeT);psi_Hzx2M=zeros(iBCM-1,jeT);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
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) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
bExM=zeros(iBCM,1);cExM=zeros(iBCM,1); 
A_ExBCM=zeros(iBCM,1); 
sigExBCM=zeros(iBCM,1); 
K_ExBCM=zeros(iBCM,1);  

bHxM=zeros(iBCM-1,1);cHxM=zeros(iBCM-1,1);
A_HxBCM=zeros(iBCM-1,1);
sigHx_BCM=zeros(iBCM-1,1);
K_HxBCM=zeros(iBCM-1,1); 

% Denominators for the update equations 
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;

% Square Cavities----------------------------------------------------------
%
% Add metal cylinder