PROGRAM RESIDENCE; 
VAR x,y,t,dst:array[0..10000] of extended;
  r,r2,rs2,ri2,u,v,a,c,s,frt,brt,d2,tf,tb,tin,tout,rt,
  tstep,sstep,d,tc,tp,xp,yp,xc,yc,dtn,dtp,dist,ra:extended;
  n,i,j,k,jout,kin,ic,max,red:integer;
  fi,fo:text; nif,nof:string;

FUNCTION arctan2(yy,xx:extended):extended;
//GEOMETRIC ARCTANGENT
  begin
  if abs(xx)<0.000000000001 then begin
    if abs(yy)<0.000000000001 then arctan2:=0
    else if yy<0 then arctan2:=-pi/2 else arctan2:=pi/2; end
  else if xx>0 then arctan2:=arctan(yy/xx)
    else arctan2:=arctan(yy/xx)+pi; end;

FUNCTION time(i1,i2,id,sign:integer):extended;
//LINEAR INTERPOLATION OF ENTER OR EXIT TIME
  begin
  a:=arctan2(y[i2]-y[i1],x[i2]-x[i1]); c:=cos(a); s:=sin(a);
  u:=(x[i]-x[i1])*c+(y[i]-y[i1])*s;
  v:=(y[i]-y[i1])*c-(x[i]-x[i1])*s;
  ra:=(u+sign*sqrt(r2-sqr(v)))/dst[id];
  time:=t[i1]*(1-ra)+t[i2]*ra; end;

BEGIN //MAIN PROGRAM
write('name of input file ?   '); readln(nif);
write('name of output file ?  '); readln(nof);
assign(fi,nif); reset(fi); assign(fo,nof); rewrite(fo);
readln(fi,xp,yp,tp); n:=0; x[n]:=xp; y[n]:=yp; t[n]:=tp;

//OPTIONAL REDISCRETIZATION
write('rediscretization: no=0, spatial=1, temporal=2  ?  ');
readln (red);
case red of
  0: begin //NO REDISCRETIZATION
    while not seekeof(fi) do begin
      inc(n); readln(fi,x[n],y[n],t[n]);
      dst[n]:=sqrt(sqr(x[n]-x[n-1])+sqr(y[n]-y[n-1]));
      end; end;
  1: begin //SPATIAL REDISCRETIZATION
    write('step length ?  '); readln(sstep);
    r2:=sqr(sstep); ri2:=0.99998*r2; rs2:=1.00002*r2;
    while not seekeof(fi) do begin
      readln(fi,xc,yc,tc); d2:=sqr(xc-x[n])+sqr(yc-y[n]);
      while d2>ri2 do begin
        inc(n); dst[n]:=sstep;
        if d2<rs2 then begin
          x[n]:=xc; y[n]:=yc; d2:=0; t[n]:=tc; end
        else begin //SEE APPENDIX IN BOVET & BENHAMOU 1988
          a:=arctan2(yc-yp,xc-xp); c:=cos(a); s:=sin(a);
          u:=(x[n-1]-xp)*c+(y[n-1]-yp)*s;
          v:=(y[n-1]-yp)*c-(x[n-1]-xp)*s;
          d:=u+sqrt(r2-sqr(v)); x[n]:=xp+d*c; y[n]:=yp+d*s;
          dist:=sqrt(sqr(xc-xp)+sqr(yc-yp)); ra:=d/dist;
          t[n]:=tp*(1-ra)+tc*ra;
          d2:=sqr(xc-x[n])+sqr(yc-y[n]); end; end;
      xp:=xc; yp:=yc; tp:=tc; end; end;
  2: begin //TEMPORAL REDISCRETIZATION
    write('step duration ?  '); readln(tstep);
    while not seekeof(fi) do begin
      readln(fi,xc,yc,tc); dtn:=tc-t[n]; dtp:=tc-tp;
      while dtn>tstep do begin
        inc(n); t[n]:=t[n-1]+tstep;
        dtn:=tc-t[n]; ra:=dtn/dtp;
        x[n]:=xp*ra+xc*(1-ra); y[n]:=yp*ra+yc*(1-ra);
        dst[n]:=sqrt(sqr(x[n]-x[n-1])+sqr(y[n]-y[n-1])); end;
      xp:=xc; yp:=yc; tp:=tc; end; end; end;

//COMPUTATION OF RESIDENCE TIME (RT)
write('radius of the virtual circle ?  '); readln(r);
write('outside step number allowed ? '); readln(max);
r2:=sqr(r); ri2:=0.99998*r2; rs2:=1.00002*r2; ic:=0;
for i:=1 to n-1 do begin
  j:=i; k:=i;
  //FORWARD RESIDENCE TIME (FRT)
  repeat inc(j); d2:=sqr(x[i]-x[j])+sqr(y[i]-y[j]);
    until (d2>ri2) or (j=n);
  if d2<ri2 then tf:=-1 else if d2<rs2 then tf:=t[j]
  else tf:=time(j-1,j,j,1); //FORWARD FPT
  frt:=tf-t[i]; jout:=j;
  while (j<n-1) and ((j-jout)<max) do begin
    repeat inc(j); d2:=sqr(x[i]-x[j])+sqr(y[i]-y[j]);
      until (d2<rs2) or (j=n-1) or ((j-jout)=max);
    if d2>rs2 then tin:=-1 else if d2>ri2 then tin:=t[j]
    else tin:=time(j-1,j,j,-1);
    if tin>0 then begin
      repeat inc(j); d2:=sqr(x[i]-x[j])+sqr(y[i]-y[j]);
        until (d2>ri2) or (j=n);
      if d2<ri2 then tout:=t[n] else if d2<rs2 then tout:=t[j]
      else tout:=time(j-1,j,j,1);
      frt:=frt+tout-tin; jout:=j; end; end;
  //BACKWARD RESIDENCE TIME (BRT)
  repeat dec(k); d2:=sqr(x[i]-x[k])+sqr(y[i]-y[k]);
    until (d2>ri2) or (k=0);
  if d2<ri2 then tb:=-1 else if d2<rs2 then tb:=t[k]
  else tb:=time(k+1,k,k+1,1); //BACKWARD FPT
  brt:=t[i]-tb; kin:=k;
  while (k>1) and ((kin-k)<max) do begin
    repeat dec(k); d2:=sqr(x[i]-x[k])+sqr(y[i]-y[k]);
      until (d2<rs2) or (k=1) or ((kin-k)=max);
    if d2>rs2 then tout:=-1 else if d2>ri2 then tout:=t[k]
    else tout:=time(k+1,k,k+1,-1);
    if tout>0 then begin
      repeat dec(k); d2:=sqr(x[i]-x[k])+sqr(y[i]-y[k]);
        until (d2>ri2) or (k=0);
      if d2<ri2 then tin:=t[0] else if d2<rs2 then tin:=t[k]
      else tin:=time(k+1,k,k+1,1);
      brt:=brt+tout-tin; kin:=k; end; end;
  //RESIDENCE TIME = FORWARD RT + BACKWARD RT
  if (tb>0) and (tf>0) then begin
    inc(ic); rt:=frt+brt;
    writeln(fo,ic:9,x[i]:9:3,y[i]:9:3,t[i]:9:3,rt:9:3);
    end; end;
close(fi); close(fo); END.