tmz_with_cpml_scan.m

FDTD TMz_with_different_PMLs

%% TMz simulation for the FDTD method

% close all;
clear all;

c0 = 2.99792458e8;
eps0 = 8.854187818e-12;
mue0 =  4*pi*1e-7;

%% user definable parameters
thickness = 10;              % thickness of the pml absorber
exp_sigma = 4;               % exponent of the polynomial grading
exp_alpha = 1;

delta = 0.001;              % distance from discretisation line to line
time = 2e-9;                    % scanned time
deltaT = delta/(c0*sqrt(2));    % timestep
nsteps = ceil(time/deltaT)+1;

xdim = 40+2*thickness;                  % dimension in x direction
ydim = 40+2*thickness;                  % dimension in y direction
centerfrequency = 12e9;
frequency = 12e9;          % frequency bandwidth for the gauss impulse
xexc = ceil((xdim)/2);        % setting the x-coordinate for the excitation
yexc = ceil((ydim)/2);        % setting the y-coordinate for the excitation

% sigma_max = 10;              % maximum value for sigma 
% sigma_max = (exp_sigma+1)/(120*pi*delta);              % optimal value for sigma 
sigma_max = 0.7*((exp_sigma+1)/(120*pi*delta));
kappa_max = 1;             % maximum value for kappa
alpha_max = 0.0;

%% allocation of memory
ez = zeros(xdim,ydim);      % ez field, electric strength
ez_x_psi = zeros(xdim,ydim);
ez_y_psi = zeros(xdim,ydim);
hx = zeros(xdim,ydim);      % hx field, magnetic strength
hx_y_psi = zeros(xdim,ydim);
hy = zeros(xdim,ydim);      % hy field, magnetic strength
hy_x_psi = zeros(xdim,ydim);

%% calculation of important parameters

idelta = 1./delta;
deltaT = (delta/(c0*sqrt(2)));

spread=((pi*frequency)^2)/log(10.);
timeshift=sqrt((5*log(10.))/spread);
endofpulse=20*timeshift;
tsteps = ceil(endofpulse/deltaT);

%% calculation of the graded kappa and sigma
xmin = thickness*delta;                 % first position of the pml absorber in x-direction
xmax = (xdim-thickness-1)*delta;          % last position of the pml absorber in x-direction
ymin = thickness*delta;                 % first position of the pml absorber in y-direction
ymax = (ydim-thickness-1)*delta;          % last position of the pml absorber in y-direction

xelength = zeros(1,xdim);   % absolute length on the electric grid in x-direction
xhlength = zeros(1,xdim);   % absolute length on the magnetic grid in x-direction
xekappa = ones(1,xdim);    % grading of kappa on the electric grid in x-direction
xhkappa = ones(1,xdim);    % grading of kappa on the magnetic grid in x-direction
ixekappa = ones(1,xdim);    % grading of kappa on the electric grid in x-direction
ixhkappa = ones(1,xdim);    % grading of kappa on the magnetic grid in x-direction
xesigma = zeros(1,xdim);    % grading of sigma on the electric grid in x-direction
xhsigma = zeros(1,xdim);    % grading of sigma of the magnetic grid in x-direction
xealpha = zeros(1,xdim);    % grading of alpha on the electric grid in x-direction
xhalpha = zeros(1,xdim);    % grading of alpha on the magnetic grid in x-direction

yelength = zeros(1,ydim);   % absolute length on the electric grid in y-direction
yhlength = zeros(1,ydim);   % absolute length on the magnetic grid in y-direction
yekappa = ones(1,ydim);    % grading of kappa on the electric grid in y-direction
yhkappa = ones(1,ydim);    % grading of kappa on the magnetic grid in y-direction
iyekappa = ones(1,ydim);    % grading of kappa on the electric grid in y-direction
iyhkappa = ones(1,ydim);    % grading of kappa on the magnetic grid in y-direction
yesigma = zeros(1,ydim);    % grading of sigma on the electric grid in y-direction
yhsigma = zeros(1,ydim);    % grading of sigma on the magnetic grid in y-direction
yealpha = zeros(1,ydim);    % grading of alpha on the electric grid in y-direction
yhalpha = zeros(1,ydim);    % grading of alpha on the magnetic grid in y-direction

