tmz_with_upml_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
exponent = 4;               % exponent of the polynomial grading

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 = 0.7*(exponent+1)/(120*pi*delta);              % optimal value for sigma 
kappa_max = 1;             % maximum value for kappa

%% allocation of memory
dz = zeros(xdim,ydim);      % dz field, electric flux
ez = zeros(xdim,ydim);      % ez field, electric strength
bx = zeros(xdim,ydim);      % bx field, magnetic flux
hx = zeros(xdim,ydim);      % hx field, magnetic strength
by = zeros(xdim,ydim);      % by field, magnetic flux
hy = zeros(xdim,ydim);      % hy field, magnetic strength

%% 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=10*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 of the electric grid in x-direction
xhlength = zeros(1,xdim);   % absolute length of the magnetic grid in x-direction
xekappa = ones(1,xdim);    % grading of kappa of the electric grid in x-direction
xhkappa = ones(1,xdim);    % grading of kappa of the magnetic grid in x-direction
xesigma = zeros(1,xdim);    % grading of sigma of the electric grid in x-direction
xhsigma = zeros(1,xdim);    % grading of sigma of the magnetic grid in x-direction

yelength = zeros(1,ydim);   % absolute length of the electric grid in y-direction
yhlength = zeros(1,ydim);   % absolute length of the magnetic grid in y-direction
yekappa = ones(1,ydim);    % grading of kappa of the electric grid in y-direction
yhkappa = ones(1,ydim);    % grading of kappa of the magnetic grid in y-direction
yesigma = zeros(1,ydim);    % grading of sigma of the electric grid in y-direction
yhsigma = zeros(1,ydim);    % grading of sigma of 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))^exponent;
            xekappa(1,i)=1+(kappa_max - 1)*(abs(xelength(1,i)-xmin)/(thickness*delta))^exponent;
        end
        for i=xdim-thickness:xdim
            xesigma(1,i)=sigma_max*(abs(xelength(1,i)-xmax)/(thickness*delta))^exponent;
            xekappa(1,i)=1+(kappa_max - 1)*(abs(xelength(1,i)-xmax)/(thickness*delta))^exponent;
        end

        for i=1:thickness-1
            xhsigma(1,i)=sigma_max*(abs(xhlength(1,i)-xmin)/(thickness*delta))^exponent;
            xhkappa(1,i)=1+(kappa_max - 1)*(abs(xhlength(1,i)-xmin)/(thickness*delta))^exponent;
        end
        for i=xdim-thickness+1:xdim
            xhsigma(1,i)=sigma_max*(abs(xhlength(1,i)-xmax)/(thickness*delta))^exponent;
            xhkappa(1,i)=1+(kappa_max - 1)*(abs(xhlength(1,i)-xmax)/(thickness*delta))^exponent;
        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))^exponent;
            yekappa(1,i)=1+(kappa_max - 1)*(abs(yelength(1,i)-ymin)/(thickness*delta))^exponent;
        end
        for i=ydim-thickness:ydim
            yesigma(1,i)=sigma_max*(abs(yelength(1,i)-ymax)/(thickness*delta))^exponent;
            yekappa(1,i)=1+(kappa_max - 1)*(abs(yelength(1,i)-ymax)/(thickness*delta))^exponent;
        end

        for i=1:thickness-1
            yhsigma(1,i)=sigma_max*(abs(yhlength(1,i)-ymin)/(thickness*delta))^exponent;
            yhkappa(1,i)=1+(kappa_max - 1)*(abs(yhlength(1,i)-ymin)/(thickness*delta))^exponent;
        end
        for i=ydim-thickness+1:ydim
            yhsigma(1,i)=sigma_max*(abs(yhlength(1,i)-ymax)/(thickness*delta))^exponent;
            yhkappa(1,i)=1+(kappa_max - 1)*(abs(yhlength(1,i)-ymax)/(thickness*delta))^exponent;
        end
        
end


%% calculation of the material operators

DeltaT = delta/(c0*sqrt(2));

dz_decay = ones(1,xdim);
dz_amp = ones(1,xdim);
ez_decay = ones(1,ydim);
ez_amp = ones(1,ydim);

bx_decay = ones(1,ydim);
bx_amp = ones(1,ydim);
hx_amp_forward = ones(1,xdim);
hx_amp_backward = ones(1,xdim);

by_decay = ones(1,xdim);
by_amp = ones(1,xdim);
hy_amp_forward = ones(1,ydim);
hy_amp_backward = ones(1,ydim);

opE = (deltaT/(eps0*delta));
opH = (deltaT/(mue0*delta));

for i=1:xdim
    
    dz_decay(1,i) = (2*eps0*xekappa(1,i) - DeltaT*xesigma(1,i))/(2*eps0*xekappa(1,i) + DeltaT*xesigma(1,i));
    dz_amp(1,i) = idelta*(2*eps0*DeltaT)/(2*eps0*xekappa(1,i) + DeltaT*xesigma(1,i));

    hx_amp_forward(1,i) = (2*eps0*xekappa(1,i) + DeltaT*xesigma(1,i))/(2*mue0*eps0);
    hx_amp_backward(1,i) = (2*eps0*xekappa(1,i) - DeltaT*xesigma(1,i))/(2*mue0*eps0);

    by_decay(1,i) = (2*eps0*xhkappa(1,i) - DeltaT*xhsigma(1,i))/(2*eps0*xhkappa(1,i) + DeltaT*xhsigma(1,i));
    by_amp(1,i) = idelta*(2*eps0*DeltaT)/(2*eps0*xhkappa(1,i) + DeltaT*xhsigma(1,i));
    
