%   FDTD1_3D FDTD - ECE529
%   First-order ABCs and super absorption
clear all;
tic
eps0      = 8.85418782e-012;               %  ANSOFT HFSS value 
mu0       = 1.25663706e-006;               %  ANSOFT HFSS value
c0        = 1/sqrt(eps0*mu0);              %  ANSOFT HFSS value

scrsz               = get(0,'ScreenSize');
figure('Position',[0.05*scrsz(3) 0.05*scrsz(4) 0.8*scrsz(3) 0.8*scrsz(4)]);
colormap(jet(128)); 
%------------------------------------------------------------------
%------------------------------------------------------------------
% 3D geometry definitions and magic time step
%------------------------------------------------------------------
%------------------------------------------------------------------
xmin = -0.4;              %   domain in meters
xmax = +0.4;              %   domain in meters 
ymin = -0.4;              %   domain in meters
ymax = +0.4;              %   domain in meters
zmin = -0.5;              %   domain in meters
zmax = +0.5;              %   domain in meters

d = 0.0100;             %   spatial step in meters                              
x = [xmin:d:xmax];      %   x variable
y = [ymin:d:ymax];      %   y variable
z = [zmin:d:zmax];      %   y variable

DUMMY   = ones(length(x), length(y), length(z));

SIGMA                   = 0.0*DUMMY;
EPS                     = ones(size(DUMMY)) + 0.0*DUMMY;    

Nx   = length(x) - 1;   %   Nx
Ny   = length(y) - 1;   %   Ny
Nz   = length(z) - 1;   %   Ny
dt   = 1/(c0*sqrt(1/d^2 + 1/d^2 +1/d^2));   %   Magic time step
%------------------------------------------------------------------
% Excitation - voltage(s)
%------------------------------------------------------------------
%------------------------------------------------------------------
f = 900e6; 
                                    %    (to travel from feed to the end)
T        = 16e-9;                   %    observation time
%   Number of time steps and time vector
NT = round(T/dt); t = [0: dt: NT*dt]; Energy = zeros(size(t)); 
%   Excitation (generator) voltage
VG = sin(2*pi*f*t);
IG = zeros(size(VG));
RG = 50;                                    %   default
bound = 4*ones(1, NT+1);                    %   graphics scale
%------------------------------------------------------------------
%------------------------------------------------------------------
% Allocate field matrices
%------------------------------------------------------------------
%------------------------------------------------------------------
ExP = zeros(Nx  , Ny+1, Nz+1);
EyP = zeros(Nx+1, Ny  , Nz+1);
EzP = zeros(Nx+1, Ny+1, Nz  );
HxP = zeros(Nx+1, Ny  , Nz  );
HyP = zeros(Nx  , Ny+1, Nz  );
HzP = zeros(Nx  , Ny  , Nz+1);
ExN = zeros(Nx  , Ny+1, Nz+1);
EyN = zeros(Nx+1, Ny  , Nz+1);
EzN = zeros(Nx+1, Ny+1, Nz  );
HxN = zeros(Nx+1, Ny  , Nz  );
HyN = zeros(Nx  , Ny+1, Nz  );
HzN = zeros(Nx  , Ny  , Nz+1);
%------------------------------------------------------------------
% Position of the excitation source
%------------------------------------------------------------------
XCT = -0.25;     %  source center 
YCT = 0;         %  source center
ZCT = 0;         %  source center
[dummy, k_e] = min(abs(x-XCT));  %   index of the excitation node  (exp. with Nx/2)  
[dummy, m_e] = min(abs(y-YCT));  %   index of the excitation node  (exp. with Nx/2)  
[dummy, p_e] = min(abs(z-ZCT));  %   index of the excitation node  (exp. with Nx/2)  
%------------------------------------------------------------------
% Position of the receiver source
%------------------------------------------------------------------
XCR = +0.25;
YCR = 0; 
ZCR = 0;   %   coordinates
[dummy, k_e1] = min(abs(x-XCR));  %   index of the excitation node  (exp. with Nx/2)  
[dummy, m_e1] = min(abs(y-YCR));  %   index of the excitation node  (exp. with Nx/2)  
[dummy, p_e1] = min(abs(z-ZCR));  %   index of the excitation node  (exp. with Nx/2) 
%------------------------------------------------------------------
%------------------------------------------------------------------
%   Antenna metal boundaries
%------------------------------------------------------------------
N = 0;          %   length of one wing (the length is N*d; the dipole length is (2*N+1)*d)
%   TX
MetalRectX      = k_e*ones(1,2*N);                  %  x-position of metal rectangles
MetalRectZ      = [p_e-N:p_e-1 p_e+1:p_e+N];        %  z-position of metal rectangles
XVR             = []; 
YVR             = []; 
for m = 1:length(MetalRectX)
    xmin_        = xmin + (MetalRectX(m)-1.5)*d; 
    xmax_        = xmin_ + d;
    zmin_        = zmin + (MetalRectZ(m)-1.5)*d; 
    zmax_        = zmin_ + d;
    XVR(:, m)   = [xmin_ xmin_ xmax_ xmax_ xmin_];
    ZVR(:, m)   = [zmin_ zmax_ zmax_ zmin_ zmin_];