if thickness~=0

        for i=1:(xdim)
            xelength(1,i) = (i-1)*delta;
            xhlength(1,i) = (i-0.5)*delta;
        end

        for i=1:thickness
            xesigma(1,i)=sigma_max*(abs(xelength(1,i)-xmin)/((thickness)*delta))^exp_sigma;
            xealpha(1,i)=alpha_max*(abs(xelength(1,i+1))/(thickness*delta))^exp_alpha;
            xekappa(1,i)=1+(kappa_max - 1)*(abs(xelength(1,i)-xmin)/(thickness*delta))^exp_sigma;
        end
        for i=xdim-thickness:xdim
            xesigma(1,i)=sigma_max*(abs(xelength(1,i)-xmax)/((thickness)*delta))^exp_sigma;
            xealpha(1,i)=alpha_max*(abs(xelength(1,i)-(xdim-1)*delta)/(thickness*delta))^exp_alpha;
            xekappa(1,i)=1+(kappa_max - 1)*(abs(xelength(1,i)-xmax)/(thickness*delta))^exp_sigma;
        end

        for i=1:thickness-1
            xhsigma(1,i)=sigma_max*(abs(xhlength(1,i)-xmin)/((thickness)*delta))^exp_sigma;
            xhalpha(1,i)=alpha_max*(abs(xhlength(1,i+1))/(thickness*delta))^exp_alpha;
            xhkappa(1,i)=1+(kappa_max - 1)*(abs(xhlength(1,i)-xmin)/(thickness*delta))^exp_sigma;
        end
        for i=xdim-thickness+1:xdim
            xhsigma(1,i)=sigma_max*(abs(xhlength(1,i)-xmax)/((thickness)*delta))^exp_sigma;
            xhalpha(1,i-1)=alpha_max*(abs(xhlength(1,i)-xdim*delta)/(thickness*delta))^exp_alpha;
            xhkappa(1,i)=1+(kappa_max - 1)*(abs(xhlength(1,i)-xmax)/(thickness*delta))^exp_sigma;
        end
        for i=1:xdim
            ixekappa(1,i)=1/xekappa(1,i);
            ixhkappa(1,i)=1/xhkappa(1,i);
        end



        for j=1:ydim
            yelength(1,j) = (j-1)*delta;
            yhlength(1,j) = (j-0.5)*delta;
        end

        for i=1:thickness
            yesigma(1,i)=sigma_max*(abs(yelength(1,i)-ymin)/((thickness)*delta))^exp_sigma;
            yealpha(1,i)=alpha_max*(abs(yelength(1,i+1))/(thickness*delta))^exp_alpha;
            yekappa(1,i)=1+(kappa_max - 1)*(abs(yelength(1,i)-ymin)/(thickness*delta))^exp_sigma;
        end
        for i=ydim-thickness:ydim
            yesigma(1,i)=sigma_max*(abs(yelength(1,i)-ymax)/((thickness)*delta))^exp_sigma;
            yealpha(1,i)=alpha_max*(abs(yelength(1,i)-(ydim-1)*delta)/(thickness*delta))^exp_alpha;
            yekappa(1,i)=1+(kappa_max - 1)*(abs(yelength(1,i)-ymax)/(thickness*delta))^exp_sigma;
        end

        for i=1:thickness-1
            yhsigma(1,i)=sigma_max*(abs(yhlength(1,i)-ymin)/((thickness)*delta))^exp_sigma;
            yhalpha(1,i)=alpha_max*(abs(yhlength(1,i+1))/(thickness*delta))^exp_alpha;
            yhkappa(1,i)=1+(kappa_max - 1)*(abs(yhlength(1,i)-ymin)/(thickness*delta))^exp_sigma;
        end
        for i=ydim-thickness+1:ydim
            yhsigma(1,i)=sigma_max*(abs(yhlength(1,i)-ymax)/((thickness)*delta))^exp_sigma;
            yhalpha(1,i-1)=alpha_max*(abs(yhlength(1,i)-ydim*delta)/(thickness*delta))^exp_alpha;
            yhkappa(1,i)=1+(kappa_max - 1)*(abs(yhlength(1,i)-ymax)/(thickness*delta))^exp_sigma;
        end
        for j=1:ydim
            iyekappa(1,j)=1/yekappa(1,j);
            iyhkappa(1,j)=1/yhkappa(1,j);
        end
        
