function gpr_reflect_v1_ex1

% Simulates normal incidence GPR reflection.

% Presently for 3 layers - easy to generalize
% No geoemetrical spreading or atanuation included
% Could add gain

% Start ray at depth to simulate receiver functions
% z_loc = 200;  idir(1) = -1   ilayer = 3

clear all; close all;

% Define velocity model in m per nanosecond
v =  [0.06, 0.15];  % Dry sand over wet

% These are node numbers
interface =  [600];           ninterface = length(interface);
plot_limit = 3;

% Maximum number of depth nodes
nzstep = 1000;   

% ilayer denotes which layer each ray is presently in

% No. of time steps;   Size of time step (ns); Total time(ns)
ntstep=1500;  	       dt =0.02;               total_t = ntstep*dt;

% Convert these parameters into true depth
z_start = 0;
for izstep = 1:interface(1)
   z(izstep) = z_start;
   v_plot(izstep)   = v(1);
   z_start = z_start + v(1)*dt;
end

if ninterface > 1
  for izstep = interface(1)+1:interface(2)
    z(izstep) = z_start;
    z_start = z_start + v(2)*dt;
    v_plot(izstep)   = v(2);
  end
end


for izstep = interface(ninterface)+1:nzstep
  z(izstep) = z_start;
  z_start = z_start + v(ninterface+1)*dt;
  v_plot(izstep)   = v(ninterface+1);
end

zmax = z(nzstep)


% Define lines to plot on model to show layers
for i=1:ninterface
  zplot(i,1) = z(interface(i));  zplot(i,2) = z(interface(i));
  xplot(i,1) = -plot_limit ;  xplot(i,2) = plot_limit ;
end

% Define source wavelet
nwavelet=20
for i=1:nwavelet
  f(i)=sin(pi*(i-1)/(nwavelet-1));
end




% Set up wave. Each wave defined by    (a) location (node number)  z_loc
%                                      (b) direction:              down idir =1
%                                                                  up   idir =-1
%                                      (c) ilayer                  which layer it is in
%                                      (d) amp                     amplitude
%                                      (e) i_terminate             Set to 1 when wave reaches
%                                                                  surface
z_loc(1)  =2;   amp(1)    =1; ilayer(1) =1; idir(1)   =1;   i_terminate(1) = 0;
nray=length(z_loc);

% ===================================================================
% Start time step loop
% ===================================================================
for it=1:ntstep;
  z_loc
  %Check if upgoing wave at surface. Send down, amplitude = 0
  for iray =1:nray
    if z_loc(iray) == 1
      if idir(iray) == -1
        amp(iray) =0;
        i_terminate(iray) = 1;
      end
    end 
 end
 
  %Stop ray at bottom of model  
  for iray =1:nray
    if z_loc(iray) == nzstep
      if idir(iray) == 1
        amp(iray) =0;
        i_terminate(iray) = 1;
      end
    end 
  end 
   
  % loop over all rays and move position according to direction 
  for iray=1:nray
    if i_terminate(iray) == 0
      z_loc(iray)=z_loc(iray)+idir(iray);
    end
  end

  % Loop over rays to see if interface has been reached.
  n_new =0;
  
  for iray=1:nray
    for i=1:ninterface
      if z_loc(iray)==interface(i)
        % Original wave transmits, generate new reflected wave. 
         
        % Down going wave
        if idir(iray) == 1
          imp1 = v(ilayer(iray));         imp2 = v(1+ilayer(iray)); 
          ref = (imp2-imp1)/(imp2+imp1);  tran= (2*imp1)/(imp2+imp1);
           
          % New reflected wave 
          z_loc(nray+1) = interface(i);
          idir(nray+1) = -idir(iray);
          ilayer(nray+1) = ilayer(iray);
          amp(nray+1) = amp(iray)*ref;
          i_terminate(nray+1) = 0
           
          % Transmit origianl wave
          %idir(iray)= -idir(iray);
          amp(iray) = amp(iray)*tran;
          ilayer(iray)=ilayer(iray)+1;       
           
          n_new = n_new+1; 
        end
        
        % Upgoing wave
        if idir(iray) == -1   
          imp2 = v(ilayer(iray)-1);        imp1 = v(ilayer(iray));   
          ref = (imp2-imp1)/(imp2+imp1);   tran= (2*imp1)/(imp2+imp1);
          % New reflected wave 
          z_loc(nray+1) = interface(i);
          idir(nray+1) = -idir(iray);
          ilayer(nray+1) = ilayer(iray);
          amp(nray+1) = amp(iray)*ref;
          i_terminate(nray+1) = 0;
           
          % Transmit original wave
          amp(iray) = amp(iray)*tran
          ilayer(iray)=ilayer(iray)-1;            
           
          n_new = n_new+1;                 
        end   
      end
    end
  end
   
  nray = nray+n_new;
      
  % Generate displacement, 
  for izstep=1:nzstep
    disp(izstep) = 0;
  end
   
  for iray =1:nray
    disp(z_loc(iray)) = disp(z_loc(iray))+ amp(iray);
  end
  
  % Extract displacement at receiver at this time step for seismic trace
  time(it) = it*dt;
  x(it) = disp(1);
   
%   Convolve wave in depth domain with source wavelet
  dcon = conv(disp,f);
  nshift = round(nwavelet/2);
  for izstep=1:nzstep
    disp(izstep)=dcon(izstep+nshift);
  end
  
% Convolve seismic trace with source wavelet  
  xcon = conv(x,f);
  for i=1:length(time)
    trace_plot(i)=xcon(i+nshift);
  end
       
%=================================================================   
% Display results for this time step   
%=================================================================
   
% Plot Earth model  with choice of waves and wavefronts (every 10 steps)

  if round (it/20) == it/20;
     
    % Plot waveform in Earth 
    subplot(1,3,2)   
      plot(disp,-z) 
      hold on
      for i=1:ninterface
        plot(xplot(i,:),-zplot(i,:),':')
      end
      hold off
      axis([-1.2,1.2,-zmax,0])
      xlabel('amplitude')

    % Plot seismic trace
    subplot(1,5,5)
      plot(trace_plot,-time)
      axis([-1.2,1.2,-total_t,0])
      ylabel('time (ns)')
      xlabel('amplitude')
      
    % Plot velocity model  
    subplot(1,3,1)
      plot(v_plot,-z)
      axis([0,1.1*max(v_plot),-zmax,0])
      xlabel('v (m/ns)') 
      ylabel('depth (m)')
      pause(0.2)
   end
% End of time step loop
end   

print -dpsc  gpr_reflect_v1_ex1.ps
print -djpeg100  gpr_reflect_v1_ex1.jpg