function [Ttable]=Mray(SW,SP,DX)
% 2D ray-tracing for migration time table
%
% usage: Ttable=Mray(SW,SP,DX)
%
% input:
%       SW: slowness model 
%       SP: shot position
%       DX: spacing
%
% code by: Ruiqing He, University of Utah, 2002

if nargin==0
    % modeling example 
%    SW=1/1000*ones(200,200); 
%    for i=1:200; SW(i,:) = 1/(1000+10*i); end;
    % call like this
    T=Mray(SW,SP,DX);     
    imagesc(T); return;
end

[Z,X]=size(SW);
def=10000;
delt=max(SW(:));
sZ=7; sX=7;
dA=subs2(sZ,sX);
ZZ=Z+2*sZ-1;XX=X+2*sX-1;
T=ones(ZZ,XX)*def;
mark=ones(ZZ,XX)*def;

Z2=Z+2*sZ-2;X2=X+2*sX-2;
S=ones(Z2,X2);

Z1=sZ:Z+sZ;X1=sX:X+sX;
mark(Z1,X1)=0;

Z2=sZ:Z+sZ-1;
X2=sX:X+sX-1;
S(Z2,X2)=SW;
S(Z+sZ,X2)=2*S(Z+sZ-1,X2)-S(Z+sZ-2,X2);
S(Z2,X+sX)=2*S(Z2,X+sX-1)-S(Z2,X+sX-2);
S(Z+sZ,X+sX)=2*S(Z+sZ-1,X+sX-1)-S(Z+sZ-2,X+sX-2);

dz=-sZ+1;dx=-sX+1;
[is,dum]=size(SP);

z=SP(1);x=SP(2);
z=z+sZ-1;x=x+sX-1;
T(z,x)=0;
mark(z,x)=def;

a=2*sZ-1;b=2*sX-1;
aa=-sZ+1:sZ-1;bb=-sX+1:sX-1;as=-sZ+1:sZ-2;bs=-sX+1:sX-2;
AS=S(as+z,bs+x);
aaa=aa+z;bbb=bb+x;
TT=T(aaa,bbb);
T(aaa,bbb)=min(reshape(dA*AS(:)+T(z,x),a,b),TT);
maxt=max(max(T(z-1:z+1,x-1:x+1)));
    
while 1
   [hz hx]=find(T+mark<=maxt+delt);[hsz,dum]=size(hz);
   for ii=1:hsz
       z=hz(ii);x=hx(ii);
       maxt=max(maxt,T(z,x));
       mark(z,x)=def;    
       AS=S(as+z,bs+x);
       aaa=aa+z;bbb=bb+x;
       TT=T(aaa,bbb);
       T(aaa,bbb)=min(reshape(dA*AS(:)+T(z,x),a,b),TT);
    end 
    if all(mark(Z2,X2)) break;end
end
Ttable=T(Z2,X2)*DX; 

%%%%%%%%%%%%%%%%%%%%%%
function [pz,px,j]=lineseg2(z0,x0,z1,x1)
% code by: Ruiqing He, University of Utah, 2002
dz=(z1-z0);dx=(x1-x0);
sgnz=sign(dz);sgnx=sign(x1-x0);
pz(1)=z0;px(1)=x0;j=2;

if sgnz~=0;for z=z0+sgnz:sgnz:z1-sgnz
    pz(j)=z;px(j)=x0+(z-z0)*dx/dz;j=j+1;
end; end

if sgnx ~=0;for x=x0+sgnx:sgnx:x1-sgnx
    px(j)=x;pz(j)=z0+(x-x0)*dz/dx;j=j+1;
end; end

pz(j)=z1;px(j)=x1;
[px,id]=sort(px);[pz,id]=sort(pz);
if(sgnx==-sgnz);px=fliplr(px);end 

%%%%%%%%%%%%%%%%%%%%%
function [L]=buildL2(L,Z,X,ind,z0,x0,z1,x1)
% code by: Ruiqing He, University of Utah, 2002
[pz,px,j]=lineseg2(z0,x0,z1,x1);
for i=1:j-1
    l=norm([pz(i+1)-pz(i),px(i+1)-px(i)]);
    a=floor((pz(i+1)+pz(i))/2);if a==Z+1;a=Z;elseif a==0;a=1;end;
    b=floor((px(i+1)+px(i))/2);if b==X+1;b=X;elseif b==0;b=1;end;
    L(ind,sub2ind([Z,X],a,b))=l;
end

%%%%%%%%%%%%%%%%%%%%%
function dA=subs2(sZ,sX)
% code by: Ruiqing He, University of Utah, 2002
z=2*sZ-1;x=2*sX-1;
z1=z-1;x1=x-1;
dA=sparse(z*x,z1*x1);
for j=1:z;for i=1:x
    dA=buildL2(dA,z1,x1,sub2ind([z,x],j,i),sZ,sX,j,i);
end;end    
