%/************************************************************************
%     1D FDTD code with PEC 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 
%   two cycles and a carrier frequency of f_0 = 700 MHz propagating in 
%   a 1-D layered medium 
%   (epsilon_r = 1.0, mu_r = 1.0, sigma_e = 0.0) 
%   (epsilon_r = 4.0, mu_r = 1.0, sigma_e = 0.04) 
%   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 the perfectly 
%   electrically conducting (PEC) boundary condition.
%//***********************************************************************

close all; 
clear all;

fp_proto  = fopen('fdtd1dpercon.txt','w');

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

debug = 0;                   %// Local debugger flag: 0=off, 1=on

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 cell number

E_x   = zeros(Nz,1);             %/* (normalized) Ex component  */
D_x   = zeros(Nz,1);             %/* (normalized) Dx component  */
I     = zeros(Nz,1);             %/* I help field               */    
H_y   = zeros(Nz,1);             %/* (normalized) Hy component  */
J_ex  = zeros(Nt,1);             %/* (normalized) Jex component */
S_z   = zeros(Nz,1);             %/* (normalized) Sz component  */
pulse = zeros(Nt,1);             %/* Pulse                      */

epsilon = zeros(Nz,1);
sigma_e = 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  */
epsilon_0 = 1.0/(c_0*c_0*mu_0);  %/* Permittivity of free space  */

c_ref = c_0;                     %/* Reference (phase) velocity  */
eps_ref = epsilon_0;             %/* Reference permittivity      */ 
mu_ref = mu_0;                   %/* Reference permeability      */ 

f_0 = 0.7e9;                     %/* Excitation frequency        */
fprintf('Excitation carrier frequency   :          f_0 = %20.10g 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 = %20.10g 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 = %20.10g m\n', Dx_ref);

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

sigma_e_ref = eps_ref / Dt_ref;

epsilon_r_1 = 1.0;
epsilon_1   = epsilon_r_1 * epsilon_0;
mu_r_1      = 1.0;
mu_1        = mu_r_1 * mu_0;
sigma_e_1   = 0.0;

Z_1 = sqrt( mu_1 / ( epsilon_0 * epsilon_r_1 + j * sigma_e_1 / ( omega_0 * epsilon_0 ) ) );

fprintf('Wave impedance                :       Z_1 = %g + j %g V/A\n', real(Z_1), imag(Z_1));

epsilon_r_2 = 4.0;
epsilon_2   = epsilon_r_2 * epsilon_0;
mu_r_2      = 1.0;
mu_2        = mu_r_2 * mu_0;
sigma_e_2   = 0.002;

Z_2 = sqrt( mu_2 / ( epsilon_0 * ( epsilon_r_2 + i * sigma_e_2 / ( omega_0 * epsilon_0 ) ) ) );

fprintf('Wave impedance              :          Z_2 = %g + j %g V/A\n', real(Z_2) , imag(Z_2));

R_12_E = (Z_2-Z_1)/(Z_2+Z_1);
T_12_E = 2*Z_2/(Z_2+Z_1);

R_12_H = -(Z_2-Z_1)/(Z_2+Z_1);
T_12_H = 2*Z_1/(Z_2+Z_1);

fprintf('Reflection coefficient   E_x: Z_1 -> Z_2  R_E = %g + j %g\n', real(R_12_E), imag(R_12_E) );
fprintf('Transmission coefficient E_x: Z_1 -> Z_2  T_E = %g + j %g\n', real(T_12_E), imag(T_12_E) );

fprintf('Reflection coefficient   H_y: Z_1 -> Z_2  R_H = %g + j %g\n', real(R_12_H), imag(R_12_H) );
fprintf('Transmission coefficient H_y: Z_1 -> Z_2  T_H = %g + j %g\n', real(T_12_H), imag(T_12_H) );

alpha_1 = omega_0 / c_0 * sqrt( epsilon_r_1/ 2 ) * sqrt( sqrt ( 1 + (sigma_e_1/(omega_0*epsilon_0*epsilon_r_1))^2 - 1 ) );
beta_1  = omega_0 / c_0 * sqrt( epsilon_r_1/ 2 ) * sqrt( sqrt ( 1 + (sigma_e_1/(omega_0*epsilon_0*epsilon_r_1))^2 + 1 ) );

1/alpha_1;

alpha_2 = omega_0 / c_0 * sqrt( epsilon_r_2/ 2 ) * sqrt( sqrt ( 1 + (sigma_e_2/(omega_0*epsilon_0*epsilon_r_2))^2 - 1 ) );
beta_2  = omega_0 / c_0 * sqrt( epsilon_r_2/ 2 ) * sqrt( sqrt ( 1 + (sigma_e_2/(omega_0*epsilon_0*epsilon_r_2))^2 + 1 ) ); 

1/alpha_2;

skin_2 = 1/sqrt( mu_0 * mu_r_2 * sigma_e_2 * pi *f_0 );

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

normDt = 1.0;                    %// Normalized time step
fprintf('Normalized delta t             :    tilde(Dt) = %20.10g\n', normDt);

T = Nt*Dt_ref*normDt;
fprintf('Simulated time interval [0-T]  :            T = %20.10g\n', T);

%// Initialize field quantities     
for n=1:N 
    E_x(n) = 0.0;
    D_x(n) = 0.0;
    H_y(n) = 0.0;
    I(n)   = 0.0;
end

%// Initialize field quantities     
for n=1:Nz 
    epsilon(n) = 1.0;
    sigma_e(n) = 0.0;
end

for n=0.5*Nz:Nz 
    epsilon(n) = 4.0;
    sigma_e(n) = 0.04 / sigma_e_ref;
