%**********************************************************************
%   FDTD_2D ver.2.01
%**********************************************************************
%
%   Aug. 27, 1998  Written by M.Fujii
%   Apr. 28, 2011  Ported into Matlab code
%   Nov. 23, 2012  Replace x  with s  (source excitation)
%                          tx with ts (time data) 
%                          fx with fs (Fourier spectrum)
%                  Expanded to 2D code for TE mode
%
%   This program calculates 2-D plane wave propagating in xz-plain.
%
%   [Variables and parameters]
%
%     v:    MAIN VARIABLE INCLUDING BOTH E- AND H-FIELDS
%     s:    SOURCE EXICITATION
%     er:   RELATIVE PERMITTIVITY
%     mr:   RELATIVE PERMEABILITY
%     sg:   ELECTRIC CONDUCTIBITY
%     freq: FREQUENCY FOR EXCITATION
%     xmax: MAXIMUM DIMENSION IN X-DIRECTION
%     zmax: MAXIMUM DIMENSION IN Z-DIRECTION
%     tt:   PULSE DURATION
%     kd:   DIMENSION IN Z-DIRECTION
%
%   [Variables for Fourier transform]
%
%     tv:   TIME SERIES DATA OF E- AND H- FIELDS
%     ts:   TIME SERIES DATA OF EXCITATION "st"
%     fv:   Fourier transform of tv
%     fs:   Fourier transform of ts
%     fcnt: CENTER FREQUENCY
%     fspn: FREQUENCY SPAN
%     frsl: FREQUENCY RESOLUTION
%
%   [Field alignment and numbering]
%
%     In this program, the computation region is set from "kl" to "ku",
%     AND "il" to "iu".
%     The H field and the E field are stored in a single 2D array "v" 
%     in a staggered Yee cell alignment as shown in the following figure. 
%
%   Structure of "v"
%   [TE grid (o:E, +:H)]
%
%        x
%        |     Hz    Hz      ...               ... Hz    Hz
%        |     Ey Hx Ey Hx   ...               ... Ey Hx Ey
%        |     o--+--o--+--o--+--o--+--o--+--o--+--o--+--o       :Ey, Hx
%        |     |     |     |     |     |     |     |     |     
%   iup  |     +     +     +     +     +     +     +     +       :Hz
%        |     |     |     |     |     |     |     |     |     
%   iu   o-----o--+--o--+--o--+--o--+--o--+--o--+--o--+--o--+--o :Ey, Hx
%        |  A  |     |     |     |     |     |     |     |     |
%   ium  +  |  +     +     +     +     +     +     +     +     + :Hz
%        |  |  |     |     |     |     |     |     |     |     |
%        o--|--o--+--o--+--o--+--o--+--o--+--o--+--o--+--o--+--o :Ey, Hx
%        |  |  |     |     |     |     |     |     |     |     |
%        |  |  |     |     |     |     |     |     |     |     |
%        .  |  .     .     .     .     .     .     .     .     .
%        . id  .     .     .     .     .     .     .     .     .
%        |  |  |     |     |     |     |     |     |     |     |
%        |  |  |     |     |     |     |     |     |     |     |
%      5 o=-|--o--+--o--+--o--+--o--+--o--+--o--+--o--+--o--+--o :Ey, Hx
%        |  |  |     |     |     |     |     |     |     |     |
%      4 +  |  +     +     +     +     +     +     +     +     + :Hz
%        |  V  |     |     |     |     |     |     |     |     |
%   il 3 o-----o--+--o--+--o--+--o--+--o--+--o--+--o--+--o--+--o :Ey, Hx
%              |     |     |     |     |     |     |     |     
%      2       +     +     +     +     +     +     +     +       :Hz
%              |     |     |     |     |     |     |     |     
%   i: 1       o--+--o--+--o--+--o--+--o--+--o--+--o--+--o-------> z
%
%     k: 1  2  3  4  5  6    ...                ...  kum ku    ku+kl-1
%              kl                                        ku
%              |<------------------kd------------------->|
% real         |                              		     |
% space index: 0  1  2  3  4  ...                        kd 
%
%
%   [TM grid] to be written
%
%
%***********************************************************************

clear all
clf reset