end
%  RX
MetalRectX1      = k_e1*ones(1,2*N);              %  x-positon of metal rectangles
MetalRectZ1      = [p_e1-N:p_e1-1 p_e1+1:p_e+N];  %  z-position of metal rectangles
for m = 1:length(MetalRectX1)
    xmin_        = xmin + (MetalRectX1(m)-1.5)*d; 
    xmax_        = xmin_ + d;
    zmin_        = zmin + (MetalRectZ1(m)-1.5)*d; 
    zmax_        = zmin_ + d;
    XVR(:, m+length(MetalRectX))   = [xmin_ xmin_ xmax_ xmax_ xmin_];
    ZVR(:, m+length(MetalRectX))   = [zmin_ zmax_ zmax_ zmin_ zmin_];
end
%------------------------------------------------------------------
%------------------------------------------------------------------
%   Difference coefficients/dielectric boundaries - Matrix DIEL
%------------------------------------------------------------------
sigmae = SIGMA;
sigmah = 0;
eps    = eps0*EPS;
%------------------------------------------------------------------
e1     = (1 - dt*sigmae./(2*eps))./(1 + dt*sigmae./(2*eps));
e2     = (dt./(d*eps))./(1 + dt*sigmae./(2*eps));
e3     = (dt./(d*eps))./(1 + dt*sigmae./(2*eps));
h1      = (1 - dt*sigmah/(2*mu0))./(1 + dt*sigmah/(2*mu0));
h2      = (dt/(d*mu0))./(1 + dt*sigmah/(2*mu0));

sigma   = d/(d*d*RG); 
es1     = (1 - dt*sigma/(2*eps0))/(1 + dt*sigma/(2*eps0));
es2     = (dt/(d*eps0))/(1 + dt*sigma/(2*eps0));
es3     = (dt*sigma/(d*eps0))/(1 + dt*sigma/(2*eps0));

