%//************************************************************************
%//     1-D finite difference solution of the scalar wave equation 
%//     with perfectly electrically conducting (PEC) boundary conditions
%//     at both sides.
%//***********************************************************************
%//    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. ++49 0561 804 6426, marklein@uni-kassel.de
%//
%//    Date of this version:  November 12, 2002
%//
%//    This MATLAB program implements the finite difference (FD) solution of 
%//    the scalar wave equation over a one-dimensional space discretized by a 
%//    1D grid of uniform grid cells.
%//      To illustrate the algorithm a broadband raised cosine (RC) pulse with
%//    two 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 the carrier frequency. The normalzied time 
%//    step (Courant number) normDt = c_0*Dt/Dx is set to normDt = 1.0.
%//      The computational domain is truncated using the perfectly 
%//    electrically conduction (PEC) boundary condition.
%//***********************************************************************
close all; clear all;

fprintf('\n');
fprintf('+---------------------------------------------------+\n');
fprintf('| 1-D finite difference (FD) solution of the scalar |\n');
fprintf('| wave equation 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 = 400;                    %// Number of time steps
fprintf('Number of time steps           :           Nt = %d\n', Nt)
N = Nz;                      %// Global cell number

Ex      = zeros(Nz,1);           %// (normalized) Ex component
Ex_old  = zeros(Nz,1);           %// (normalized) Ex_old component
Ex_new  = zeros(Nz,1);           %// (normalized) Ex_old component
Hy      = zeros(Nz,1);           %// (normalized) Hy component
Jex     = zeros(Nt,1);           %// (normalized) Jex componen
Sz      = zeros(Nz,1);           %// (normalized) Sz componentt
pulse   = zeros(Nt,1);           %// Pulse
    
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
Z_0   = sqrt(mu_0/eps_0);    %// Wave impednace of free space

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

f_0 = 10e9;                  %// Excitation frequency
fprintf('Excitation 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;

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 
  Ex(n)     = 0.0;
  Ex_old(n) = 0.0;
  Ex_new(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 = %20.10g 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');  %// 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
  fprintf('Time step: n_t = %5d; Time = %20.10g s\n', nt, t);
    
  %// Compute wave FD equation          
  for n=2:N-1
    Ex_new(n) = 2.0*(1.0-normDt^2) * Ex(n) - Ex_old(n) + normDt^2 * ( Ex(n+1) + Ex(n-1) );
  end
  
  %// Apply PEC boundary condition at left and right boundary
  Ex_new(1) = 0;
  Ex_new(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
  
  %// Apply excittion pulse as an impress electric current density
  Jex(nt) = pulse(nt);
  
  %// Apply excittion pulse as an impress electric current density
  if nt > 1 
    Ex_new(N_Source) = Ex_new(N_Source) - normDt * ( Jex(nt) - Jex(nt-1) );
  else
    Ex_new(N_Source) = Ex_new(N_Source) - normDt * Jex(nt);
  end
  
  if debug == 1    
    fprintf('t= %8.4f pulse= %8.4f  Ex[N_Source]= %8.4f\n', t, pulse(nt), Ex_new(N_Source));
  end
  
  for n=1:N
    Ex_old(n) = Ex(n);
    Ex(n)     = Ex_new(n);
  end
  
  for n=2:N-1
    Hy(n) = Hy(n) - 0.5*normDt * ( Ex(n+1) - Ex(n-1) );
  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
     
  for n=1:N
    Sz(n) = Ex(n) * Hy(n);
  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_sz_bin,Sz,'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(4,1,1)
  plot(Jex,'b','LineWidth',1.5)
  t = Dt_ref*normDt*(nt-1+0.5); %// Actual time
  title(['J_{ex} component at t = ',num2str(t*1e12),' ps']);
  axis([1 Nt -1 1]);   
  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('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 J_ex');
  grid on 
  
  subplot(4,1,2)
  plot(Ex,'b','LineWidth',1.5)
  t = Dt_ref*normDt*(nt-1+0.5); %// Actual time
  title(['E_x component at t = ',num2str(t*1e12),' ps']);
  axis([1 Nz -0.5 0.5]);   
  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',[-0.5 -0.25 0 0.25 0.5],'FontSize',8);
  set(gca, 'YTickLabel',[-0.5 -0.25 0 0.25 0.5],'FontSize',8);
  ylabel('Amplitude of E_x');
  grid on 
  
  subplot(4,1,3)
  plot(Hy,'r','LineWidth',1.5)
  t = Dt_ref*normDt*(nt-1); %// Actual time
  title(['H_y component at t = ',num2str(t*1e12),' ps']);
  axis([1 Nz -0.5 0.5]);   
  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',[-0.5 -0.25 0 0.25 0.5],'FontSize',8);
  set(gca, 'YTickLabel',[-0.5 -0.25 0 0.25 0.5],'FontSize',8);
  ylabel('Amplitude of H_y');
  grid on   
  
  subplot(4,1,4)
  plot(Sz,'g','LineWidth',1.5)
  t = Dt_ref*normDt*(nt-1); %// Actual time
  title(['S_z component at t = ',num2str(t*1e12),' ps']);
  axis([1 Nz -0.25 0.25]);   
  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',[-0.25 -0.125 0 0.125 0.25],'FontSize',8);
  set(gca, 'YTickLabel',[-0.25 -0.125 0 0.125 0.25],'FontSize',8);
  ylabel('Amplitude of S_z');
  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 fd1dvac_rcn_poyn.eps %// Generate a EPS output file of figure 1