%**********************************************************************
%   FDTD_2D ver.2.11
%**********************************************************************
%
%   Aug. 27, 1998  1D Fortran code written by M.Fujii
%   Apr. 28, 2011  Ported into Matlab code
%   Nov. 23, 2012  Expanded to 2D TE mode
%                  Replace x  with s  (source excitation)
%                          tx with ts (time data) 
%                          fx with fs (Fourier spectrum)
%   Jan. 15, 2013  APML implemented for 2D TE mode
%
%   This program calculates 2-D plane wave propagating in xz-plain.
%
%   [Variables and parameters]
%
%     v1:   MAIN VARIABLE INCLUDING BOTH E- AND H-FIELDS
%     v2:   MAIN VARIABLE INCLUDING BOTH D- AND B-FLUXES
%     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 "v1" 
%     in a staggered Yee cell alignment as shown in the following figure.
%     Storage of B and D fluxes follows those of H and E fields, respectively.
%
%   Structure of "v"
%   [TE grid (o:E,D, +:H,B)]
%
%        x
%        |     Hz    Hz      ...               ... Hz    Hz
%        |     Ey Hx Ey Hx   ...               ... Ey Hx Ey
%   im   |     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    km=ku+kl-1
%              kl                                        ku    
%              |<------------------kd------------------->|
% real         |                              		     |
% space index: 0  1  2  3  4  ...                        kd 
%
%***********************************************************************

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   = 0.7e+9;                %frequency of source excitation (under cut-off)
freq   = 1.0e+9;                %frequency of source excitation (above cut-off)
%(Note: Xmax = 0.2m, Cut-off frequency = 0.75GHz) 

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

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

id = 40; % 2x(number of Yee cells) in x-direction
kd = 100; % 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;
im  = iu+il-1; %i_max

kl  = 3;
klp = kl+1;    %kl plus 1
ku  = kl+kd;
kum = ku-1;    %ku minus 1
kum2= ku-2;
km  = ku+kl-1; %k_max

%***********************************************************************
%     Number of time steps
%***********************************************************************

nmax = 250;

%***********************************************************************
% set PML parameters 
% npml: 2x(Number of PML layers), npml=1 no-PML
% mpml: Order of PML variation, mpml=3~5
% nswpml?? = 0:off, 1:on
%***********************************************************************
npml = 8;     %8 to 20 is sufficient
nswpmlxl = 0; %1:PML on, 0:off
nswpmlxu = 0; %1:PML on, 0:off
nswpmlzl = 1; %1:PML on, 0:off
nswpmlzu = 1; %1:PML on, 0:off

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

xmax = 0.2; %(Cut-off freq = 0.75 GHz) 
zmax = 1.0; 

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

er(1:im,1:km) = 1.0;
mr(1:im,1:km) = 1.0;
sg(1:im,1:km) = 0.0; %Not used in the PML formulation

% Averaging for inhomogeneous material boundary
% Ey node
er(il:2:iu,kl:2:ku) = 0.25*(er(il-1:2:iu-1,kl-1:2:ku-1) ...
                           +er(il-1:2:iu-1,kl+1:2:ku+1) ...
                           +er(il+1:2:iu+1,kl-1:2:ku-1) ...
                           +er(il+1:2:iu+1,kl+1:2:ku+1));
sg(il:2:iu,kl:2:ku) = 0.25*(sg(il-1:2:iu-1,kl-1:2:ku-1) ...
                           +sg(il-1:2:iu-1,kl+1:2:ku+1) ...
                           +sg(il+1:2:iu+1,kl-1:2:ku-1) ...
                           +sg(il+1:2:iu+1,kl+1:2:ku+1));

% Hx node
mr(il:2:iu,klp:2:kum) = 2.0*mr(il-1:2:iu-1,klp:2:kum) ...
                          .*mr(il+1:2:iu+1,klp:2:kum) ...
                         ./(mr(il-1:2:iu-1,klp:2:kum) ...
                           +mr(il+1:2:iu+1,klp:2:kum));