end

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

%// Simulation time
fprintf('Simulation time                :            T = %20.10g s\n', normDt*Dt_ref*Nt);
        
%// Source position at the center
Nz_Source = 0.25*Nz;
Nz_Source = 1;
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;

f1 = figure('Units','centimeters','Position',[2 0 16 21])
colormap(jet(256))

fp_ex_bin = fopen('Ex.bin','wb');  %// Open all output files
fp_hy_bin = fopen('Hy.bin','wb');
fp_sz_bin = fopen('Sz.bin','wb');
fp_pulse_bin = fopen('Pulse.bin','wb');

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

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

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

    %// Compute Ampere-Maxwell's grid equation          
    for n=2:N-1
        D_x(n) = D_x(n) - normDt * ( H_y(n+1) - H_y(n) );
    end
    
    for n=2:N-1
        E_x(n) = ( D_x(n) - I(n) ) / ( epsilon(n) + normDt * sigma_e(n) );
    end
    
    %// Apply PEC boundary condition at left and right boundary
    Ex(1) = 0; Ex(N) = 0;

    for n=2:N-1
        I(n) = I(n) + normDt * sigma_e(n) * E_x(n);
    end    
    I(1) = 0; I(N) = 0;
    
    %// Apply excitation pulse
    %// Raised cosine pulse with two cycles - (RC2) pulse
    t = Dt_ref*normDt*(nt-1); %// Actual time
    
    %// RCn pulse
    ncycle = 2;
    if ( t > 0 ) & ( t < ncycle*2.0*M_PI/(omega_0) )  
        pulse(nt) = 0.5 * ( 1.0 - cos(omega_0*t/ncycle) ) * cos(omega_0*t);
    else 
        pulse(nt) = 0;
    end
    
    % pulse(nt) = sin(omega_0*t);
               
    % /* Apply excittion pulse as an impress electric current density */
    J_ex(nt) = pulse(nt);

    % /* Apply excittion pulse as an impress electric current density */
    E_x(N_Source) = Ex(N_Source) - normDt * J_ex(nt);
        
    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 = %20.10g Ex_max = %20.10g\n', Ex_min, Ex_max);
    end
    
    Hy_max = max(max(Hy)); Hy_min = min(min(Hy));
    if debug == 1
        fprintf('Hy_min = %20.10g Hy_max = %20.10g\n', Hy_min, Hy_max);
    end 
    
    for n=1:N
        S_z(n) = E_x(n) * H_y(n);
    end
    
    nwritten = fwrite(fp_ex_bin,E_x,'float32');
    %// fprintf('%d byte written\n', nwritten);
    nwritten = fwrite(fp_hy_bin,H_y,'float32');
    %// fprintf('%d byte written\n', nwritten);
    nwritten = fwrite(fp_sz_bin,S_z,'float32');
    %// fprintf('%d byte written\n', nwritten);
    nwritten = fwrite(fp_pulse_bin,pulse(nt),'float32');
    %// fprintf('%d byte written\n', nwritten);
    
    %// Generate output figures
    
    set(gca,'FontSize',8);
    set(gca,'LineWidth',1);

    subplot('Position',[0.1 0.8 0.8 0.15])
    plot(J_ex,'b')
    t = Dt_ref*normDt*(nt-1+0.5); %// Actual time
    title(['J_{ex} component']);
    text(0.8*Nt,1.05,['t = ',num2str(t*1e12),' ps'],'FontSize',8);
    axis([1 Nt -2 2]);   
    set(gca, 'XTick',[1 Nt/4 Nt/2 3*Nt/4 Nt],'FontSize',8);
    set(gca, 'XTickLabel',[0 T/4 T/2 3*T/4 T],'FontSize',8);
    xlabel('t in ps');    
    set(gca, 'YTick',[-2 -1.0 0 1.0 2],'FontSize',8);
    set(gca, 'YTickLabel',[-2 -1.0 0 1.0 2],'FontSize',8);
    ylabel('Amplitude of J_ex');
    grid on 
    
    subplot('Position',[0.1 0.55 0.8 0.15])
    plot(E_x,'b');
    hold on
    X = [Nz/2 Nz/2];
    Y = [-1 1];
    plot(X,Y,'k');
    X = [(skin_2/Z*Nz+Nz/2) (skin_2/Z*Nz+Nz/2)];
    Y = [-1 1];
    plot(X,Y,'r');
    hold off
    t = Dt_ref*normDt*(nt-1+0.5); %// Actual time
    title(['E_x component']);
    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.5 0 0.5 1],'FontSize',8);
    set(gca, 'YTickLabel',[-1 -0.5 0 0.5 1],'FontSize',8);
    ylabel('Amplitude of E_x');
    
    grid on 

    subplot('Position',[0.1 0.3 0.8 0.15])
    plot(H_y,'r')
    t = Dt_ref*normDt*(nt-1); %// Actual time
    title(['H_y component']);
    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.5 0 0.5 1],'FontSize',8);
    set(gca, 'YTickLabel',[-1 -0.5 0 0.5 1],'FontSize',8);
    ylabel('Amplitude of H_y');
    grid on   

    subplot('Position',[0.1 0.05 0.8 0.15])
    plot(S_z,'g')
    t = Dt_ref*normDt*(nt-1); %// Actual time
    title(['S_z component']);
    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.5 0 0.5 1],'FontSize',8);
    set(gca, 'YTickLabel',[-1 -0.5 0 0.5 1],'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 stopped\n');

fclose(fp_ex_bin);     %// Close all output files
fclose(fp_hy_bin);
fclose(fp_sz_bin);
fclose(fp_pulse_bin);

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

fclose(fp_proto);