end


%% calculation of the material operators

DeltaT = delta/(c0*sqrt(2));
ic0 = 1/c0;
iZ0 = sqrt(eps0/mue0);
Z0 = sqrt(mue0/eps0);
Z0p2 = mue0/eps0;

opE = (DeltaT)/(eps0);
ez_x_psi_decay = zeros(1,xdim);
ez_x_psi_amp = zeros(1,xdim);
ez_y_psi_decay = zeros(1,ydim);
ez_y_psi_amp = zeros(1,ydim);

opH = (DeltaT)/(mue0);
hx_decay = ones(1,xdim);
hx_amp = ones(1,xdim);
hx_y_psi_decay = zeros(1,ydim);
hx_y_psi_amp = zeros(1,ydim);

hy_decay = ones(1,ydim);
hy_amp = ones(1,ydim);
hy_x_psi_decay = zeros(1,xdim);
hy_x_psi_amp = zeros(1,xdim);

for i=1:xdim
    
    if i <= thickness
        ez_x_psi_decay(1,i) = exp(-(DeltaT/eps0)*((xesigma(1,i)/xekappa(1,i))+xealpha(1,i)));
        ez_x_psi_amp(1,i) = (xesigma(1,i)*( ez_x_psi_decay(1,i) - 1 ))/(xekappa(1,i) * ( xesigma(1,i) + xealpha(1,i)*xekappa(1,i)));
    end

    if i >= xdim-thickness+1
        ez_x_psi_decay(1,i) = exp(-(DeltaT/eps0)*((xesigma(1,i)/xekappa(1,i))+xealpha(1,i)));
        ez_x_psi_amp(1,i) = (xesigma(1,i)*( ez_x_psi_decay(1,i) - 1 ))/(xekappa(1,i) * ( xesigma(1,i) + xealpha(1,i)*xekappa(1,i)));
    end
    
%     hx_decay(1,i) = (2*mue0*xekappa(1,i) - Z0p2*xesigma(1,i)*DeltaT)/(2*mue0*xekappa(1,i) + Z0p2*xesigma(1,i)*DeltaT);
%     hx_amp(1,i) = (2*Z0*DeltaT)/(2*mue0*xekappa(1,i) + Z0p2*xesigma(1,i)*DeltaT);
%     hx_decay(1,i) = 1;
%     hx_amp(1,i) = (Z0*DeltaT)/(mue0);
    
    if i < thickness
        hy_x_psi_decay(1,i) = exp(-(DeltaT/eps0)*( xhsigma(1,i)./xhkappa(1,i) + xhalpha(1,i) ));
        hy_x_psi_amp(1,i) = (xhsigma(1,i)*( hy_x_psi_decay(1,i) - 1 ))/(xhkappa(1,i) * ( xhsigma(1,i) + xhalpha(1,i)*xhkappa(1,i)));
    end

    if (i >= xdim-thickness+1) && (i~=xdim)
        hy_x_psi_decay(1,i) = exp(-(DeltaT/eps0)*( xhsigma(1,i)./xhkappa(1,i) + xhalpha(1,i) ));
        hy_x_psi_amp(1,i) = (xhsigma(1,i)*( hy_x_psi_decay(1,i) - 1 ))/(xhkappa(1,i) * ( xhsigma(1,i) + xhalpha(1,i)*xhkappa(1,i)));
    end
    
end

for j=1:ydim
    
    if j <= thickness
        ez_y_psi_decay(1,j) = exp(-(DeltaT/eps0)*(yesigma(1,j)/yekappa(1,j)+yealpha(1,j)));
        ez_y_psi_amp(1,j) = (yesigma(1,j)*( ez_y_psi_decay(1,j) - 1 ))/(yekappa(1,j)*( yesigma(1,j) + yealpha(1,j)*yekappa(1,j)));
    end
    
    if j >= ydim-thickness+1
        ez_y_psi_decay(1,j) = exp(-(DeltaT/eps0)*(yesigma(1,j)/yekappa(1,j)+yealpha(1,j)));
        ez_y_psi_amp(1,j) = (yesigma(1,j)*( ez_y_psi_decay(1,j) - 1 ))/(yekappa(1,j)*( yesigma(1,j) + yealpha(1,j)*yekappa(1,j)));
    end
    
