%% Finite-Difference in Time Domain Method
% This is a minimum working example of an implementation of competition topic.
clc, clear, close all;

%% Constants and Properties
% Constants
c0 = 299792458; % speed of light
mu0 = 4*pi*1e-7; % vacuum permeability
eps0 = 1 / (c0^2*mu0); % vacuum permitivity

% Space dimensions and discretization
Lx = 1; % length of x-dimension
Ly = 1; % length of y-dimension
Nx = 101; % number of discretization points in x
Ny = 101; % number of discretization points in y
dx = Lx / (Nx - 1); % delta-x
dy = Ly / (Ny - 1); % delta-y
x = linspace(0, Lx, Nx); % space discretization in x
y = linspace(0, Ly, Ny); % space discretization in y

% Time discretization
Nt = 500; % number of discretization points in t
dt = 1 / c0 / sqrt(dx^-2 + dy^-2); % delta-t

% Excitation definition
f = 1 / (16*dt*pi); % excitation frequency
G = sin(2*pi*f*dt*(0:Nt-1)); % excitation signal
xG = floor(Nx / 2) + 1; % index of position of feeding in x
yG = floor(Nx / 2) + 1; % index of position of feeding in y

%% Wave equation solution
% Allocation of variables
g = zeros(Nx, Ny);
u = zeros(Nx, Ny, 3);

% Run solver
for k = 1:Nt
   % Excitation update
   g(xG,yG) = G(k);
   
   % Update in time
   u(:,:,3) = u(:,:,2);
   u(:,:,2) = u(:,:,1);
   
   % Update in space
   for i = 2:Nx-1
      for j = 2:Nx-1
         u(i,j,1) = dt^2*g(i,j) + 2*u(i,j,2) - u(i,j,3) + ...
         (dt*c0/dx)^2*(u(i-1,j,2) - 2*u(i,j,2) + u(i+1,j,2)) + ...
         (dt*c0/dy)^2*(u(i,j-1,2) - 2*u(i,j,2) + u(i,j+1,2));
      end 
   end
   
   % Boundary condition
   for i = 2:Nx-1
      u(i,[1 end],1) = u(i,[3 end-2],1);
      u([1 end],i,1) = u([3 end-2],i,1);
   end
   
   % Draw results
   surf(abs(u(:,:,1).'));
   xlim([2 Nx - 1]);
   ylim([2 Ny - 1]);
   view(0, 90);
   drawnow;
end