% Plots traces and positions of sources
%
% Input parameters:
%   section - trace section on which the semblance is measured
%           This section can be returned by 'reflection()' or 'reflection_CMP()'
%   offsets - source receiver offsets for each trace (column) in 'section'
%   V - array of stacking velocities to measure semblance
% Returns:
%   semb - resulting semblance, also in the form of a trace section

function semb = semblance(section,offsets,V)

global  time dt

  t0 = time';     % column of intercept times
  
    % form a 40-ms smoothing filter to remove ripples in semblance below
    
  dfilt = dt/20;
  filt = [ 0:dfilt:1, (1-dfilt):-dfilt:0 ]';
       
  semb = zeros(size(section,1),length(V));
  
     % for each trial stacking velocity, evaluate eq. (1)

  for i=1:length(V)
    
    u_stack = zeros(size(section,1),1);   %# initialize zero stack
    
    for j=1:length(offsets)
      t = trefl(t0,V(i),offsets(j));    % reflection times at this offset
        
      ind = find(t >= time(1) & t <= time(end));  % indices of samples lying within the given trace
      
      tt = [ ind t0(ind,1) t(ind,1) ];
      
      u = interp1(t0,section(:,j),t(ind,1));  % wave displacements along stacking hyperbola
      
      u_stack(ind,1) = u_stack(ind,1) + u;   % simply stack along reflection hyperbola
    end
          
    semb(:,i) = conv(u_stack.*u_stack,filt,'same');   % square and smooth stacked waveform    

  end  
  
  % scale the semblance for plotting using V as "distances"
  
  if length(V) > 1
    dV = range(V)/(length(V) - 1);      % average velocity increment 
  else
    dV = 1.0;
  end
  
  semb = semb *(2*dV/max(max(semb)));   % factor 2 to visually enhance the peak
  
end