%constants
c_0 = 3E8;                  % Speed of light in free space
mu_0 = 4.0*pi*1.0e-7;        % Permeability of free space
eps_0 = 1.0/(c_0*c_0*mu_0);  % Permittivty of free space

%to be set
Nx = 100;                    % Number of cells in x-direction
Ny = 100;                    % Number of cells in y-direction
Nt = 500;                    % Number of time steps

c_ref = c_0;                % Reference velocity
eps_ref = eps_0;
mu_ref = mu_0;

f_0 = 10e9;                  % Excitation frequency
f_ref = f_0;                % Reference frequency 

omega_0 = 2.0*pi*f_0;        % Excitation circular frequency
omega_ref = omega_0;

lambda_ref = c_ref/f_ref;    % Reference wavelength

Dx_ref = lambda_ref/20;      % Reference cells width
Dy_ref = lambda_ref/20;

X = Nx*Dx_ref;
Y = Ny*Dy_ref;

r = 0.5;                    % Normalized factor
Dt_ref = r/c_ref*Dx_ref;    % Reference time step
Dt = Dt_ref;

% initialization
Ez0 = zeros(Nx, Ny);
Ez1 = zeros(Nx, Ny);
Hx0 = zeros(Nx, Ny);
Hx1 = zeros(Nx, Ny);
Hy0 = zeros(Nx, Ny);
Hy1 = zeros(Nx, Ny);

% Source position
Nx_Source = int16(0.5*(Nx+1));
Ny_Source = int16(0.5*(Nx+1));

pulse = 0;