% Hz node
mr(ilp:2:ium,kl:2:ku) = 2.0*mr(ilp:2:ium,kl-1:2:ku-1) ...
                          .*mr(ilp:2:ium,kl+1:2:ku+1) ...
                         ./(mr(ilp:2:ium,kl-1:2:ku-1) ...
                           +mr(ilp:2:ium,kl+1:2:ku+1));
     
%***********************************************************************
%    Excitation: specify at E-node, i.e. il+(even num.), kl+(even num.)
%***********************************************************************

s = zeros(im,km); %initialization of source excitation (J) array
%s(il+id/2,kl+kd/2) = 1.0; %at center
s(il+id/2,kl+npml+2) = 1.0; %at left-end

%***********************************************************************
%     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)

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

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

%***********************************************************************
%     PML parameters and update equation coefficients
%***********************************************************************

% PML material parameters
mpml  = 4; % 4 to 5 is recommended
sgopt = (mpml+1)/(150*pi);
sgmxl = sgopt/dx*nswpmlxl;
sgmxu = sgopt/dx*nswpmlxu;
sgmzl = sgopt/dz*nswpmlzl;
sgmzu = sgopt/dz*nswpmlzu;
akmxl = 2.0; % 1.0 - 20.0 recommended
akmxu = 2.0;
akmzl = 2.0;
akmzu = 2.0;

% PML space variations
fpos = zeros(1,npml); %increasing order
fneg = zeros(1,npml); %decreasing order
fpos(1:npml) = ((1:npml)./npml).^mpml;
fneg(1:npml) = ((npml-(0:npml-1))./npml).^mpml;
akx = ones (1,im);
sgx = zeros(1,im);
aky = 1.0;
sgy = 0.0;
akz = ones (1,km);
sgz = zeros(1,km);

% PML coefficients
% X
if (nswpmlxl == 1)
  sgx(il:il+npml-1) = sgmxl*fneg(1:npml);
  akx(il:il+npml-1) = (akmxl-1)*fneg(1:npml)+1;
end
if (nswpmlxu ==1)
  sgx(iu-npml+1:iu) = sgmxu*fpos(1:npml);
  akx(iu-npml+1:iu) = (akmxu-1)*fpos(1:npml)+1;
end
% Z
if (nswpmlzl == 1)
  sgz(kl:kl+npml-1) = sgmzl*fneg(1:npml);
  akz(kl:kl+npml-1) = (akmzl-1)*fneg(1:npml)+1;
end
if (nswpmlzu ==1)
  sgz(ku-npml+1:ku) = sgmzu*fpos(1:npml);
  akz(ku-npml+1:ku) = (akmzu-1)*fpos(1:npml)+1;
end

% 1D coefficients for update equations

cxp(il:iu) = 2.0*ep0*akx(il:iu) + dt*sgx(il:iu);
cxm(il:iu) = 2.0*ep0*akx(il:iu) - dt*sgx(il:iu);
c0x(il:iu) = cxm(il:iu)./cxp(il:iu);
c1x(il:iu) = 2.0*ep0*dt./cxp(il:iu);

cyp = 2.0*ep0*aky + dt*sgy;
cym = 2.0*ep0*aky - dt*sgy;
c0y = cym/cyp;
c1y = 2.0*ep0*dt/cyp;

czp(kl:ku) = 2.0*ep0*akz(kl:ku) + dt*sgz(kl:ku);
czm(kl:ku) = 2.0*ep0*akz(kl:ku) - dt*sgz(kl:ku);
c0z(kl:ku) = czm(kl:ku)./czp(kl:ku);
c1z(kl:ku) = 2.0*ep0*dt./czp(kl:ku);

