%% TMz simulation for the FDTD method% close all;clear all;c0 = 2.99792458e8;eps0 = 8.854187818e-12;mue0 =  4*pi*1e-7;%% user definable parametersthickness = 10;              % thickness of the pml absorberexp_sigma = 4;               % exponent of the polynomial gradingdelta = 0.001;              % distance from discretisation line to linetime = 2e-9;                    % scanned timedeltaT = delta/(c0*sqrt(2));    % timestepnsteps = ceil(time/deltaT)+1;xdim = 40+2*thickness;                  % dimension in x directionydim = 40+2*thickness;                  % dimension in y directioncenterfrequency = 12e9;frequency = 12e9;          % frequency bandwidth for the gauss impulsexexc = ceil((xdim)/2);        % setting the x-coordinate for the excitationyexc = 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));%% allocation of memoryez = zeros(xdim,ydim);      % ez field, electric strengthez_x_psi = zeros(xdim,ydim);ez_y_psi = zeros(xdim,ydim);hx = zeros(xdim,ydim);      % hx field, magnetic strengthhx_y_psi = zeros(xdim,ydim);hy = zeros(xdim,ydim);      % hy field, magnetic strengthhy_x_psi = zeros(xdim,ydim);%% calculation of important parametersidelta = 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 sigmaxmin = thickness*delta;                 % first position of the pml absorber in x-directionxmax = (xdim-thickness-1)*delta;          % last position of the pml absorber in x-directionymin = thickness*delta;                 % first position of the pml absorber in y-directionymax = (ydim-thickness-1)*delta;          % last position of the pml absorber in y-directionxelength = zeros(1,xdim);   % absolute length on the electric grid in x-directionxhlength = zeros(1,xdim);   % absolute length on the magnetic grid in x-directionxesigma = zeros(1,xdim);    % grading of sigma on the electric grid in x-directionxhsigma = zeros(1,xdim);    % grading of sigma of the magnetic grid in x-directionyelength = zeros(1,ydim);   % absolute length on the electric grid in y-directionyhlength = zeros(1,ydim);   % absolute length on the magnetic grid in y-directionyesigma = zeros(1,ydim);    % grading of sigma on the electric grid in y-directionyhsigma = zeros(1,ydim);    % grading of sigma on the magnetic grid in y-directionif 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;        end        for i=xdim-thickness:xdim            xesigma(1,i)=sigma_max*(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;        end        for i=xdim-thickness+1:xdim            xhsigma(1,i)=sigma_max*(abs(xhlength(1,i)-xmax)/((thickness)*delta))^exp_sigma;        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;        end        for i=ydim-thickness:ydim            yesigma(1,i)=sigma_max*(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;        end        for i=ydim-thickness+1:ydim            yhsigma(1,i)=sigma_max*(abs(yhlength(1,i)-ymax)/((thickness)*delta))^exp_sigma;        end        end%% calculation of the material operatorsDeltaT = 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*xesigma(1,i))/eps0));        ez_x_psi_amp(1,i) = ( ez_x_psi_decay(1,i) - 1 );    end    if i >= xdim-thickness+1        ez_x_psi_decay(1,i) = exp(-((DeltaT*xesigma(1,i))/eps0));        ez_x_psi_amp(1,i) = ( ez_x_psi_decay(1,i) - 1 );    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*xhsigma(1,i))/eps0));        hy_x_psi_amp(1,i) = ( hy_x_psi_decay(1,i) - 1 );    end    if (i >= xdim-thickness+1) && (i~=xdim)        hy_x_psi_decay(1,i) = exp(-((DeltaT*xhsigma(1,i))/eps0));        hy_x_psi_amp(1,i) = ( hy_x_psi_decay(1,i) - 1 );    end    endfor j=1:ydim        if j <= thickness        ez_y_psi_decay(1,j) = exp(-((DeltaT*yesigma(1,j))/eps0));        ez_y_psi_amp(1,j) = ( ez_y_psi_decay(1,j) - 1 );    end        if j >= ydim-thickness+1        ez_y_psi_decay(1,j) = exp(-((DeltaT*yesigma(1,j))/eps0));        ez_y_psi_amp(1,j) = ( ez_y_psi_decay(1,j) - 1 );    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*yhsigma(1,j))/eps0));        hx_y_psi_amp(1,j) = ( hx_y_psi_decay(1,j) - 1 );    end        if (j >= xdim - thickness + 1)  && (j~=ydim)        hx_y_psi_decay(1,j) = exp(-((DeltaT*yhsigma(1,j))/eps0));        hx_y_psi_amp(1,j) = ( hx_y_psi_decay(1,j) - 1 );    endendclear xelength xhlength xekappa xhkappa xesigma xhsigma xealpha xhalpha;clear yelength yhlength yekappa yhkappa yesigma yhsigma yealpha yhalpha;%% allocate array for examinationstime = zeros(1,nsteps);excitation = zeros(1,nsteps);uxdirect = zeros(1,nsteps);uydirect = zeros(1,nsteps);uxycorner = zeros(1,nsteps);%% calculation of the time looptmp_dx = 0;tmp_dy = 0;opExc = opE*idelta;tw = 26.53e-12;t0 = 4*tw;ticfor 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*( tmp_dx - 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 * ( tmp_dy + hx_y_psi(i,j) );%             hx(i,j) = hx_decay(1,i)*hx(i,j) - hx_amp(1,i) * ( tmp_dy +%             hx_y_psi(i,j) );%             hx(i,j) = hx_decay(1,i)*hx(i,j) - hx_amp(1,i) * 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 * ( tmp_dx + hy_x_psi(i,j) );%             hy(i,j) = hy_decay(1,j)*hy(i,j) + hy_amp(1,j) * ( tmp_dx + hy_x_psi(i,j) );%             hy(i,j) = hy_decay(1,j)*hy(i,j) + hy_amp(1,j) * 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 Nearly perfect PML, step = %i', n));%         title(sprintf('FDTD Time (NPML) = %.3f nsec',n*deltaT*1e9),'Color',[1 0 0],'FontSize', 22); drawnow;%     endendtocclear 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;