function[v]=ray(xl,yl,xr,yr,dx,dy,kx,ky)
% [v]=ray(xl,yl,xr,yr,dx,dy,kx,ky)
% gives a line in the tomography corresponds to transmitter in xl,yl, 
% and a receiver in xr,yr. The media is divided into kx times ky cells
% of size dx,dy
% Input: xl - x position of transmitter
%        yl - y position of transmitter
%        xr - x position of receiver
%        yr - y position of receiver
%        dx - size of a cell (x direction)
%        dy - size of a cell (y direction)
%        kx - How many cells in the x direction
%        ky - How many cells in the y direction
%
% Output: v - A vector size (kx*ky) of ray lengths.
% Remark The transmitters should always be on the left and the
%        receivers on the right as follows:
%
%        |                |
%       s*                |
%        |                |
%        |                >r
%        |                |
%


% Check that transmitters are on the left and receivers on the 
% right an if not change between them

% For plotting the raypaths, remove the comment signs on the
%    the code at the end.


%====== Preparations =====================================================
 if xr<xl, a=xr;xr=xl;xl=a;a=yr;yr=yl;yl=a; end;

% ==== First possibility transmitter and receiver have the same y ===========

  if yl==yr
       xpoint=[0:dx:kx*dx];
       ypoint=ones(1,kx+1)*yl;
  else,end;

%=====================Step 1 ============================================= 
% Calculating the slope of each ray and its intersection

  a=(yr-yl)/(xr-xl);
  b=yl-a*xl;

% added Feb. 2009 by reso, check if a=0, if so, then add a small
% perturbation to the location of the receiver location.  However, this 
% problem is caused when a=0, is not trapped by any of the if/elseif logic in
% Step 2 of the code (see below).  Simply changing one of  a>0 or a<0 to a>=0 or a<=0
% should trap this condition. However, there may be implications to doing
% this.  This "perturbation" method will always place the ray into the
% "lower cell" when the ray runs along a boundary, and has logic to ensure
% that the receiver isn't positioned.  Note that this "hack" is only
% activated when yl == yr so it shouldn't greatly affect the solution to
% the problem.
if(a == 0)
    yr = yr + 1e-6;
    ind=find(yr > dy*ky);
    yr(ind) = yr - 2*1e-6;  % check if it is below the grid

    if xr<xl, a=xr;xr=xl;xl=a;a=yr;yr=yl;yl=a; end;
    if yl==yr
        xpoint=[0:dx:kx*dx];
        ypoint=ones(1,kx+1)*yl;
    else,end;
    a=(yr-yl)/(xr-xl);
    b=yl-a*xl;
 end;  
%======================Step 2 ============================================ 
% Finding and Sorting all the intersections
     
  x=[ceil(xl/dx)*dx:dx:floor(xr/dx)*dx];
  if a>0,
        y=[ceil(yl/dy)*dy:dy:floor(yr/dy)*dy];
  elseif a<0,
        y=[floor(yl/dy)*dy:-dy:ceil(yr/dy)*dy];
  else, end;
  
  xy=(y-b)/a;
  xpoint=[x,xy];
  xpoint=sort(xpoint);

  
  ypoint=a*xpoint+b;
 
%===================== Step 3 =============================================
% Assigning cell number for each intersection
 
  xmidpnt=xpoint+0.5*[diff(xpoint),0];
  ymidpnt=ypoint+0.5*[diff(ypoint),0];

  cellnum=floor(ymidpnt/dy)*kx+floor(xmidpnt/dx)+1;
%==================== Step 4 ==============================================  
% Putting all this in a vector
 
  v=zeros(1,kx*ky);
  
  for i=1:length(cellnum)-1, 
       
       if (cellnum(i)<=kx*ky & cellnum(i)>0)
             v(1,cellnum(i))=((ypoint(i+1)-ypoint(i))^2+...
               (xpoint(i+1)-xpoint(i))^2)^0.5;
       else
       end;
   end;


%================= Done Calculation ======================================   
% If there are any trouble take the remark of the next part  
% to see each ray.

             %figure(1)
              %[X,Y]=meshgrid(0:dx:dx*kx,0:dy:dy*ky); 
              %plot(X,Y); 
              %hold on 
              %plot(Y,X); 
              %axis([0,kx*dx,0,ky*dy]);
              %plot(xpoint,ypoint) 
             %plot(xmidpnt,ymidpnt,'o');
              %axis('ij') 
 
	      %print -dps raypath.ps;