%***********************************************************************
% 1-D FDTD code with initial condition E(z,0) = E0
% Illustrates the staggered grid i FDTD
%***********************************************************************
%
% Program author: Mats Gustafsson
% Department of Electroscience
% Lund Institute of Technology
% Sweden
%
% Date of this version: October 2002
%
%
%
% PEC PEC
% ----------------------------------------------->
% z=0 z=Lz
% r=1 r=Nz+1
%
%***********************************************************************
% Physical constants
mu0 = 4e-7 * pi; % Permeability of vacuum
c0 = 299792458; % Speed of light in vacuum
eps0 = 1/c0^2/mu0; % Permittivity of vacuum
eta0 = sqrt(mu0/eps0); % Wave impedance in vacuum
GHz = 10^9;
mm = 10^(-3);
%***********************************************************************
% Parameter initiation
%***********************************************************************
Lz = 1; % Length in meters
Tmax = Lz/c0; % Time interval [0,Tmax]
Dz = 20*mm; % spatial grid size
R = 1; % stabillity number (grid cells per time step)
freq = 1*GHz; % center frequency of source excitation
n_pict = 2; % plot each n_pict step
Nz = round(Lz/Dz); % Number of cells in the z-direction
Dt = Dz*R/c0; % Time step
Nt = round(Tmax/Dt); % Number of time steps
z = linspace(0,Lz,Nz+1); % z-coordinates
t = Dt*[0:Nt-1]'; % time
lambda = c0/freq; % center wavelength of source excitation
ppw = lambda/Dz % sample points per wave length at the center frequency
omega = 3.0*pi*freq; % angular frequency
% colormap(ho)
hott=hot(32);
hottt=hott(size(hott,1):-1:1,:);
ho=[1-hottt; hottt];
%***********************************************************************
% Allocate field vectors
%***********************************************************************
Ex = zeros(Nz+1 ,1); % Electric field
Hy = zeros(Nz,1); % Magnetic field
res = zeros(2*Nz+1,2*Nt);
%***********************************************************************
% Wave excitation, initial condition
%***********************************************************************
E0 = sin(omega*(z-Lz/2)/c0).*exp(-omega^2*((z-Lz/2).^2/c0^2));
Ex = E0';
%***********************************************************************
% Time stepping
%***********************************************************************
for n = 1:Nt;
res1 = zeros(2*Nz+1,2*Nt);
% Update H from n*Dt-Dt/2 to n*Dt+Dt/2
Hy = Hy - (Dt/mu0/Dz) * diff(Ex,1); % main grid
res(2:2:2*Nz,2*n-1) = eta0*Hy;
res1(2:2:2*Nz,2*n-1) = eta0*Hy;
% 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
res(1:2:2*Nz+1,2*n) = Ex;
res1(1:2:2*Nz+1,2*n) = Ex;
% Sample the electric field at chosen points
if mod(n,n_pict) == 0
figure(2);
subplot(2,1,1);
ma = max(max(res));
imagesc(res',[-ma ma])
colormap(ho)
colorbar
axis xy
title(['time is ',num2str(n*Dt*10^9),'ns'])
axis tight;
xlabel('2*space/Dz');
ylabel('2*time/Dt');
subplot(2,1,2);
ma = max(max(res1));
imagesc(res1',[-ma ma])
colormap(ho)
colorbar
axis xy
title(['time is ',num2str(n*Dt*10^9),'ns'])
axis tight;
xlabel('2*space/Dz');
ylabel('2*time/Dt');
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');
xlabel('z/m');
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');
xlabel('z/m');
pause(0.05);
end
end