tmz_with_spml_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
delta = 0.001;              % distance from discretisation line to line
thickness = 10;              % thickness of the pml absorber
exponent = 4;               % exponent of the polynomial grading
sigma_max = 0.7*(exponent+1)/(120*pi*delta);             % maximum value for kappa

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


%% allocation of memory
ez = zeros(xdim,ydim);      % ez field, electric strength
ezx = zeros(xdim,ydim);      % ezx field, electric strength
ezy = zeros(xdim,ydim);      % ezy field, electric strength
hx = zeros(xdim,ydim);      % hx field, magnetic strength
hy = zeros(xdim,ydim);      % hy field, magnetic strength

%% calculation of important parameters

idelta = 1./delta;

spread=((pi*frequency)^2)/log(10.);
timeshift=sqrt((5*log(10.))/spread);
endofpulse=10.0*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 = zeros(1,xdim);    % grading of kappa of the electric grid in x-direction
xhkappa = zeros(1,xdim);    % grading of kappa 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 = zeros(1,ydim);    % grading of kappa of the electric grid in y-direction
yhkappa = zeros(1,ydim);    % grading of kappa 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
            xekappa(1,i)=sigma_max*(abs(xelength(1,i)-xmin)/(thickness*delta))^exponent;
        end
        for i=xdim-thickness:xdim
            xekappa(1,i)=sigma_max*(abs(xelength(1,i)-xmax)/(thickness*delta))^exponent;
        end

        for i=1:thickness-1
            xhkappa(1,i)=sigma_max*(abs(xhlength(1,i)-xmin)/(thickness*delta))^exponent;
        end
        for i=xdim-thickness+1:xdim
            xhkappa(1,i)=sigma_max*(abs(xhlength(1,i)-xmax)/(thickness*delta))^exponent;
        end



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

        for i=1:thickness
            yekappa(1,i)=sigma_max*(abs(yelength(1,i)-ymin)/(thickness*delta))^exponent;
        end
        for i=ydim-thickness:ydim
            yekappa(1,i)=sigma_max*(abs(yelength(1,i)-ymax)/(thickness*delta))^exponent;
        end

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


%% calculation of the material operators

z0p2 = mue0/eps0;
DeltaT = delta/(c0*sqrt(2));

ezy_decay = ones(1,xdim);
ezy_amp = ones(1,xdim);
ezx_decay = ones(1,ydim);
ezx_amp = ones(1,ydim);

hx_decay = ones(1,ydim);
hx_amp = ones(1,ydim);

hy_decay = ones(1,xdim);
hy_amp = ones(1,xdim);

opE = ((DeltaT)/(delta*eps0));

for i=1:xdim
    
    ezy_decay(1,i) = (2*eps0 - DeltaT*xekappa(1,i))/(2*eps0 + DeltaT*xekappa(1,i));
    ezy_amp(1,i) = idelta*((2*DeltaT)/(2*eps0 + DeltaT*xekappa(1,i)));

    hy_decay(1,i) = (2*eps0 - DeltaT*xhkappa(1,i))/(2*eps0 + DeltaT*xhkappa(1,i));
    hy_amp(1,i) = idelta*((2*DeltaT)/(2*mue0 + DeltaT*z0p2*xhkappa(1,i)));
    
end

for j=1:ydim
    
    ezx_decay(1,j) = (2*eps0 - DeltaT*yekappa(1,j))/(2*eps0 + DeltaT*yekappa(1,j));
    ezx_amp(1,j) = idelta*((2*DeltaT)/(2*eps0 + DeltaT*yekappa(1,j)));

    hx_decay(1,j) = (2*eps0 - DeltaT*yhkappa(1,j))/(2*eps0 + DeltaT*yhkappa(1,j));
    hx_amp(1,j) = idelta*((2*DeltaT)/(2*mue0 + DeltaT*z0p2*yhkappa(1,j)));

end


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

ezy_decay(1,1)=0;
ezy_decay(1,end)=0;
ezy_amp(1,1)=0;
ezy_amp(1,end)=0;
ezx_decay(1,1)=0;
ezx_decay(1,end)=0;
ezx_amp(1,1)=0;
ezx_amp(1,end)=0;

hy_decay(1,end)=0;
hx_decay(1,end)=0;
hy_amp(1,end)=0;
hx_amp(1,end)=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
                ezy(i,j) = ezy_decay(1,i)*ezy(i,j) + ezy_amp(1,i)*( hy(i,j) - hy(i-1,j) );
                ezx(i,j) = ezx_decay(1,j)*ezx(i,j) - ezx_amp(1,j)*( hx(i,j) - hx(i,j-1) );
                ez(i,j) = ezx(i,j) + ezy(i,j);
            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
            
            %    hx(i,j) = hx(i,j) - opH*( ez(i,j+1) - ez(i,j) );
            hx(i,j) = hx_decay(1,j)*hx(i,j) - hx_amp(1,j)*( ez(i,j+1) - ez(i,j) );
        end
    end

    
    for j=2:ydim-1
        for i=1:xdim-1
            
            %    hy(i,j) = hy(i,j) + opH*( ez(i+1,j) - ez(i,j) );
            hy(i,j) = hy_decay(1,i)*hy(i,j) + hy_amp(1,i)*( ez(i+1,j) - ez(i,j) );
        end
    end

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

toc

clear iZ0 Z0 Z0p2;
clear c0 eps0 mue0;
clear xdim ydim frequency xexc yexc;
clear ez ezx ezy hx hy idelta spread timeshift endofpulse tsteps;
clear xmin xmax ymin ymax;
clear xelength xhlength xekappa xhkappa;
clear yelength yhlength yekappa yhkappa;
clear DeltaT ezy_decay ezy_amp ezx_decay ezx_amp hx_decay hx_amp hy_decay hy_amp opE;
clear i j n centerfrequency delta deltaT exponent nsteps sigma_max t0 tw thickness z0p2 time;