tmz_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
time = 2e-9;                    % scanned time
deltaT = delta/(c0*sqrt(2));    % timestep
nsteps = ceil(time/deltaT)+1;

radius = c0*time;                   % distance which light propagtes for the time t
xdim = ceil((2*radius)/delta) + 20; % dimension in x direction
ydim = xdim;                        % dimension in y direction


% frequency = 24e9;          % frequency bandwidth for the gauss impulse
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
hx = zeros(xdim,ydim);      % hx field
hy = zeros(xdim,ydim);      % hy field

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


%% calculation of important parameters

iZ0 = sqrt(eps0/mue0);
Z0 = sqrt(mue0/eps0);
idelta = 1./delta;
deltaT = delta/(c0*sqrt(2));
opH = (deltaT/(mue0*delta));
opE = (deltaT/(eps0*delta));

spread=((pi*frequency)^2)/log(10.);
timeshift=sqrt((5*log(10.))/spread);
endofpulse = 3*timeshift;

%% calculation of the time loop

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
    %             ez(i,j) = ez(i,j) + opE*( hy(i,j) - hy(i-1,j) - hx(i,j) + hx(i,j-1) );
    %         end
    %     end
    ez(2:end-1,2:end-1) = ez(2:end-1,2:end-1) + opE*( hy(2:end-1,2:end-1) - hy(1:end-2,2:end-1) - hx(2:end-1,2:end-1) + hx(2:end-1,1:end-2) );

    
    % soft excitation
    % ez(xexc,yexc) = ez(xexc,yexc) - opE*exp(-((((n-1)*deltaT-timeshift)^2)*spread));
    % 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*sin(2*pi*frequency*(n-1)*deltaT);

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

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

    
    % visualisation
%     if (mod(n,5)==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, nx = %i, ny = %i, step = %i',xdim,ydim,n));
%         title(sprintf('FDTD Time = %.3f nsec = %.2f %%',n*deltaT*1e9,((n*deltaT)/(nsteps*deltaT))*100),'Color',[1 0 0],'FontSize', 22); drawnow;
%     end
    
end

toc

clear hx hy ez Z0 c0 endofpulse eps0 mue0 iZ0 idelta n opE opH radius spread xdim xexc ydim yexc;
clear delta deltaT t0 tw frequency centerfrequency time timeshift nsteps;