%//************************************************************************
%//     1D FDTD code with absorbing boundary conditions
%//***********************************************************************
%//    Program author: Dr.-Ing. René Marklein
%//                    University of Kassel, Department of Electrical and 
%//                    Computer Engineering, Electromagnetic Field Theory
%//                    Wilhelmshöher Allee 71
%//                    D-34121 Kassel
%//                    Tel. 0561 804 6426, marklein@uni-kassel.de
%//
%//    Date of this version:  November 15, 2001
%//
%//      This MATLAB program implements the finite-difference time-domain 
%//   (FDTD) solution of Maxwell's eqautions in differential form over a 
%//    one-dimensional space comprised of a staggered 1D grid of uniform
%//    grid cells.
%//      To illustrate the algorithm, a broadband raised cosine (RC) pulse 
%//    2 cycles and a carrier frequency of f_0 = 10 GHz propagating in vaccum 
%//    (epsilon_r=1.0, mu_r=1.0) is modeled.  
%//      The cell size of Dx = lambda_0/20, G = 20, is chosen to provide 20 
%//    grid cells per wavelength at center frequency. The normalzied time 
%//    step (Courant factor) normDt = c_0*Dt/Dx is set to normDt = 0.5.
%//      The computational domain is truncated using an absorbing boundary 
%//    condition.
%//
%//***********************************************************************

close all;
clear all;

fprintf('\n');
fprintf('+------------------------------------------------------+\n');
fprintf('| 1D FDTD code with absorbing boundary condition       |\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

%// Material parameters    
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;

%// Excitation pulse
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);

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

%// 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.5*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;

figure(1);

fp_ex_bin = fopen('Ex.bin','wb');
fp_hy_bin = fopen('Hy.bin','wb');
fp_pulse_bin = fopen('Pulse.bin','wb');

%// Initialization of help variables for the open boundary condition
Ex_low_m1 = 0;   Ex_low_m2 = 0;
Ex_high_m1 = 0;  Ex_high_m2 = 0;

%// Marching-on-in-time loop
fprintf('Start marching-on-in-time loop ...\n');
for nt=1:Nt         
    
    t = Dt_ref*normDt*(nt-1); %// Actual time
    fprintf('Time step: n_t = %5d; Time = %20.10g s\n', nt, t)

    %// 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) - normDt * ( 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
    t = Dt_ref*normDt*(nt-1); %// Actual time
    ncycles = 2;
    if t < ncycles*2.0*M_PI/(omega_0)  
        pulse(nt) = 0.5*( 1.0 - cos(omega_0*t/ncycles) ) * 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 * 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 
    
    nwritten = fwrite(fp_ex_bin,Ex,'float32');
    %// fprintf('%d byte written\n', nwritten);
    nwritten = fwrite(fp_hy_bin,Hy,'float32');
    %// fprintf('%d byte written\n', nwritten);
    nwritten = fwrite(fp_pulse_bin,pulse(nt),'float32');
    %// fprintf('%d byte written\n', nwritten);
    
    %// Generate output figures
    
    subplot(2,1,1)
    plot(Ex,'b')
    t = Dt_ref*normDt*(nt-1+0.5); %// Actual time
    title(['E_x component at t = ',num2str(t*1e12),' ps']);
    axis([1 Nz -1 1]);   
    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.0 -0.5 0 0.5 1.0],'FontSize',8);
    set(gca, 'YTickLabel',[-1 -0.5 0 0.5 1.0],'FontSize',8);
    ylabel('Amplitude of E_x');
 
    grid on 

    subplot(2,1,2)
    plot(Hy,'r')
    t = Dt_ref*normDt*(nt-1); %// Actual time
    title(['H_y component at t = ',num2str(t*1e12),' ps']);
    axis([1 Nz -1.0 1.0]);   
    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.0 -0.5 0 0.5 1.0],'FontSize',8);
    set(gca, 'YTickLabel',[-1 -0.5 0 0.5 1.0],'FontSize',8);
    ylabel('Amplitude of H_y');
    grid on   
    
    pause(0.005)
   
end %// End of marching-on-in-time loop

fprintf('Marching-on-in-time finished\n');

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

fprintf('Generate an EPS output file of Figure 1\n');
print -depsc2 -r300 fdtd1dvacab.eps %// Generate a EPS output file of figure 1
