%% Traces a straight ray in the given model between arbitrary given source and receiver
%% ======================================================================================
%%
%% Parameters:
%%	model - structure describing the 2D rectangular model. 
%%		The structure should contain the following fields:
%%		model.xz0  - array of (x,z) coordinates of the upper-left corner
%%		model.dxz  - array of (dx,dz) increments in the grid (cell dimensions)
%%		model.nx and model.nz - numbers of grid cells 
%%	src - 2-element array containing source (x,z)	
%%	rec - 2-element array containing receiver (x,z)
%%
%% Returns:
%%	ray - matrix, each column of which contains parameters of one model
%%		cell crossed by the ray from 'src' to 'rec' 
%%		The column contains:
%%			row 1 - cell number, counting by scanlines 
%%				(1:nx for the top model row, (nx+1):2*nx for the second row, etc.)
%%			row 2 - length of ray segment in that row;
%%			rows 3 and 4: coordinates (x,z) of the entry point into this cell
%%				(closer to source 'src');
%%			rows 5 and 6: coordinates (x,z) of the exit point from this cell
%%				(closer to receiver 'rec')

function ray = trace_ray(model,src,rec)

  ray = [];	% start with empty output

  raypath = [src; rec];	% matrix containing two ends of the ray as rows

	% for each model cell, if it is crossed by the ray,
	% add parameters of intersection to output

  k = 0;

  raydist = [];		% this matrix will contain distances from the source for the ray

  for i=1:model.nz
    for j=1:model.nx

	k = k + 1;		% model cell number

	% determine four corners of the cell:      p1 ----- p2
	%                                          |        |
	%                                          p3 ----- p4

	p4 = model.xz0 + [ j*model.dxz(1), i*model.dxz(2) ];	% lower-right corner	
	p1 = p4 - model.dxz;					% upper-left 
	p2 = p1 + [ model.dxz(1), 0 ];				% upper-right
	p3 = p1 + [ 0, model.dxz(2) ];				% lower-left
		
	% determine whether source and/or receiver are inside this cell
	
	src_inside = src(1) >= p1(1) && src(2) >= p1(2) && src(1) <= p4(1) && src(2) <= p4(2);
	rec_inside = rec(1) >= p1(1) && rec(2) >= p1(2) && rec(1) <= p4(1) && rec(2) <= p4(2);
	
	% if both source and receiver are inside the cell, the segment has been found!

	if src_inside && rec_inside
		segment = [ src; rec, distance(src,rec)];
	else

		% otherwise locate intersections of the ray with four sides of the cell:
		% the lambda_... outputs are only for checking; they are not needed for output

		[ top lambda_top] = intersection(raypath,[p1; p2]);
		[ left lambda_left] = intersection(raypath,[p1; p3]);
		[ bottom lambda_bottom ] = intersection(raypath,[p3; p4]);
		[ right lambda_right ] = intersection(raypath,[p2; p4]);

		segment = [];		% ray segment within cell

		% if intersections exist,
		% add to 'segment' rows consisting of (x,z) coordinates of intersections
		% and distances from 'src'

		if ~isempty(top)
			segment = [segment; top, distance(src,top) ];
		end

		if ~isempty(bottom)
			segment = [segment; bottom, distance(src,bottom) ];
		end

		if ~isempty(left)
			segment = [segment; left, distance(src,left) ];
		end

		if ~isempty(right)
			segment = [segment; right, distance(src,right) ];
		end

		% if the ray does not cross the edge of the cell, try another cell

 		if size(segment,1) < 1
			continue
		elseif size(segment,1) < 2	% only one crossing

			% if source is inside this cell, count it as entry, and the intersection is exit

			if src_inside
				segment = [ 	src, 0.0 
						segment(1,1:2), distance(src,segment(1,1:2))
						segment ]; 

				% otherwise if receiver is inside this cell, count it as exit
			elseif rec_inside
				segment = [ segment; rec, distance(segment(1,1:2),rec) ]; 
			else	% this should not really happen; just continue for safety
				continue
			end
		end
	end

		% test printout for a cell which seems to cause trouble

%	if k == 17
%		top
%		left
%		bottom
%		right
%		segment
%	end

	% if we got here, sort 'segment' by increasing distance from 'src'

	[s ind] = sort(segment(:,3));	% third column in 'segment' is distance

	segment = segment(ind,:);	% sorted

	% output 'segment' as an additional column into 'ray' 

	entry_p = segment(1,:);		% entry point
	exit_p = segment(end,:);	% exit point

	raycol = [	k				% cell number
			exit_p(3) - entry_p(3)		% length of ray segment
			transpose(entry_p(1:2))		% entry point as a column
			transpose(exit_p(1:2))		% exit point as a column
		];

	ray = [ ray, raycol ];

	% also accumulate a list of entry and exit distances from the source

	raydist = [ raydist, [ entry_p(3); exit_p(3) ] ];

    end	% for loop nx
  end	% for loop nz

	% In the resulting output, some ray segments may still be shared by different 
	% model cells (this happens if the ray passes along a boundary between two cells.
	% In such cases, average the output ray-segment lengths among these cells.

  nsegm = size(ray,2);

  for i=1:nsegm
    for j=(i+1):nsegm
	
	l = max( [ raydist(1,i), raydist(1,j)] );	% max of left ends
	r = min( [ raydist(2,i), raydist(2,j)] );	% min of right ends
 
	if l < r					% have an overlap
	  o2 = 0.5*(r - l);

#o2
#i
#rdi = raydist(:,i)
#rayi = ray(:,i)

#j
#rdj = raydist(:,j)
#rayj = ray(:,j)

		% remove half of the overlap from each of the competing cells

		ray(2,i) = ray(2,i) - o2;
		ray(2,j) = ray(2,j) - o2;
	end

    end	% loop j

  end % loop i

end