function TIME = rfltim(D,V,DP,NOL,XDIST,TEST)
% RFLTIM calculates reflection travel times for all layers of an 
% NOL-layered model (D = thicknesses, V = velocities, DP = dips)
% at distance XDIST by iteration on angle of incidence.
% Modified from FORTRAN code in TXCURV.
% TEST is distance tolerance for convergence (in units of distance)
% L. Braile 2/6/03
if (TEST < 0.0000001) TEST=0.001; end
NL=NOL-1;
PI2=pi/2;
STANGL=0.0;
for I = 1:NL
    DIP(I)=DP(I)*0.017453; 
end
for I = 1:NL  % 1002 loop for layers        
ITNUM=0;  % set iteration counter
ITERFAIL = 0;
ANGLE=STANGL*0.017453;
DELANG=0.5000001*0.017453;
XXDIST = -1;
while abs(XDIST-XXDIST) > TEST
NS = 0;
XXDIST=0.0;
YYDIST=0.0;
TIMED=0.0;
ITNUM=ITNUM+1;
if ITNUM > 1000
    ITERFAIL = 1;
    break
end
for JJ = 1:I  % 1004 loop for angles
if JJ == 1; A1(I,JJ)=ANGLE; end
if JJ ~= 1 
if V(JJ-1) <= V(JJ);
    ALPHB=asin(V(JJ-1)/V(JJ)); 
    if (A2(I,JJ-1) > ALPHB); 
        NS = 1; 
        break 
    end  %  exit while loop, equivalent to got to 85
end  %  3003 continue
A1(I,JJ)=asin(V(JJ)*sin(A2(I,JJ-1))/V(JJ-1));
end  %  100 continue
if JJ == 1; YYDIST=D(JJ); end
if JJ ~= 1
YYDIST=D(JJ)+XXDIST*sin(DIP(JJ-1))/cos(DIP(JJ-1))-XXDIST*sin(DIP(JJ))/cos(DIP(JJ)); 
end  %  101 continue
if JJ == 1 A2(I,JJ)=A1(I,JJ)-DIP(JJ); end
if JJ ~= 1
    A2(I,JJ)=A1(I,JJ)-DIP(JJ)+DIP(JJ-1);
end
XDOWN(JJ)=YYDIST*sin(PI2-DIP(JJ))/sin(PI2-A2(I,JJ));
XXDIST=XXDIST+XDOWN(JJ)*sin(A2(I,JJ)+DIP(JJ));
TIMED=TIMED+XDOWN(JJ)/V(JJ);
if JJ ~= I
    continue
end
for JK = 1:I  % 1005 loop
    JKK=I-JK+1;
if JK <= 1 
    B2(I,JKK)=A2(I,JKK);
end
if JK == I
    YYDIST=D(JKK)-XXDIST*sin(DIP(JKK))/cos(DIP(JKK));
end
if JK ~= I 
YYDIST=D(JKK)-XXDIST*sin(DIP(JKK))/cos(DIP(JKK))+XXDIST*sin(DIP(JKK-1))/cos(DIP(JKK-1)); 
end
if JK ~= 1 
     B2(I,JKK)=asin(V(JKK)*sin(B1(I,JKK+1))/V(JKK+1));
end
if JKK ~= 1 
      B1(I,JKK)=B2(I,JKK)-DIP(JKK)+DIP(JKK-1); 
            else 
      B1(I,JKK)=B2(I,JKK)-DIP(JKK);
end
if JKK ~= 1 
      XUP(JKK)=YYDIST*sin(PI2+DIP(JKK-1))/sin(PI2-B1(I,JKK));
            else 
      XUP(JKK)=YYDIST/sin(PI2-B1(I,JKK));
end
XXDIST=XXDIST+XUP(JKK)*sin(B2(I,JKK)-DIP(JKK));
TIMED=TIMED+XUP(JKK)/V(JKK);
  end  % 1005 CONTINUE
end  % 1004 CONTINUE
if NS == 1
    ANGLE=ANGLE+DELANG;
    DELANG = DELANG/2;
continue
end
if XXDIST <= XDIST
    ANGLE=ANGLE+DELANG;
else
    ANGLE=ANGLE-DELANG;
    DELANG = DELANG/2;
end
end
TIME(I)=TIMED;
if ITERFAIL == 1
    TIME(I) = NaN;
end
% ITNUM  %  to check max number of iterations per layer per reflection
end % 1002 CONTINUE