%//************************************************************************
%//     2-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
excitation = 2;

Nx = 200;                    %// Number of cells in z-direction
fprintf('Number of cells in x direction :           Nx = %d\n', Nx)
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 = Nx*Nz;                      %// Global cell number

Ey      = zeros(Nx,Nz);         %// (normalized) Ey component
Ey_old  = zeros(Nx,Nz);         %// (normalized) Ey_old component
Ey_new  = zeros(Nx,Nz);         %// (normalized) Ey_old component
Jey     = zeros(Nx,Nt);          %// (normalized) Jex componen
pulse   = zeros(Nt,1);          %// Pulse

Nrec=5;
fprintf('Number of receivers     :           Nrec = %d\n', Nrec)
Ey_rec  = zeros(Nt,Nrec);       %// (normalized) Ey component
    
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;                  %// Eycitation frequency
fprintf('Eycitation frequency           :          f_0 = %20.10g Hz\n', f_0)
f_ref = f_0;                 %// Reference frequency 

theta_in = 30;               %// Angle of incidence   

omega_0 = 2.0*M_PI*f_0;      %// Eycitation 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
Dx = Dx_ref;
Dz = Dx_ref;
nDz = Dz/Dx_ref;

X = Nx*Dz;
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 = 0.5;                    %// 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 nz=1:Nz 
  for nx=1:Nx 
    Ey(nx,nz)     = 0.0;
    Ey_old(nx,nz) = 0.0;
    Ey_new(nx,nz) = 0.0;
  end
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
Nx_Source = round(0.5*Nx);
Nz_Source = round(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);

figure(1);

fp_ex_bin = fopen('Ey.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');


nt_snapshot       = 40;
nt_snapshot_count = 0;

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 nz=2:Nz-1
      for nx=2:Nx-1
        Ey_new(nx,nz) = 2.0*(1.0-2.0*normDt^2) * Ey(nx,nz) - Ey_old(nx,nz)...
          + normDt^2 * ( Ey(nx,nz+1) + Ey(nx+1,nz) + Ey(nx-1,nz) + Ey(nx,nz-1) );
      end
  end
  
  %// Apply PEC boundary condition at left and right boundary
  for nz=1:Nz
    Ey_new(1,nz) = 0;
  end
  for nz=1:Nz
    Ey_new(Nx,nz) = 0;
  end
  for nx=1:Nx
    Ey_new(nx,1) = 0;
  end
  for nz=1:Nx
    Ey_new(nx,Nz) = 0;
  end
    
  switch excitation
    case 1
      
      %// 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 excitation pulse as an impress electric current density
      Jex(nt) = pulse(nt);

      %// Apply excitation pulse as an impress electric current density
      if nt > 1 
        Ey_new(Nx_Source,Nz_Source) = Ey_new(Nx_Source,Nz_Source) - normDt * ( Jex(nt) - Jex(nt-1) );
      else
        Ey_new(Nx_Source,Nz_Source) = Ey_new(Nx_Source,Nz_Source) - normDt * Jex(nt);
      end
    case 2
      
      nz = 2;
      for nx=round(0.1*Nx):round(0.5*Nx)
        t_ret = Dx*(nx-round(0.1*Nx))*sin(theta_in*M_PI/180.0)/c_0;
        
        %// Apply excitation pulse
        %// Raised cosine pulse with two cycles - (RC2) pulse
        t = Dt_ref*normDt*(nt-1) - t_ret; %// 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 excitation pulse as an impress electric current density
        Jex(nx,nt) = pulse(nt);
        
        %// Apply excitation pulse as an impress electric current density
        if nt > 1 
          Ey_new(nz,nx) = Ey_new(nz,nx) - normDt * ( Jex(nx,nt) - Jex(nx,nt-1) );
        else
          Ey_new(nz,nx) = Ey_new(nz,nx) - normDt * Jex(nx,nt);
        end
      end
  end
  
  if debug == 1    
    fprintf('t= %8.4f pulse= %8.4f  Ey[N_Source]= %8.4f\n', t, pulse(nt), Ey_new(Nx_Source,Nz_Source));
  end

  for nz=2:Nz-1
    for nx=2:Nx-1
      Ey_old(nx,nz) = Ey(nx,nz);
      Ey(nx,nz)     = Ey_new(nx,nz);
    end
  end
  
  for n=1:Nrec
    Eyrec(nt,n)  = Ey(Nx_Source+round((n-1)*Nx/(2*Nrec)),Nz_Source);
  end
   
  Ey_max = max(max(Ey));  Ey_min = min(min(Ey));
  %  if debug == 1
  fprintf('Ey_min = %20.10g Ey_max = %20.10g\n', Ey_min, Ey_max);
  % end
      
  nwritten = fwrite(fp_ex_bin,Ey,'float32');
  %// fprintf('%d byte written\n', nwritten);
  %// fprintf('%d byte written\n', nwritten);
  nwritten = fwrite(fp_pulse_bin,pulse(nt),'float32');
  %// fprintf('%d byte written\n', nwritten);
  
  %// Generate output figures
  
  nt_snapshot_count = nt_snapshot_count + 1;
  
  if nt_snapshot_count >= nt_snapshot
        
    figure(1)
    plot(Jex,'b','LineWidth',1.5)
    t = Dt_ref*normDt*(nt-1+0.5); %// Actual time
    title(['J_{ey} 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('t in s');    
    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_ey');
    grid on 
    
    figure(3)
    %  for n=1:Nrec
    plot(Eyrec(:,2),'b','LineWidth',1.5);
    hold on;
    plot(Eyrec(:,3),'r','LineWidth',1.5);
    hold on;
    plot(Eyrec(:,4),'g','LineWidth',1.5);
    hold on;
    
    %  end
  hold off;
  t = Dt_ref*normDt*(nt-1+0.5); %// Actual time
  title(['J_{ey} component at t = ',num2str(t*1e12),' ps']);
  axis([1 Nt -0.25 0.25]);   
  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 s');    
  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_ey');
  grid on 
  
    
    figure(2)
    pcolor(Ey)
    caxis([-0.07 0.07]);
    shading flat;  
    colorbar;
    t = Dt_ref*normDt*(nt-1+0.5); %// Actual time
    title(['E_x component at t = ',num2str(t*1e12),' ps']);
    axis square
    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 Nz/4 Nz/2 3*Nz/4 Nz],'FontSize',8);
    set(gca, 'YTickLabel',[0 Z/4 Z/2 3*Z/4 Z],'FontSize',8);
    ylabel('z in m');
    grid on  
    
    nt_snapshot_count = 0;
    
    s = sprintf('fd2dtm_vac_rcn_poyn_%d.tif',nt);
       
    fprintf('Generate a TIFF output file of Figure 1\n');
    print('-f2','-dtiff','-r150',s); %// Generate a EPS output file of figure 1
    
  end
    
  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 a TIFF output file of Figure 1\n');
print('-f1','-dtiff','-r150','rc2.tif'); %// Generate a EPS output file of figure 1
