%//************************************************************************
%//    1D FDTD code modelling a pulse hitting a dielectric medium 
%//    with absorbing boundary conditions at the left and right side
%//***********************************************************************
%//
%//    Program author: Dr.-Ing. René Marklein
%//                    Department of Electrical and Computer Engineering
%//                    University of Kassel
%//                    Wilhelmshöher Allee 71
%//                    D-34121 Kassel
%//                    Tel. 0561 804 6426 
%//                    marklein@uni-kassel.de
%//
%//    Date of this version:  January 15, 2002
%//
%//    This MATLAB program implements the finite-difference time-domain (FDTD)
%//    solution of Maxwell's curl equations over a one-dimensional space
%//    comprised of a uniform staggered grid.
%//
%//    To illustrate the algorithm, a broadband raised cosine (RC) pulse 
%//    propagating hitting an interface between vaccum (epsilon_r=1.0, mu_r=1.0) 
%//    and a dielectric medium with a relative permittivity of epsilon_r=4.  
%//
%//    The cell size of Dx = lambda/20, G=20, is chosen to provide 20 grid 
%//    cells per wavelength. The normalzied time step (Courant factor)
%//    normDt= c*Dt/Dx is set to normDt=0.5.
%//
%//    The computational domain is truncated using a absorbing boundary 
%//    condition at the lower and higher boundary.
%//
%//***********************************************************************

close all;
clear all;

fprintf('\n');
fprintf('+-------------------------------------------------+\n');
fprintf('| 1D FDTD code with absorbing boundary conditions |\n');
fprintf('+-------------------------------------------------+\n\n');

debug = 0;

Nz = 200;                    %// Number of cells in z-direction
fprintf('Number of cells in z-direction :           Nz = %d\n', Nz)

Nt = 500;                    %// Number of time steps
fprintf('Number of time steps           :           Nt = %d\n', Nt)

N = Nz;                      %// Global number of cells

Ex         = zeros(Nz,1);      %// Ex component
Hy         = zeros(Nz,1);      %// Hy component
pulse      = zeros(Nt,1);    %// Pulse
rel_permittivity = zeros(Nz,1);      %// Relative permittivity 
beta             = zeros(Nz,1);
     
c_0 = 2.99792458E8;          %// Speed of light in free space

M_PI = pi;
mu_0 = 4.0*M_PI*1.0e-7;      %// Permeability of free space
eps_0 = 1.0/(c_0*c_0*mu_0);  %// Permittivty of free space

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

f_0 = 10e9;                  %// Excitation frequency
fprintf('Excitation frequency           :          f_0 = %g Hz\n', f_0)
f_ref = f_0;                 %// Reference frequency 

omega_0 = 2.0*M_PI*f_0;      %// Excitation circular frequency
omega_ref = omega_0;

lambda_ref = c_ref/f_ref;    %// Reference wavelength

fprintf('Reference wavelength           :   lambda_ref = %g m\n', lambda_ref)

Dx_ref = lambda_ref/20;          %// Reference cells width
Dz = Dx_ref;
nDz = Dz/Dx_ref;

Z = Nz*Dz;

fprintf('Reference cell size            :       Dx_ref = %g m\n', Dx_ref);

Dt_ref = Dx_ref/c_ref;       %// Reference time step
Dt = Dt_ref;

fprintf('Reference time step            :       Dt_ref = %g s\n', Dt_ref);

normDt = 0.5;                    %// Normalized time step

fprintf('Normalized delta t             :    tilde(Dt) = %g\n', normDt);

%// Initialize field quantities     
for n=1:N 
    Ex(n) = 0.0;
    Hy(n) = 0.0;
end

eps_r_1 = 1.0;
mu_r_1  = 1.0;

Z_1 = sqrt(mu_0*mu_r_1/(eps_0*eps_r_1));

eps_r_2 = 4.0;
mu_r_2  = 1.0;

Z_2 = sqrt(mu_0*mu_r_2/(eps_0*eps_r_2));

fprintf('Wave impedance                 :          Z_1 = %g\n', Z_1);
fprintf('Wave impedance                 :          Z_2 = %g\n', Z_2);

fprintf('Reflection coefficient   E_x: Z_1 -> Z_2  R_E = %g\n', (Z_2-Z_1)/(Z_2+Z_1));
fprintf('Transmission coefficient E_x: Z_1 -> Z_2  T_E = %g\n', 2*Z_2/(Z_2+Z_1));

fprintf('Reflection coefficient   H_y: Z_1 -> Z_2  R_H = %g\n', -(Z_2-Z_1)/(Z_2+Z_1));
fprintf('Transmission coefficient H_y: Z_1 -> Z_2  T_H = %g\n', 2*Z_1/(Z_2+Z_1));

%// Define the relative permittivity distribution
for n=1:N 
    rel_permittivity(n) = 1.0;
end
for n=N/2:N 
    rel_permittivity(n) = 4.0;
end

% Calculate the Define the 
for n=1:N-1 
    beta(n) = normDt*2/(rel_permittivity(n) + rel_permittivity(n+1));
