%----------  finite difference time domain method  -------------------
clear;
% ----------- define the domain of calculation ------------------
NumobjectX = 5;                                  % ---- the numbers of object on the X(Y) direction
NumobjectY = 5;
a = 0.2;                                          % ------ the lattice constant (unit: um )
Gridsize = 20;                                    % -------  the number of grids between two lattice point
DX = a / Gridsize;
DY = a / Gridsize;
domainLX = 0 : DX : (NumobjectX+1)*a;
domainLY = 0 : DY : (NumobjectX+1)*a;
% ------------  metiral params --------------------------
e0=4*pi*1.0e-1; %Epsilon Zero, if using Gauss Unit, it equals to 1.
mu0=8.85*1e-6;   %Mu Zero, if using Gauss Unit, it equals to 1.
c=1/sqrt(mu0*e0);
% ----------------  the dielectric objects ----------------------
%           ---------  the substrate ----------
R = 0.2 * Gridsize;
Eps_back = 1;
Eps_back = Eps_back * e0;
Eps_object = 12;
Eps_object = Eps_object * e0;
Eps_distribute = ones( length(domainLX),length(domainLY) ) * Eps_back;
eps_show = Eps_distribute;
for n = 1 : NumobjectX
    for m = 1 : NumobjectY
        for p = 1 : length(domainLX)
            for q = 1 : length(domainLY)
                if (p-(n*Gridsize+1))^2+(q-(m*Gridsize+1))^2 < R^2
                    Eps_distribute(p,q) = Eps_object;
                end
            end
        end
    end
end
%           ----------  the defect -----------
for p = 1: length(domainLX)+1
    for q = 1 : length(domainLX)+1
        if (p-(3*Gridsize+1))^2+(q-(3*Gridsize+1))^2 < R^2
            Eps_distribute(p,q) = Eps_back;
        elseif (p-(4*Gridsize+1))^2+(q-(3*Gridsize+1))^2 < R^2
            Eps_distribute(p,q) = Eps_back;
        elseif (p-(5*Gridsize+1))^2+(q-(3*Gridsize+1))^2 < R^2
            Eps_distribute(p,q) = Eps_back;
        end
    end
end
%   ------------------ show the origin ------------------------         
surf(domainLX,domainLY,Eps_distribute);
shading interp;
view(0,90);
axis([min(domainLX), max(domainLX),min(domainLY), max(domainLY)])
axis equal
disp('Press any key to continue...');
pause
%  -------------initial the field -----------------------
Hx = zeros(length(domainLX)+1,length(domainLY));
Hy = zeros(length(domainLX),length(domainLY)+1);
Ez = zeros(length(domainLX),length(domainLY));
% --------------the source ----------------------------
lambda = 532*1e-3;
omega=2*pi*c/lambda;
%---------------- the time params --------------------------
NTimeSteps=1000; %Total number of Time Steps
DT = 1/(c*sqrt(1/DX^2+1/DY^2));
for T=1:NTimeSteps
    Hx(2:length(domainLX),:) = Hx(2:length(domainLX),:) + ...
        DT*(Ez(1:length(domainLX)-1,:)-Ez(2:length(domainLX),:))/(mu0*DY);
   
    Hy(:,2:length(domainLY)) = Hy(:,2:length(domainLY)) + ...
        DT*(Ez(:,2:length(domainLY))-Ez(:,1:length(domainLY)-1))/(mu0*DX);
   
    Ez = Ez +DT*( (Hy(:,2:length(domainLY)+1)-Hy(:,1:length(domainLY)))/DX ...
        -( Hx(2:length(domainLX)+1,:)-Hx(1:length(domainLX),:))/DY)./Eps_distribute;
   
    Ez(round(length(domainLY)/2),round(length(domainLY)/2)) = ...
        Ez(round(length(domainLY)/2),round(length(domainLY)/2)) + sin(omega*DT*T);
   
    clf
   
    surf(Ez);
    Colormax=0.1;
    caxis([-Colormax Colormax]);
    shading interp;
   
    axis equal;
    axis off;
    view(0,90);
    colorbar
    pause(0.001)
end