%     hy_decay(1,j) = (2*mue0*yekappa(1,j) - Z0p2*yesigma(1,j)*DeltaT)/(2*mue0*yekappa(1,j) + Z0p2*yesigma(1,j)*DeltaT);
%     hy_amp(1,j) = (2*Z0*DeltaT)/(2*mue0*yekappa(1,j) + Z0p2*yesigma(1,j)*DeltaT);
%     hy_decay(1,j) = 1;
%     hy_amp(1,j) = (Z0*DeltaT)/(mue0);
    
    if j < thickness
        hx_y_psi_decay(1,j) = exp(-(DeltaT/eps0)*( yhsigma(1,j)./yhkappa(1,j) + yhalpha(1,j)));
        hx_y_psi_amp(1,j) = (yhsigma(1,j)*( hx_y_psi_decay(1,j) - 1 ))/(yhkappa(1,j)*(yhsigma(1,j) + yhalpha(1,j)*yhkappa(1,j)));
    end
    
    if (j >= xdim - thickness + 1)  && (j~=ydim)
        hx_y_psi_decay(1,j) = exp(-(DeltaT/eps0)*( yhsigma(1,j)./yhkappa(1,j) + yhalpha(1,j)));
        hx_y_psi_amp(1,j) = (yhsigma(1,j)*( hx_y_psi_decay(1,j) - 1 ))/(yhkappa(1,j)*(yhsigma(1,j) + yhalpha(1,j)*yhkappa(1,j)));
    end

end


clear xelength xhlength xekappa xhkappa xesigma xhsigma xealpha xhalpha;
clear yelength yhlength yekappa yhkappa yesigma yhsigma yealpha yhalpha;

%% allocate array for examination

stime = zeros(1,nsteps);
excitation = zeros(1,nsteps);
uxdirect = zeros(1,nsteps);
uydirect = zeros(1,nsteps);
uxycorner = zeros(1,nsteps);


%% calculation of the time loop

tmp_dx = 0;
tmp_dy = 0;

opExc = opE*idelta;
tw = 26.53e-12;
t0 = 4*tw;

tic
for n=1:nsteps

    stime(1,n) = (n-1)*deltaT;
    excitation(1,n) = delta*ez(xexc,yexc);
    uxdirect(1,n) = delta*ez(xexc-18,yexc);
    uydirect(1,n) = delta*ez(xexc,yexc-18);
    uxycorner(1,n) = delta*ez(xexc-18,yexc-18);

    for j=2:ydim-1
        for i=2:xdim-1
            
            tmp_dx = idelta*(hy(i,j) - hy(i-1,j));
            tmp_dy = idelta*(hx(i,j) - hx(i,j-1));
            
            ez_x_psi(i,j) = ez_x_psi_decay(1,i) * ez_x_psi(i,j) + ez_x_psi_amp(1,i) * tmp_dx;
            ez_y_psi(i,j) = ez_y_psi_decay(1,j) * ez_y_psi(i,j) + ez_y_psi_amp(1,j) * tmp_dy;

%             ez(i,j) = ez(i,j) + opE*( tmp_dx - tmp_dy + ez_x_psi(i,j) - ez_y_psi(i,j));
            ez(i,j) = ez(i,j) + opE*( ixekappa(1,i)*tmp_dx - iyekappa(1,j)*tmp_dy + ez_x_psi(i,j) - ez_y_psi(i,j));
