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 d2tstep 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 d2rs2 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 d2ri2) or (k=0); if d21) and ((kin-k)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 d20) 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.