for n = 1:Nt            % Sources' definitions
    t = Dt_ref*r*n;    % Actual time
    Source_type = 1;    % Choice of source type
    
    switch Source_type 
        case 1          % Modified source function
            ncycles = 1;
            if t < ncycles*2.0*pi/(omega_0)  
                pulse = -0.5*( 1.0 - cos(omega_0*t/ncycles) ) * sin(omega_0*t);
            else 
                pulse = 0;
            end
        case 2          % Sigle cos source function
            if t < 2.0*pi/(omega_0)
                pulse = 8*c_0^2*Dt_ref^2*mu_ref*omega_0*cos(omega_0*t);
            else 
                pulse = 0;
            end
        case 3          % Gaussian pulse
            if t < Dt_ref*r*50
                pulse = -40*c_0*(t-t*25/(n-1))*exp(-(t-t*25/(n-1))^2/2/(50/2.3548)^2/(t/(n-1))^2);
            else
                pulse = 0;
            end
        otherwise      % For debug
            pulse = 1;
    end
    Ez0(Nx_Source,Ny_Source) = Ez0(Nx_Source,Ny_Source) - r*pulse;
    
    CHy = Dt_ref/mu_ref/Dx_ref;        % Coefficients used below
    CHx = Dt_ref/mu_ref/Dy_ref;
    CEzHy = Dt_ref/eps_ref/Dx_ref;
    CEzHx = Dt_ref/eps_ref/Dy_ref;
    
    for i = 2:Nx
            % H update
            Hx1(i,1:Ny-1) = Hx0(i,1:Ny-1) - CHx.*(Ez0(i,2:Ny)-Ez0(i,1:Ny-1));
    end
    for i = 1:Nx-1
            Hy1(i,2:Ny) = Hy0(i,2:Ny) + CHy.*(Ez0(i+1,2:Ny)-Ez0(i,2:Ny));
    end
    
    % Boundary conditions **************************************
    boundary = 3;      % Choice: '1'=Mur ABC; '2'=Dirichlet; '3'=Neumann
    switch boundary
        case 1  % For H Mur ABC
            Hx1(1,1:Ny-1) = Hx0(2, 1:Ny-1) + (r-1)/(r+1).*(Hx1(2, 1:Ny-1)-Hx0(1,1:Ny-1));    % Mur ABC @ left side x = 0
            Hx1(1:Nx, Ny) = Hx0(1:Nx,Ny-1) + (r-1)/(r+1).*(Hx1(1:Nx,Ny-1)-Hx0(1:Nx, Ny));    % Mur ABC @ left side y = Ny
            Hy1(1:Nx-1,1) = Hy0(1:Nx-1, 2) + (r-1)/(r+1).*(Hy1(1:Nx-1, 2)-Hy0(1:Nx-1,1));    % Mur ABC @ left side y = 0
            Hy1(Nx, 1:Ny) = Hy0(Nx-1,1:Ny) + (r-1)/(r+1).*(Hy1(Nx-1,1:Ny)-Hy0(Nx, 1:Ny));    % Mur ABC @ left side x = Nx
        case 2  % Dirichlet
            
            Hx1(1:Nx, 1) = 0;
            Hx1(1:Nx,Ny) = 0;
            Hx1(1, 1:Ny) = 0;
            Hx1(Nx,1:Ny) = 0;
            
            Hy1(1:Nx, 1) = 0;
            Hy1(1:Nx,Ny) = 0;
            Hy1(1, 1:Ny) = 0;
            Hy1(Nx,1:Ny) = 0;
        case 3  % Neumann
            Hx1(1, 1:Ny-1) = Hx0(1, 1:Ny-1);
            Hx1(Nx,1:Ny-1) = Hx1(Nx,1:Ny-1);
            Hx1(1:Nx,  Ny) = Hx0(1:Nx,  Ny);
            Hx1(1:Nx,  1) = Hx1(1:Nx,  1);
            
            Hy1(1:Nx-1, 1) = Hy0(1:Nx-1, 1);
            Hy1(1:Nx-1,Ny) = Hy0(1:Nx-1,Ny);
            Hy1(Nx,  1:Ny) = Hy0(Nx,  1:Ny);
            Hy1(1,  1:Ny) = Hy0(1,  1:Ny);
            
    end
    for i = 2:Nx
            % E update
            Ez1(i,2:Ny) = Ez0(i,2:Ny) + CEzHy.*(Hy1(i,2:Ny)-Hy1(i-1,2:Ny)) - CEzHx.*(Hx1(i,2:Ny)-Hx1(i,1:Ny-1));
    end
    switch boundary
        case 1  % For E Mur ABC
            Ez1(1,1:Ny) = Ez0(2,1:Ny) + (r-1)/(r+1).*(Ez1(2,1:Ny)-Ez0(1,1:Ny));            % Mur ABC @ right side x = 0
            Ez1(1:Nx,1) = Ez0(1:Nx,2) + (r-1)/(r+1).*(Ez1(1:Nx,2)-Ez0(1:Nx,1));            % Mur ABC @ right side y = 0
        case 2  % Dirichlet
            Ez1(1, 1:Ny) = 0;
            Ez1(1:Nx, 1) = 0;
            Ez1(Nx,1:Ny) = 0;
            Ez1(1:Nx,Ny) = 0;
        case 3  % Neumann
            Ez1(1, 1:Ny) = Ez0(1, 1:Ny);
            Ez1(Nx,1:Ny) = Ez0(Nx,1:Ny);
            Ez1(1:Nx, 1) = Ez0(1:Nx, 1);
            Ez1(1:Nx,Ny) = Ez0(1:Nx,Ny);
    end
            Hx0 = Hx1;
            Hy0 = Hy1; 
            Ez0 = Ez1;
    
    % Display*********************************************
    i = 1:Nx;
    j = 1:Ny;
    display = 1;        % Choice: '1' = Ez, '2' = Hx, '3' = Hy, 'Otherwise' = three components
    switch display
        case 1
            surf(i,j,Ez0);
            axis([0 Nx 0 Ny -0.03 0.03]);
            set(gca, 'XTick',[1 Nx/4 Nx/2 3*Nx/4 Nx],'FontSize',8);
            set(gca, 'XTickLabel',[0 X/4 X/2 3*X/4 X],'FontSize',8);
            xlabel('x in m');  
            set(gca, 'YTick',[1 Ny/4 Ny/2 3*Ny/4 Ny],'FontSize',8);
            set(gca, 'YTickLabel',[0 Y/4 Y/2 3*Y/4 Y],'FontSize',8);
            ylabel('y in m');
            zlabel('Amplitude of Ez');
        case 2
            surf(i,j,Hx0);
            axis([0 Nx 0 Ny -1e-4 1e-4]);
            set(gca, 'XTick',[1 Nx/4 Nx/2 3*Nx/4 Nx],'FontSize',8);
            set(gca, 'XTickLabel',[0 X/4 X/2 3*X/4 X],'FontSize',8);
            xlabel('x in m');  
            set(gca, 'YTick',[1 Ny/4 Ny/2 3*Ny/4 Ny],'FontSize',8);
            set(gca, 'YTickLabel',[0 Y/4 Y/2 3*Y/4 Y],'FontSize',8);
            ylabel('y in m');
            zlabel('Amplitude of Hx');
        case 3
            surf(i,j,Hy0);
            axis([0 Nx 0 Ny -1e-4 1e-4]);
            set(gca, 'XTick',[1 Nx/4 Nx/2 3*Nx/4 Nx],'FontSize',8);
            set(gca, 'XTickLabel',[0 X/4 X/2 3*X/4 X],'FontSize',8);
            xlabel('x in m');  
            set(gca, 'YTick',[1 Ny/4 Ny/2 3*Ny/4 Ny],'FontSize',8);
            set(gca, 'YTickLabel',[0 Y/4 Y/2 3*Y/4 Y],'FontSize',8);
            ylabel('y in m');
            zlabel('Amplitude of Hy');
        otherwise 
            subplot(2,2,1);
            surf(i,j,Ez0);
            title('Ex');
            axis([0 Nx 0 Ny -0.03 0.03]);
            set(gca, 'XTick',[1 Nx/4 Nx/2 3*Nx/4 Nx],'FontSize',8);
            set(gca, 'XTickLabel',[0 X/4 X/2 3*X/4 X],'FontSize',8);
            xlabel('x in m');  
            set(gca, 'YTick',[1 Ny/4 Ny/2 3*Ny/4 Ny],'FontSize',8);
            set(gca, 'YTickLabel',[0 Y/4 Y/2 3*Y/4 Y],'FontSize',8);
            ylabel('y in m');
            zlabel('Amplitude of Ez');
            
            subplot(2,2,2);
            surf(i,j,Hx0);
            title('Hx');
            axis([0 Nx 0 Ny -1e-4 1e-4]);
            set(gca, 'XTick',[1 Nx/4 Nx/2 3*Nx/4 Nx],'FontSize',8);
            set(gca, 'XTickLabel',[0 X/4 X/2 3*X/4 X],'FontSize',8);
            xlabel('x in m');  
            set(gca, 'YTick',[1 Ny/4 Ny/2 3*Ny/4 Ny],'FontSize',8);
            set(gca, 'YTickLabel',[0 Y/4 Y/2 3*Y/4 Y],'FontSize',8);
            ylabel('y in m');
            zlabel('Amplitude of Hx');
            
            subplot(2,2,3);
            surf(i,j,Hy0);
            title('Hy');
            axis([0 Nx 0 Ny -1e-4 1e-4]);
            set(gca, 'XTick',[1 Nx/4 Nx/2 3*Nx/4 Nx],'FontSize',8);
            set(gca, 'XTickLabel',[0 X/4 X/2 3*X/4 X],'FontSize',8);
            xlabel('x in m');  
            set(gca, 'YTick',[1 Ny/4 Ny/2 3*Ny/4 Ny],'FontSize',8);
            set(gca, 'YTickLabel',[0 Y/4 Y/2 3*Y/4 Y],'FontSize',8);
            ylabel('y in m');
            zlabel('Amplitude of Hy');
    end

    %surf(i,j,Ez0);
    %axis([0 Nx 0 Ny -0.03 0.03])
    %A = rot90(Ez0);
    %imagesc(A);
    pause(0.005);
end