%             ez(i,j) = ez(i,j) + opE*( tmp_dx - tmp_dy);
%             ez(i,j) = ez(i,j) + opE * ez_x_psi(i,j);
%             ez(i,j) = ez(i,j) - opE * ez_y_psi(i,j);
        end
    end

    
    % soft excitation
    % ez(xexc,yexc) = ez(xexc,yexc) - opExc*2.0*(((n-1)*deltaT-t0)/tw)*exp(-(((n-1)*deltaT-t0)/tw)*(((n-1)*deltaT-t0)/tw));
    ez(xexc,yexc) = ez(xexc,yexc) - opExc*cos(2*pi*centerfrequency*(n-0.5)*deltaT)*exp(-((((n-0.5)*deltaT-timeshift)^2)*spread));
    % ez(xexc,yexc) = ez(xexc,yexc) - opExc*exp(-((((n-1)*deltaT-timeshift)^2)*spread));
    % ez(xexc,yexc) = ez(xexc,yexc) - opExc*sin(2*pi*frequency*(n-1)*deltaT);

    
    for j=1:ydim-1
        for i=2:xdim-1
            
            tmp_dy = idelta*(ez(i,j+1) - ez(i,j));
            
            hx_y_psi(i,j) = hx_y_psi_decay(1,j) * hx_y_psi(i,j) + hx_y_psi_amp(1,j) * tmp_dy;
            
%             hx(i,j) = hx(i,j) - opH * ( tmp_dy + hx_y_psi(i,j) );
            hx(i,j) = hx(i,j) - opH * ( iyhkappa(1,j) * tmp_dy + hx_y_psi(i,j) );
%             hx(i,j) = hx_decay(1,i)*hx(i,j) - hx_amp(1,i) * ( iyhkappa(1,j) * tmp_dy + hx_y_psi(i,j) );
%             hx(i,j) = hx_decay(1,i)*hx(i,j) - hx_amp(1,i) * iyhkappa(1,j) * tmp_dy;
%             hx(i,j) = hx(i,j) - hx_amp(1,i) * hx_y_psi(i,j);
            
        end
    end

    
    for j=2:ydim-1
        for i=1:xdim-1
            
            tmp_dx = idelta*(ez(i+1,j) - ez(i,j));
            
            hy_x_psi(i,j) = hy_x_psi_decay(1,i) * hy_x_psi(i,j) + hy_x_psi_amp(1,i) * tmp_dx;
            
%             hy(i,j) = hy(i,j) + opH * ( tmp_dx + hy_x_psi(i,j) );
            hy(i,j) = hy(i,j) + opH * ( ixhkappa(1,i) * tmp_dx + hy_x_psi(i,j) );
%             hy(i,j) = hy_decay(1,j)*hy(i,j) + hy_amp(1,j) * ( ixhkappa(1,i) * tmp_dx + hy_x_psi(i,j) );
%             hy(i,j) = hy_decay(1,j)*hy(i,j) + hy_amp(1,j) * ixhkappa(1,i) * tmp_dx;
%             hy(i,j) = hy(i,j) + hy_amp(1,j) * hy_x_psi(i,j);
            
        end
    end

    
    % visualisation
%     if (mod(n,2)==0)
%         surf(ez); axis on; zlim([-5 +5]); shading interp; colormap jet; caxis([-5 +5]); lighting phong;  
%         set(gcf,'Color', 'white', 'Number', 'off', 'Name', sprintf('Tiny FDTD with Convolutional PML, step = %i', n));
%         title(sprintf('FDTD Time (CPML) = %.3f nsec',n*deltaT*1e9),'Color',[1 0 0],'FontSize', 22); drawnow;
%     end

end

toc

clear ic0 iZ0 Z0;
clear ez ez_x_psi ez_y_psi hx hx_y_psi hy hy_x_psi;
clear Z0p2 opE ez_x_psi_decay ez_x_psi_amp ez_y_psi_decay ez_y_psi_amp;
clear opH hx_decay hx_amp hx_y_psi_decay hx_y_psi_amp hy_decay hy_amp hy_x_psi_decay hy_x_psi_amp;
clear xmin xmax ymin ymax ixekappa ixhkappa iyekappa iyhkappa;
clear i j endofpulse eps0 mue0 c0 idelta n opExc spread timeshift tmp_dx tmp_dy xdim xexc ydim yexc DeltaT time;
clear alpha_max centerfrequency delta deltaT exp_alpha exp_sigma frequency kappa_max nsteps sigma_max t0 thickness tw tsteps;