%***********************************************************************
%     Fundamental constants
%***********************************************************************

cvac = 2.99792458e8;            %speed of light in vacuum free space 
mu0  = 4.0*pi*1.0e-7;           %permeability of free space
ep0  = 1.0/(cvac*cvac*mu0);     %permittivity of free space

freq   = 1.0e+9;                %frequency of source excitation
lambda = cvac/freq;             %wavelength of source excitation
omega  = 2.0*pi*freq;           %angular frequency

%***********************************************************************
%     Grid parameters
%***********************************************************************

id = 80; % 2x(number of Yee cells) in x-direction
kd = 80; % 2x(number of Yee cells) in z-direction

% indices
il  = 3;
ilp = il+1;   %il plus 1
iu  = il+id;
ium = iu-1;   %iu minus 1
ium2= iu-2;
kl  = 3;
klp = kl+1;   %kl plus 1
ku  = kl+kd;
kum = ku-1;   %ku minus 1
kum2= ku-2;

%***********************************************************************
%     Length of the space [m]
%***********************************************************************

xmax = 1.0;
zmax = 1.0; 

%***********************************************************************
%     Initialization of arrays for fields and excitation
%***********************************************************************

v = zeros(iu+il-1,ku+kl-1); %row vector of fields (E and H)
s = zeros(iu+il-1,ku+kl-1); %row vector of source excitation (J) at E-node

%***********************************************************************
%     Material parameters
%***********************************************************************

er = 1.0;
mr = 1.0;
sg = 1.0e-3;

%***********************************************************************
%    Excitation (specified at E-node, i.e. kl+(even number))
%***********************************************************************

s(il+id/2,kl+kd/2) = 1.0; %at center
%s(kl+10)  = 1.0; %at left-hand side

%***********************************************************************
%     Field detection position
%***********************************************************************

idtct = il+id/2;
kdtct = kl+kd/2;

%***********************************************************************
%     Time and space conditions
%***********************************************************************