%------------------------------------------------------------------
%   Main loop - "bootstrapping" (initial conditions are zeros)
%------------------------------------------------------------------
n = 2; hr = []; hrd = []; VLoad = zeros(size(t));
tic
while n < NT+1   
    %------------------------------------------------------------------
    % E-update (everywhere except on boundary; (40% of time))
    ExN(:,2:Ny,2:Nz) = e1(1:Nx,2:Ny,2:Nz).*ExP(:,2:Ny,2:Nz) + e2(1:Nx,2:Ny,2:Nz).*(diff(HzP(:,:,2:Nz),1,2) - diff(HyP(:,2:Ny,:),1,3));
    EyN(2:Nx,:,2:Nz) = e1(2:Nx,1:Ny,2:Nz).*EyP(2:Nx,:,2:Nz) + e2(2:Nx,1:Ny,2:Nz).*(diff(HxP(2:Nx,:,:),1,3) - diff(HzP(:,:,2:Nz),1,1));
    EzN(2:Nx,2:Ny,:) = e1(2:Nx,2:Ny,1:Nz).*EzP(2:Nx,2:Ny,:) + e2(2:Nx,2:Ny,1:Nz).*(diff(HyP(:,2:Ny,:),1,1) - diff(HxP(2:Nx,:,:),1,2));       
    %------------------------------------------------------------------
    %------------------------------------------------------------------
    %   Boundary conditions
    abs_mur1; 
    %------------------------------------------------------------------
    %------------------------------------------------------------------
    %   Feed model - TX -lumped port 
    EzN(k_e, m_e,p_e) =      es1 * EzP(k_e, m_e,p_e)+ ...
                             es2 *(HyN(k_e, m_e, p_e) - HyN(k_e-1, m_e, p_e) - ...
                                   HxN(k_e, m_e, p_e) + HxN(k_e, m_e-1, p_e))- ...
                             es3 *(VG(n) + VG(n-1))/2;

    
    IG(n)         =   (d/RG) * EzN(k_e, m_e, p_e)  + VG(n)/RG ; 
                           
    EzN(MetalRectX, m_e, MetalRectZ)     =  0;           %   top/bottom wings
    %------------------------------------------------------------------
    %------------------------------------------------------------------
    %   Feed model - RX -lumped load 
    EzN(k_e1, m_e1, p_e1) =       es1 * EzP(k_e1, m_e1, p_e1)+ ...
                                  es2 *(HyN(k_e1, m_e1, p_e1) - HyN(k_e1-1, m_e1, p_e1) - ...
                                        HxN(k_e1, m_e1, p_e1) + HxN(k_e1, m_e1-1, p_e1));
    IG(n)         =   (d/RG) * EzN(k_e1, m_e1, p_e1); 
    VLoad(n)      = RG*IG(n);
  
    EzN(MetalRectX1, m_e, MetalRectZ)     =  0;           %   top/bottom wings
    %------------------------------------------------------------------
    %------------------------------------------------------------------
    %  H-update (40% of time)
    HxN =  h1*HxP + h2*(diff(EyN,1,3)- diff(EzN,1,2));   
    HyN =  h1*HyP + h2*(diff(EzN,1,1)- diff(ExN,1,3));      
    HzN =  h1*HzP + h2*(diff(ExN,1,2)- diff(EyN,1,1));      
    %------------------------------------------------------------------
    %------------------------------------------------------------------
    %   Boundary conditions
    abs_super;
    %------------------------------------------------------------------
    Energy(n) = TotalEnergy3D(ExN, EyN, EzN, HxN, HyN, HzN, d, eps0, mu0); 
    %------------------------------------------------------------------
    %   Next step
    n   = n + 1;      
    ExP = ExN; EyP = EyN; EzP = EzN; HxP = HxN; HyP = HyN; HzP = HzN;   
    %------------------------------------------------------------------ 
    %   Scale/Plot fields           
    output              = squeeze(EzP(:, m_e, :));    
    output              = abs(output).^0.25.*sign(output); 
    output(1,1)         =   -bound(n); 
    output(end,end)     =   +bound(n);     
    
    if (~isempty(hr)) delete(hr); end;
    if (~isempty(hrd)) delete(hrd); end;
    
    subplot(2,2,[1 3])
        imagesc( [x(1) x(end)], [z(1) z(end-1)], output');          %   main field
        %hr  = patch(XVR, ZVR, [0.5 0.5 0.5]);        
        string      =   strcat(num2str(1e9*t(n)), ' ns');  
        axis 'equal';    axis 'tight', set(gca,'YDir','normal'); 
        xlabel('x, m'); ylabel('z, m');
        title(strcat('FDTD Simulation: Electric field at t=',string), 'FontSize', 13);
    subplot(2,2,2)
        plot(t*1e9, VG,'b', t(1:n)*1e9, VLoad(1:n)*1e3, 'r'); grid on;
        title ('Received voltage (red, mV) vs. TX voltage (blue, V)', 'FontSize', 13); 
        xlabel('time, ns');
    subplot(2,2,4)
        plot(t(2:n-1)*1e9, Energy(2:n-1) ); grid on;
        title ('Total energy, linear', 'FontSize', 13); 
        xlabel('time, ns');
    drawnow;     
    %------------------------------------------------------------------
end
Energy(end)