%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Photonic_crystal_resonant_cavities
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
clear;
clear all;
clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Total number of time steps
MaxTime  = 250;
%--------------------------------    
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;         
%--------------------------------  
dxyz = .25E-3;
%dxyz = lambda/8;
dt =   dxyz/(c0);
%--------------------------------  
% Specify the Impulsive Source
 %  tw = 53.0E-12; tO = 4.0*tw;
   tw=6; t0 = 20;
% PML thickness in each direction
kBC = 40; 
%--------------------------------  
% Size of main grid 
ke=160; 
kh=ke+1;   
% Size of total computational domain      
keT=ke+2*kBC;        
khT=keT+1;
%--------------------------------  
Ex =zeros(khT,1);  
Hy =zeros(khT,1); 
CA =zeros(khT,1);
CB =zeros(khT,1);
sig =zeros(khT,1);
epsi=zeros(khT,1);
%--------------------------------  
ksrc=round(keT/2); 
%--------------------------------  
% Specify the CPML Order and Other Parameters 
m = 3; ma = 1 ;
sigZmax =  (0.8*(m+1)/(dxyz*(mu0/eps0*epsR)^0.5));
aZmax = 0.05;
kZmax = 1.0;
%--------------------------------      
record_grid = MaxTime;
NoFig=99;
Nplt_jmp=5;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
psi_Exz1=zeros(khT,1);psi_Exz2=zeros(khT,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
psi_Hyz1=zeros(khT,1);psi_Hyz2=zeros(khT,1);       
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
bEz=zeros(kBC,1);cEz=zeros(kBC,1) ;
A_EzBC=zeros(kBC,1) ;
sigEzBC=zeros(kBC,1) ;
K_EzBC=zeros(kBC,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 
F_ez=zeros(keT,1);
F_hz=zeros(keT,1);
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  INITIALIZE VARIABLES
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
sig(:) = sigM1;
epsi(:) = epsR*eps0;  
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  SET CPML PARAMETERS IN EACH DIRECTION
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  z-dir
%
   for k = 1:kBC
      sigEzBC(k) = sigZmax * ( (kBC - k ) / (kBC - 1.0) )^m;
      K_EzBC(k) = 1.0+(kZmax-1.0)*((kBC - k) / (kBC - 1.0))^m;
      A_EzBC(k) = aZmax*((k-1)/(kBC-1.0))^ma;
      bEz(k) = exp(-(sigEzBC(k) / K_EzBC(k) + A_EzBC(k))*dt/eps0);
      if ((sigEzBC(k) == 0.0) &&    ...
         (A_EzBC(k) == 0.0) && (k == kBC)) 
         cEz(k) = 0.0;
      else
         cEz(k) = sigEzBC(k)*(bEz(k)-1.0)/ (sigEzBC(k)+K_EzBC(k)*A_EzBC(k)) / K_EzBC(k);
      end
   end
   for k = 1:kBC-1
      sigHz_BC(k) = sigZmax * ( (kBC - k - 0.5)/(kBC-1.0))^m;
      K_HzBC(k) = 1.0+(kZmax-1.0)*((kBC - k - 0.5) / (kBC - 1.0))^m;
      A_HzBC(k) = aZmax*((k-0.5)/(kBC-1.0))^ma;
      bHz(k) = exp(-(sigHz_BC(k) / K_HzBC(k) + A_HzBC(k))*dt/eps0);
      cHz(k) = sigHz_BC(k)*(bHz(k)-1.0)/ (sigHz_BC(k)+K_HzBC(k)*A_HzBC(k)) / K_HzBC(k);
   end
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  FILL IN UPDATING COEFFICIENTS
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   DA = 1.0;
   DB = dt/mu0; 
%-------------------------
for k = 1:keT 
    CA(k) = (1.0 - sig(k)*dt / (2.0*epsi(k))) / (1.0 + sig(k) * dt / (2.0*epsi(k)));
    CB(k) = (dt/epsi(k)) / (1.0 + sig(k)*dt / (2.0*epsi(k))); 
end
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  FILL IN DENOMINATORS FOR FIELD UPDATES
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  CPML_Psi
   kk = kBC-1;
   for k = 1:keT
      if k <= kBC-1
         F_hz(k) = 1.0/(K_HzBC(k)*dxyz);
      elseif k >= khT+1-kBC
         F_hz(k) = 1.0/(K_HzBC(kk)*dxyz);
         kk = kk - 1;
      else
         F_hz(k) = 1.0/dxyz;
      end
   end 
   kk = kBC;
   for k = 1:keT
      if  k <= kBC  
         F_ez(k) = 1.0/(K_EzBC(k)*dxyz);
      elseif  k >= khT+1-kBC 
         F_ez(k) = 1.0/(K_EzBC(kk)*dxyz);
         kk = kk - 1;
      else
         F_ez(k) = 1.0/dxyz;
      end
   end
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  BEGIN TIME STEP
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
for nTimeStep = 1:MaxTime 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  UPDATE Ex
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
for k = 2:keT 
	Ex(k) = CA(k) * Ex(k) + CB(k) *  (Hy(k-1) - Hy(k))*F_ez(k) ;   
end 
%.....................................................................
%  CPML for bottom Ex, k-direction
%.....................................................................
         for k = 2:kBC
     	     psi_Exz1(k) = bEz(k)*psi_Exz1(k) + cEz(k) *(Hy(k-1) - Hy(k))/dxyz;
             Ex(k) = Ex(k) + CB(k)*psi_Exz1(k);
         end
%.....................................................................
%  CPML for top Ex, k-direction
%.....................................................................
         kk = kBC;
         for k = khT+1-kBC:keT
  	         psi_Exz2(kk) = bEz(kk)*psi_Exz2(kk) + cEz(kk) *(Hy(k-1) - Hy(k))/dxyz;
             Ex(k) = Ex(k) + CB(k)*psi_Exz2(kk);
             kk = kk-1;
         end 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-----------------------------------------------------------------------
%   SOURCE
%-----------------------------------------------------------------------
source =  exp(-.5*((t0-nTimeStep)/tw)^2.0) ;
%source = sin(2*(nTimeStep*dt/t0^2.0));
Ex(ksrc) = Ex(ksrc)+source ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  UPDATE Hy
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
for k = 1:keT 
    Hy(k) = DA * Hy(k) + DB * (Ex(k) - Ex(k+1))*F_hz(k)  ;  
end 
%.....................................................................
%  CPML for bottom Hy, k-direction
%.....................................................................
         for k = 1:kBC-1
	        psi_Hyz1(k) = bHz(k)*psi_Hyz1(k) + cHz(k)*(Ex(k) - Ex(k+1))/dxyz;
	        Hy(k) = Hy(k) + DB*psi_Hyz1(k);
         end
%.....................................................................
%  CPML for top Hy, k-direction
%.....................................................................
         kk = kBC-1;
         for k = khT+1-kBC:keT
	         psi_Hyz2(kk) = bHz(kk)*psi_Hyz2(kk) + cHz(kk)*(Ex(k) - Ex(k+1))/dxyz;
	         Hy(k) = Hy(k) + DB*psi_Hyz2(kk);
             kk = kk-1;
         end 
%--------------------------------------------------------------------------
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%   VISUALIZATION
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  if mod(nTimeStep,Nplt_jmp)==0
      NoFig=NoFig+1 ;
      timestep=int2str(nTimeStep); 
    %----------------------------------------------------
    %----------------------------------------------------
    %
      subplot(3,1,1),plot(Ex);   
      axis([1 keT -1 1]);
      title(['Ex at time step = ',timestep]);  
    %}
    %------------------------------  
    %
      subplot(3,1,2),plot(Hy);   
      axis([1 keT -1.5e-3 1.5e-3]); 
      title(['Hy']); 
    %}
    %----------------------------------------------------
    %----------------------------------------------------
    % saveas(gcf,int2str(NoFig), 'jpg') 
   pause(0.5)
   end
%-----------------------------------------------------------------------
end
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  end TIME STEP
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~