dx = xmax/(id/2);
dz = zmax/(kd/2); %space step  (Yee's cell size)
dt = 1.0/cvac/sqrt((1/dx)^2+(1/dz)^2);     %time step (maximum Courant limit)

nmax = 100;       %number of time steps

%***********************************************************************
%     Excitation pulse period
%***********************************************************************

tt     = 1.0/freq;       %period
npulse = round(tt/dt/2); %number of pulse excitation steps

%***********************************************************************
%     Updating coefficients for space region with nonpermeable media
%***********************************************************************

dnm = 2*er*ep0 + sg*dt;
c0 = (2*er*ep0 - sg*dt)/dnm;
ce = 2*dt/dnm;
%cs = 2*dt/dnm;          %excitation by current (soft source)
cs = 1;                  %excitation by E-field (soft source)
cm = dt/(mr*mu0);

%***********************************************************************
%     Plot excitation
%***********************************************************************

x = linspace(0,xmax,1+id); %x-coordinates
z = linspace(0,zmax,1+kd); %z-coordinates

figure(1)

% Ey
subplot(2,2,1),pcolor(z(1:2:kd+1),x(1:2:id+1),s(il:2:iu,kl:2:ku));
shading interp;
axis([0 zmax 0 xmax]);
caxis([-1 1]);
title('Excitation Jx [A/m]');

% Hx
subplot(2,2,2),pcolor(z(2:2:kd),x(1:2:id+1),s(il:2:iu,klp:2:kum));
shading interp;
axis([0 zmax 0 xmax]);
caxis('auto');
title('Hx [A/m]');

% Hz
subplot(2,2,3),pcolor(z(1:2:kd+1),x(2:2:id),s(ilp:2:ium,kl:2:ku));
shading interp;
axis([0 zmax 0 xmax]);
caxis('auto');
title('Hz [A/m]');

%***********************************************************************
%***********************************************************************
%     Start FDTD time-marching loop
%***********************************************************************
%***********************************************************************

%***********************************************************************
%     Initialization of time data arrays
%***********************************************************************

tv = zeros(1,nmax);
ts = zeros(1,nmax);

for n=1:nmax

%***********************************************************************
%     Excitation time signal (choose "pulse" or "cw")
%***********************************************************************

%******
% Pulse
%******
%if n <= npulse                 %half-sine period
%    st = -sin(pi*n/npulse)^2;  %raised cosine pulse
%    %st = 1;                   %not work (too abrupt)
%else
%    st = 0;
%end

%*********************
% Continuous wave (cw)
%*********************
if n <= npulse/2               %quarter-sine period
    st = -sin(pi*n/npulse)^2;  %raised-cosine for smooth wave front
else  
    st = -sin(pi*n/npulse);    %continuous sine wave
end

%***********************************************************************
%     Update electric fields
%***********************************************************************

%Ey (Ampere's law)
v(il:2:iu,kl:2:ku) ...
  = c0*v(il:2:iu,kl:2:ku) ...
  + ce*(  (v(il:2:iu,kl+1:2:ku+1)-v(il:2:iu,kl-1:2:ku-1))/dz ...
        - (v(il+1:2:iu+1,kl:2:ku)-v(il-1:2:iu-1,kl:2:ku))/dx ...
        - cs*st*s(il:2:iu,kl:2:ku));

% PEC boundary
v(il,kl:2:ku) = zeros; 
v(iu,kl:2:ku) = zeros;
v(il:2:iu,kl) = zeros;
v(il:2:iu,ku) = zeros;

%***********************************************************************
%     Update magnetic fields
%***********************************************************************

%Hx (Faraday's law)
v(il:2:iu,klp:2:kum) ...
  = v(il:2:iu,klp:2:kum) ...
  + cm*(v(il:2:iu,klp+1:2:kum+1)-v(il:2:iu,klp-1:2:kum-1))/dz; 

%Hz (Faraday's law)
v(ilp:2:ium,kl:2:ku) ...
  = v(ilp:2:ium,kl:2:ku) ...
  - cm*(v(ilp+1:2:ium+1,kl:2:ku)-v(ilp-1:2:ium-1,kl:2:ku))/dx; 

%***********************************************************************
%     Visualizing fields
%***********************************************************************

figure(2)

% Ey
subplot(2,2,1),pcolor(z(1:2:kd+1),x(1:2:id+1),v(il:2:iu,kl:2:ku));
shading interp;
axis([0 zmax 0 xmax]);
caxis([-3 3]);
title(['Ey [V/m] at Time step = ',num2str(n)]);

% Hx
subplot(2,2,2),pcolor(z(2:2:kd),x(1:2:id+1),v(il:2:iu,klp:2:kum));
shading interp;
axis([0 zmax 0 xmax]);
caxis([-5e-3 5e-3]);
title('Hx [A/m]');

% Hz
subplot(2,2,3),pcolor(z(1:2:kd+1),x(2:2:id),v(ilp:2:ium,kl:2:ku));
shading interp;
axis([0 zmax 0 xmax]);
caxis([-5e-3 5e-3]);
title('Hz [A/m]');

M(:,n)=getframe; %get frame for movie

%***********************************************************************
%     End of the time-marching loop
%***********************************************************************

time(n) = n*dt;
tv(n)   = v(idtct,kdtct);
ts(n)   = st;

end

%***********************************************************************
%     Plot time and frequency data
%***********************************************************************

%******************
% Fourier transform
%******************

tmax = time(nmax); %Maximum time
df   = 1/tmax;     %Frequency resolution
fmax = df*nmax;    %Maximum frequency
fs = abs(fft(ts,nmax));
fv = abs(fft(tv,nmax));

%*****
% Plot
%*****

figure(3)

subplot(2,1,1),plot(time,ts,'-r', time,tv,':b');
title('Time data')
xlabel('Time [s]');ylabel('Ex [V/m] or Hy [A/m]');
legend ('Excitation','Response')
legend boxoff
grid on

subplot(2,1,2),plot(df:df:fmax/2,fs(1:1:nmax/2),'-r', df:df:fmax/2,fv(1:1:nmax/2),':b');
%subplot(2,1,2),plot(df:df:fmax,fs,'r',df:df:fmax,fv,'b'); %Full range of frequency
title('Frequency spectra')
xlabel('Frequency [Hz]');ylabel('Ex or Hy [arb. unit]');
legend ('Excitation','Response')
legend boxoff
grid on