end
beta(N) = normDt/rel_permittivity(N);


%// Number of time steps
fprintf('Number of time steps           :           Nt = %d\n', Nt);

%// Simulation time
fprintf('Simulation time                :            T = %g s\n', normDt*Dt_ref*Nt);
        
%// Source position at the center
Nz_Source = 0.25*Nz;
fprintf('Jex position [%d-%d]           :    Nz_Source = %d\n', 1, Nz, Nz_Source);
fprintf('                               :     z_Source = %d m\n', Dz*Nz_Source);

N_Source = Nz_Source;

%// Marching-on-in-time loop

figure(1);

fp_ex_bin = fopen('ex.bin','w');
fp_hy_bin = fopen('hy.bin','w');
fp_pulse_bin = fopen('pulse.bin','w');

Ex_low_m1 = 0;   Ex_low_m2 = 0;
Ex_high_m1 = 0;  Ex_high_m2 = 0;

for nt=1:Nt
    
    t = Dt_ref*normDt*(nt-1); %// Actual time
    if debug == 1
      fprintf('Time = %g s\n', t);
    end

    %// Compute Faraday's induction grid equation
    for n=2:N
        Hy(n) = Hy(n) - normDt * ( Ex(n) - Ex(n-1) ); 
    end

    %// Initialize all Hy components outside the simulation area to zero
    Hy(1) = 0.0;

    %// Compute Ampere-Maxwell's grid equation          
    for n=2:N-1
        Ex(n) = Ex(n) - beta(n) * ( Hy(n+1) - Hy(n) );
    end

    %// Apply PEC boundary condition at left and right boundary
    Ex(1) = Ex_low_m2;
    Ex_low_m2 = Ex_low_m1;
    Ex_low_m1 = Ex(2);
    
    Ex(N) = Ex_high_m2;
    Ex_high_m2 = Ex_high_m1;
    Ex_high_m1 = Ex(N-1);
    
    %// Apply excitation pulse
    %// Raised cosine pulse with two cycles - (RC2) pulse
    if t < 2.0*2.0*M_PI/(omega_0)  
        pulse(nt) = 0.5*( 1.0 - cos(omega_0*t/2.0) ) * cos(omega_0*t);
    else 
        pulse(nt) = 0;
    end
            
    %// Apply excittion pulse as an impress electric current density
    Jex = pulse(nt);

    %// Apply excittion pulse as an impress electric current density
    Ex(N_Source) = Ex(N_Source) - normDt * 2.0 * Jex;
        
    if debug == 1    
        fprintf('t= %8.4f pulse= %8.4f  Ex[N_Source]= %8.4f\n', t, pulse(nt), Ex(N_Source));
    end
    
    Ex_max = max(max(Ex));  Ex_min = min(min(Ex));
    if debug == 1
        fprintf('Ex_min = %g Ex_max = %g\n', Ex_min, Ex_max);
    end
    
    Hy_max = max(max(Hy)); Hy_min = min(min(Hy));
    if debug == 1
        fprintf('Hy_min = %g Hy_max = %g\n', Hy_min, Hy_max);
    end 
    
    fwrite(fp_ex_bin,Ex,'float32');
    fwrite(fp_hy_bin,Hy,'float32');
    fwrite(fp_pulse_bin,pulse(nt),'float32');
    
    %// Generate output figures
    
    subplot(2,1,1)
    plot(Ex,'b')
    title(['E_x component at t = ',num2str(t*1e12),' ps']);
    axis([1 Nz -1.33 1.33]);   
    set(gca, 'XTick',[1 Nz/4 Nz/2 3*Nz/4 Nz],'FontSize',8);
    set(gca, 'XTickLabel',[0 Z/4 Z/2 3*Z/4 Z],'FontSize',8);
    xlabel('z in m');    
    set(gca, 'YTick',[-1.33 -1.0 -0.5 0 0.5 1.0 1.33],'FontSize',8);
    set(gca, 'YTickLabel',[-1.33 -1 -0.5 0 0.5 1.0 1.33],'FontSize',8);
    ylabel('Amplitude of E_x');
 
    grid on 

    subplot(2,1,2)
    plot(Hy,'r')
    title(['H_y component at t = ',num2str(t*1e12),' ps']);
    axis([1 Nz -1.33 1.33]);   
    set(gca, 'XTick',[1 Nz/4 Nz/2 3*Nz/4 Nz],'FontSize',8);
    set(gca, 'XTickLabel',[0 Z/4 Z/2 3*Z/4 Z],'FontSize',8);
    xlabel('z in m');    
    set(gca, 'YTick',[-1.33 -1.0 -0.5 0 0.5 1.0 1.33],'FontSize',8);
    set(gca, 'YTickLabel',[-1.33 -1 -0.5 0 0.5 1.0 1.33],'FontSize',8);
    ylabel('Amplitude of H_y');
    grid on   
    
    pause(0.005)
   
end %// End of marching-on-in-time loop

fclose(fp_ex_bin);
fclose(fp_hy_bin);
fclose(fp_pulse_bin);

%// print -depsc2 -r300 fdtd1d_fig_1.eps