%***********************************************************************%     1-D FDTD code with E-field boundary condition%***********************************************************************%%     Program author: Mats Gustafsson%                     Department of Electroscience%                     Lund Institute of Technology%                     Sweden%                     %     Date of this version:  October 2002%%%   Ein%   --->                Eut                         PEC%     ----------------------------------------------->%    z=0               z=Lut                        z=Lz%    r=1                                            r=Nz+1%%***********************************************************************% Physical constantsmu0  = 4e-7 * pi;             % Permeability of vacuumc0   = 299792458;             % Speed of light in vacuumeps0 = 1/c0^2/mu0;            % Permittivity of vacuumeta0 = sqrt(mu0/eps0);        % Wave impedance in vacuum  GHz = 10^9;mm = 10^(-3);%***********************************************************************% Parameter initiation%***********************************************************************Lz = 3;                     % Length in metersTmax = Lz/c0;                 % Time interval [0,Tmax]Dz = 20*mm;                    % spatial grid sizeR = 0.9;                  % stabillity number (grid cells per time step)freq = 1*GHz;                  % center frequency of source excitationfmax = 4*GHz;                  % Max frequency to plotLut = Lz/2;                        % Output pointn_pict = 20;                    % plot each n_pict step Nfft = 4*1024;                    % Number of points for the fftNz = round(Lz/Dz);                     % Number of cells in the z-directionNut = round(Lut/Dz);           % Coordinate for outputDt = Dz*R/c0;                  % Time stepNt = round(Tmax/Dt);           % Number of time stepsz = linspace(0,Lz,Nz+1);         % z-coordinatest = Dt*[0:Nt-1]';               % timelambda = c0/freq;               % center wavelength of source excitationppw = lambda/Dz                 % sample points per wave length at the center frequencyomega = 2.0*pi*freq;            % angular frequency%***********************************************************************% Allocate field vectors%***********************************************************************Ex = zeros(Nz+1  ,1);          % Electric fieldHy = zeros(Nz,1);          % Magnetic fieldEut = zeros(Nt,1);           % Outputres = zeros(Nz+1,2*Nt);%***********************************************************************%     Wave excitation, Incident puls from the left%***********************************************************************% Broad band pulstau = 250.0e-12;delay = 2.5*tau;Ein = sin(omega*(t-delay)).*exp(-((t-delay).^2/tau^2));% Incident time harmonic wave from the left%Ein = sin(omega*(t));Ehin = fft(Ein,Nfft)/Nt*2;Df = 1/(Dt*Nfft);freqN = round(fmax/Df);f = linspace(0,Df*freqN,freqN);% analytic solutionpropdelay = (Nut-1)*Dz/c0;Eut_ana = sin(omega*(t-delay-propdelay)).*exp(-((t-delay-propdelay).^2/tau^2)).*(sign(t-propdelay)+1)/2;%Eut_ana = sin(omega*(t-propdelay)).*(sign(t-propdelay)+1)/2;Ehut_ana = fft(Eut_ana,Nfft)/Nt*2;%***********************************************************************% Time stepping%***********************************************************************for n = 1:Nt;   % Update H from n*Dt-Dt/2 to n*Dt+Dt/2  Hy = Hy - (Dt/mu0/Dz) * diff(Ex,1);                     % main grid    % Update E from n*Dt to (n+1)*Dt  Ex(2:Nz) = Ex(2:Nz) - (Dt/eps0/Dz) * diff(Hy,1);    % main grid except boundary    Ex(1) = Ein(n);          % E-field on the boundary      % Sample the electric field at chosen points  if mod(n,n_pict) == 0      figure(1);      subplot(2,1,1);      plot(z,Ex(1:Nz+1),'LineWidth',1);      title(['time is ',num2str(n*Dt*10^9),'ns'])      axis tight;      ylabel('E-field');            subplot(2,1,2);      plot(z(1:Nz),Hy(1:Nz),'LineWidth',1);      title(['time is ',num2str(n*Dt*10^9),'ns'])      axis tight;      ylabel('H-field');      pause;%(0.05);  end  Eut(n) = Ex(Nut);end%***********************************************************************%     Post processing%***********************************************************************Eerror = norm(Eut-Eut_ana)*Dt   % L2 errorfigure(2);clf;subplot(3,1,1);plot(t,Eut,t,Eut_ana,'r','LineWidth',1)axis([0.9*propdelay 1.9*propdelay -1 1])xlabel('time');ylabel('E-field');subplot(3,1,2);warning off;A=plotyy(f,abs(Ehin(1:freqN)),f,c0./f/Dz); warning on;set(A(2),'yLim',[0 30],'yTick',[0:5:30],'xlim',[0 Df*freqN],'xGrid','on','yGrid','on');set(A(1),'xlim',[0 Df*freqN]);set(get(A(2),'Ylabel'),'String','points/wavelength');set(get(A(1),'Ylabel'),'String','input');set(get(A(1),'Xlabel'),'String','frequency');subplot(3,1,3);Ehut = fft(Eut,Nfft)/Nt*2;%plot(f,20*log10(abs(abs(Ehut(1:freqN))-abs(Ehut_ana(1:freqN)))./abs(Ehut_ana(1:freqN))));Eherror = abs(Ehut(1:freqN)-Ehut_ana(1:freqN))./abs(Ehut_ana(1:freqN));Eherror1 = interp1(f,Eherror,freq)plot(f,20*log10(Eherror),'LineWidth',1);axis tightax=axis;axis([ax(1:2) -120 20]);xlabel('frequency');ylabel('Error in dB');