% 3D coefficients for update equations
% initialization
ce=zeros(im,km);
cm=zeros(im,km);
% Ey node
ce(il:2:iu, kl:2:ku)  = 1.0./(ep0*er(il:2:iu, kl:2:ku)...
          .* repmat(cxp(il:2:iu)',1,(ku-kl)/2+1));  %repeat matrix: see help
% Hx node
cm(il:2:iu, klp:2:ku) = 1.0./(mu0*mr(il:2:iu, klp:2:kum)...
          .* repmat(czp(klp:2:kum),(iu-il)/2+1,1)); %repeat matrix: see help
% Hz node
cm(ilp:2:ium,kl:2:ku) = 1.0./(mu0*mr(ilp:2:ium, kl:2:ku)*cyp);

%***********************************************************************
%     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]');

%***********************************************************************
%     Initialization of field arrays
%***********************************************************************

v1= zeros(im,km); %E and H fields for PML
v2= zeros(im,km); %D and B fluxes for PML

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

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

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

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 of the standard FDTD equations)
%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));

%Ey for PML
vd = v2(il:2:iu, kl:2:ku);
v2(il:2:iu, kl:2:ku) = repmat(c0z(kl:2:ku),(iu-il)/2+1,1)...  %Dy
                    .* v2(il:2:iu, kl:2:ku)...
                     + repmat(c1z(kl:2:ku),(iu-il)/2+1,1)...  %repeat matrix
                    .* ((v1(il:2:iu,kl+1:2:ku+1)-v1(il:2:iu,kl-1:2:ku-1))/dz...
                      - (v1(il+1:2:iu+1,kl:2:ku)-v1(il-1:2:iu-1,kl:2:ku))/dx...
                      - st*s(il:2:iu,kl:2:ku));
v1(il:2:iu, kl:2:ku) = repmat(c0x(il:2:iu)',1,(ku-kl)/2+1)... %Ey
                    .* v1(il:2:iu, kl:2:ku)...
                     + ce(il:2:iu, kl:2:ku)...
                    .* (cyp*v2(il:2:iu, kl:2:ku) - cym*vd);
                
% PEC boundary
v1(il,kl:2:ku) = zeros; 
v1(iu,kl:2:ku) = zeros;
v1(il:2:iu,kl) = zeros;
v1(il:2:iu,ku) = zeros;

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

%Hx (Faraday's law of the standard FDTD equations)
%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; 

%Hx for PML
vb = v2(il:2:iu, klp:2:kum);
v2(il:2:iu, klp:2:kum) = c0y...                                  %Bx
                       * v2(il:2:iu, klp:2:kum)...
                       + c1y...
                       * (v1(il:2:iu,klp+1:2:kum+1)-v1(il:2:iu,klp-1:2:kum-1))/dz;
v1(il:2:iu, klp:2:kum) = repmat(c0z(klp:2:kum),(iu-il)/2+1,1)... %Hx
                      .* v1(il:2:iu, klp:2:kum)...
                       + cm(il:2:iu, klp:2:kum)...
                      .* (repmat(cxp(il:2:iu)',1,(kum-klp)/2+1).*v2(il:2:iu, klp:2:kum)...
                       -  repmat(cxm(il:2:iu)',1,(kum-klp)/2+1).*vb);

%Hz (Faraday's law of the standard FDTD equations)
%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; 

%Hz for PML
vb = v2(ilp:2:ium,kl:2:ku);
v2(ilp:2:ium,kl:2:ku) = repmat(c0x(ilp:2:ium)',1,(ku-kl)/2+1)... %Bz
                     .* v2(ilp:2:ium,kl:2:ku)...
                      - repmat(c1x(ilp:2:ium)',1,(ku-kl)/2+1)...
                     .* (v1(ilp+1:2:ium+1,kl:2:ku)-v1(ilp-1:2:ium-1,kl:2:ku))/dx;
v1(ilp:2:ium,kl:2:ku) = c0y...                                   %Hz
                      * v1(ilp:2:ium,kl:2:ku)...
                      + cm(ilp:2:ium,kl:2:ku)...
                     .* (repmat(czp(kl:2:ku),(ium-ilp)/2+1,1).*v2(ilp:2:ium,kl:2:ku)...
                      -  repmat(czm(kl:2:ku),(ium-ilp)/2+1,1).*vb);

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

figure(2)

% Ey
subplot(2,2,1),pcolor(z(1:2:kd+1),x(1:2:id+1),v1(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),v1(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),v1(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
getframe; %get frame for movie

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

time(n) = n*dt;
tv(n)   = v1(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

