% Returns reflection cross-section corresponding to all reflections 
% and diffractions from one boundary, assuming reflection amplitude equal 1.
% In contrast to reflection(), the first parameter is the coordinates of
% midpoints.
%
% Input parameters:
%   mp - midpoint positions (can be an array)
%   R - array of receiver positions
%   layer - matrix of boundary points returned by picklayer()
%   velocity - wave velocity

%   Common model parameters are provided by 'global' 

% Returns:
%   section - the resulting seismic section. n-th columns of the
%         matrix contains reflection amplitudes at n-th receiver
%         at times sampled in array 'time'

function [section] = reflection_CMP(mp,R,layer,velocity)

global  time dt
global  line_length depth spacing

  nt = length(time);      % number of time samples
  
  % first create zero traces for each receiver and time
  
  section = zeros(nt,length(R));
  
  if size(layer,1) < 2      % single reflection point
      xp = layer(1,1);
      yp = layer(1,2);
  else        % distribute many (200 in this example) reflection points 
              % along the whole boundary
  
    xp = linspace(layer(1,1),layer(end,1),200)';   % x coordinates (column)
    yp = interp1(layer(:,1),layer(:,2),xp);        % y cordinates (column)
  end
  
%[ xp yp]

  % add reflections from each reflection point to the section
  
  for i=1:length(mp)
    
    for j=1:length(R)
      
      % source coordinate for this midpoint and receiver:
      
      S = 2*mp(i) - R(j);
      
      % distances to the source for each point on the reflector:
          
      dS = sqrt( (xp - S).^2 + yp.^2 );  
      
      % distances to the receiver for each point on the reflector
          
      dR = sqrt( (xp - R(j)).^2 + yp.^2 );  
      
      % two-way geometrical spreading 1/dR/dS:
          
      A = (dS.*dR).^(-1); 

      % two-way travel time:
      
      t = (dS +dR)/velocity;
      
      % place amplitude A at time t in receiver record:
      
      k = round(t/dt);                    % indices in time array
      ind = find( k >= 1 & k <= nt);     % indices of points
                              % producing reflections within time range
      k = k(ind);
      
  %    section(k,j) =  section(k,j) + A(ind);
  
      for k=1:length(t)
        section(:,j) =  section(:,j) + wavelet(A(k),t(k),time)';
      end
    end   % for j
  end     % for i
  
  % finally, scale the whole section so that it has peak amplitude equals 
  % receiver spacing (for plotting)
  
  section = section *(spacing/ max(max(section)));
  
end