end

for j=1:ydim
    
    ez_decay(1,j) = (2*eps0*yekappa(1,j) - DeltaT*yesigma(1,j))/(2*eps0*yekappa(1,j) + DeltaT*yesigma(1,j));
    ez_amp(1,j) = (2 / (2*eps0*yekappa(1,j) + DeltaT*yesigma(1,j)));

    hy_amp_forward(1,j) = (2*eps0*yekappa(1,j) + DeltaT*yesigma(1,j))/(2*mue0*eps0);
    hy_amp_backward(1,j)  = (2*eps0*yekappa(1,j) - DeltaT*yesigma(1,j))/(2*mue0*eps0);

    bx_decay(1,j) = (2*eps0*yhkappa(1,j) - DeltaT*yhsigma(1,j))/(2*eps0*yhkappa(1,j) + DeltaT*yhsigma(1,j));
    bx_amp(1,j) = idelta*(2*eps0*DeltaT)/(2*eps0*yhkappa(1,j) + DeltaT*yhsigma(1,j));

end


clear xmin xmax ymin ymax;
clear xelength xhlength xekappa xhkappa xesigma xhsigma;
clear yelength yhlength yekappa yhkappa yesigma yhsigma;

%% 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

dz_decay(1,1)=0;
dz_decay(1,end)=0;
dz_amp(1,1)=0;
dz_amp(1,end)=0;
ez_decay(1,1)=0;
ez_decay(1,end)=0;
ez_amp(1,1)=0;
ez_amp(1,end)=0;

by_amp(1,end)=0;
by_decay(1,end)=0;
bx_amp(1,end)=0;
bx_decay(1,end)=0;
hx_amp_forward(1,end)=0;
hx_amp_backward(1,end)=0;
hy_amp_forward(1,end)=0;
hy_amp_backward(1,end)=0;

tmp = 0;
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
            
            if ((i>thickness) && (j>thickness) && (i<xdim-thickness-1) && (j<ydim-thickness-1))
                ez(i,j) = ez(i,j) + opE*( hy(i,j) - hy(i-1,j) - hx(i,j) + hx(i,j-1) );
            else
                tmp = dz(i,j);
                dz(i,j) = dz_decay(1,i)*dz(i,j) + dz_amp(1,i)*( hy(i,j) - hy(i-1,j) - hx(i,j) + hx(i,j-1) );
                ez(i,j) = ez_decay(1,j)*ez(i,j) + ez_amp(1,j)*( dz(i,j) - tmp );
           end
        end
    end

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

    
    for j=1:ydim-1
        for i=2:xdim-1
            if ((i>thickness) && (j>thickness) && (i<xdim-thickness-1) && (j<ydim-thickness-1))
                hx(i,j) = hx(i,j) - opH*( ez(i,j+1) - ez(i,j) );
            else
                tmp = bx(i,j);
                bx(i,j) = bx_decay(1,j)*bx(i,j) - bx_amp(1,j)*( ez(i,j+1) - ez(i,j) );
                hx(i,j) = hx(i,j) + hx_amp_forward(1,i)*bx(i,j) - hx_amp_backward(1,i)*tmp;
            end
        end
    end

    
    for j=2:ydim-1
        for i=1:xdim-1
            if ((i>thickness) && (j>thickness) && (i<xdim-thickness-1) && (j<ydim-thickness-1))
                hy(i,j) = hy(i,j) + opH*( ez(i+1,j) - ez(i,j) );
            else
                tmp = by(i,j);
                by(i,j) = by_decay(1,i)*by(i,j) + by_amp(1,i)*( ez(i+1,j) - ez(i,j) );
                hy(i,j) = hy(i,j) + hy_amp_forward(1,j)*by(i,j) - hy_amp_backward(1,j)*tmp;
            end
        end
    end

    
    % visualisation
%     if (mod(n,5)==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 Unaxial PML, step = %i', n));
%         title(sprintf('FDTD Time (UPML) = %.3f nsec',n*deltaT*1e9),'Color',[1 0 0],'FontSize', 22); drawnow;
%     end
    
end

toc

clear dz_decay dz_amp ez_decay ez_amp;
clear bx_decay bx_amp hx_amp_forward hx_amp_backward;
clear by_decay by_amp hy_amp_forward hy_amp_backward;
clear opE opH;
clear dz ez bx hx by hy;
clear idelta spread timeshift endofpulse tsteps;
clear DeltaT c0 eps0 i j mue0 n tmp xdim xexc ydim yexc;
clear centerfrequency delta deltaT exponent frequency kappa_max nsteps sigma_max t0 thickness tw time;