% Outputs depth section containing scattering ellipse (actually, circle)
% consisting of points that could produce reflection at time 't' \
% at a single source-receiver point x_SR 
%
% Input parameters:
%   x - grid of coordinate (x) values at which the output is procuced
%   z - grid of depth values at which the output is procuced
%   x_SR - collocated source and receiver coordinate in reflection zero-offset section 
%   t - reflection times at x_SR 
%   u- is the array of amplitude at coordinates x_SR and times t
%   V - wave velocity in the medium
% Returns:
%   section - depth section in the format similar to other seismic sections in the labs
%       (columns correspond to different values of x, rows - to values of z)

function section = scattering_ellipse(x,z,x_SR,t,u,V)

  % Note that this code is very similar to function reflection() in lab #5
  % This is just an inverse of the procedure of finding reflection times from
  % arbitrary reflector positions.
  % ========================================================================
  
  % create a helper array z2 or pointz z shifted upward by half-spacing
 
  z2 = z - 0.5*(z(2)-z(1));
  
  nz = length(z);      % number of depth samples
  
  % first create zero records for each receiver and depth
  
  section = zeros(nz,length(x));
  
  % add reflections at each x_SR, t, and u in the section
  
  for i=1:length(x_SR)
          
    xs = x_SR(i);     % source coordinate
    R = t(i)*V/2.0;      % distance to the source/receiver for this time t

      % find the range of points within x that are closer than R to the source:
           
    ix = find(x < xs + R & x > xs - R);
    
      % in each vertical column, add response to the output section
      
    for j=1:length(ix)
      
      k = ix(j);      % index of source/receiver to use
      
      z1 = sqrt(R^2 - (x(k) - xs)^2);  % reflection depth under this x(k)
      
      % find the first depth level below z1 and assign value 1.0 to it
      
      iz = find(z2 > z1);
      
      if ~isempty(iz)
        section(iz(1),k) = u(i);
      end

    end   % for j
  end     % for i
  
end   % function