{$N+,E+} {* Link with full 8087 emulator, uses 8087 if present *} UNIT SPLINES; {***********************************************} {* *} {* SPLINE ROUTINES *} {* *} {* J S Perry *} {* *} {* Filename - SPLINES.PAS Pascal v6 *} {* *} {***********************************************} INTERFACE uses graphgp; const coord_array_size = 60; {size of array to store series of points on curves} type coord_array=array[0..coord_array_size] of real; var nsteps:integer; {number of intemediate steps on spline curves} procedure CURVE(x,y:coord_array; i,j,coordsystem:integer); Procedure INTERCEPT_CURVE(x,y:coord_array; i,j:integer;xx:real; var nintercept:integer; var yy:real); procedure SCURVE(x,y:coord_array; i,j,coordsystem:integer); procedure OCURVE(x,y:coord_array; i,j,coordsystem:integer); IMPLEMENTATION type array5 = array[1..5] of real; procedure ERRORMESSAGE; begin graphgpoff;writeln('curve with only one point');gpause;halt; end; procedure SPLINE(a,b:array5;coordsystem:integer); var j:integer; t,xx,yy:real; begin movto(coordsystem,a[1],b[1]); for j:=1 to nsteps do begin t:=j/nsteps; xx:=a[1]+t*(a[2]+t*(a[3]+t*a[4])); yy:=b[1]+t*(b[2]+t*(b[3]+t*b[4])); linto(coordsystem,xx,yy); end; end; procedure CURVE(x,y:coord_array; i,j,coordsystem:integer); var a,b:array5; xa,ya,sa,xb,yb,sb,q,xgradb,ygradb,xgrada,ygrada:real; k:integer; begin movto(coordsystem,x[i],y[i]); if (j-i)=0 then errormessage; if (j-i)=1 then linto(coordsystem,x[j],y[j]); if (j-i)>1 then begin xa:=x[i+1]-x[i]; ya:=y[i+1]-y[i]; sa:=sqrt(xa*xa+ya*ya); for k:=i to j-1 do begin if k<>j-1 then begin xb:=x[k+2]-x[k+1]; yb:=y[k+2]-y[k+1]; sb:=sqrt(xb*xb+yb*yb); q:=sb/sa; xgradb:=(xb/q+q*xa)/(1.0+q); ygradb:=(yb/q+q*ya)/(1.0+q); end; if k=j-1 then begin a[1]:=x[k]; a[2]:=xgrada; a[3]:=xa-xgrada; a[4]:=0; b[1]:=y[k]; b[2]:=ygrada; b[3]:=ya-ygrada; b[4]:=0; end; if k=i then begin a[1]:=x[k]; a[2]:=2*xa-xgradb; a[3]:=xgradb-xa; a[4]:=0; b[1]:=y[k]; b[2]:=2*ya-ygradb; b[3]:=ygradb-ya; b[4]:=0; end; if (k<>i) and (k<>j-1) then begin a[1]:=x[k]; a[2]:=xgrada; a[3]:=3*xa-2*xgrada-xgradb; a[4]:=xgrada+xgradb-2*xa; b[1]:=y[k]; b[2]:=ygrada; b[3]:=3*ya-2*ygrada-ygradb; b[4]:=ygrada+ygradb-2*ya; end; spline(a,b,coordsystem); xa:=xb; ya:=yb; sa:=sb; xgrada:=q*xgradb; ygrada:=q*ygradb; end; end; end; function INBETWEEN(x,x1,x2:real):boolean; var a:real; begin if x1>x2 then begin a:=x1; x1:=x2; x2:=a; end; inbetween:=false; if ((x>=x1) and (x2>x)) then inbetween:=true; end; procedure INTERPOLATE(a,b:array5;xx:real; var nintercept:integer; var yy:real); var j:integer; t,x1,y1,x2,y2:real; begin x1:=a[1]; y1:=b[1]; for j:=1 to nsteps do begin t:=j/nsteps; x2:=a[1]+t*(a[2]+t*(a[3]+t*a[4])); y2:=b[1]+t*(b[2]+t*(b[3]+t*b[4])); if ((xx>=x1) and (xx0 then begin if (j-i)=1 then begin if inbetween(xx,x[i],x[j]) then begin nintercept:=1; yy:=y[i]+(xx-x[i])*(y[j]-y[i])/(x[j]-x[i]); end; end; if (j-i)>1 then begin xa:=x[i+1]-x[i]; ya:=y[i+1]-y[i]; sa:=sqrt(xa*xa+ya*ya); for k:=i to j-1 do begin if k<>j-1 then begin xb:=x[k+2]-x[k+1]; yb:=y[k+2]-y[k+1]; sb:=sqrt(xb*xb+yb*yb); q:=sb/sa; xgradb:=(xb/q+q*xa)/(1.0+q); ygradb:=(yb/q+q*ya)/(1.0+q); end; if k=j-1 then begin a[1]:=x[k]; a[2]:=xgrada; a[3]:=xa-xgrada; a[4]:=0; b[1]:=y[k]; b[2]:=ygrada; b[3]:=ya-ygrada; b[4]:=0; end; if k=i then begin a[1]:=x[k]; a[2]:=2*xa-xgradb; a[3]:=xgradb-xa; a[4]:=0; b[1]:=y[k]; b[2]:=2*ya-ygradb; b[3]:=ygradb-ya; b[4]:=0; end; if (k<>i) and (k<>j-1) then begin a[1]:=x[k]; a[2]:=xgrada; a[3]:=3*xa-2*xgrada-xgradb; a[4]:=xgrada+xgradb-2*xa; b[1]:=y[k]; b[2]:=ygrada; b[3]:=3*ya-2*ygrada-ygradb; b[4]:=ygrada+ygradb-2*ya; end; interpolate(a,b,xx,nintercept,yy); xa:=xb; ya:=yb; sa:=sb; xgrada:=q*xgradb; ygrada:=q*ygradb; end; end; end; end; procedure SCURVE(x,y:coord_array; i,j,coordsystem:integer); var a,b:array5; xa,ya,sa,xb,yb,sb,q,xgradb,ygradb,xgrada,ygrada:real; k:integer; begin movto(coordsystem,x[i],y[i]); if (j-i)=0 then errormessage; if (j-i)=1 then linto(coordsystem,x[j],y[j]); if (j-i)>1 then begin xa:=x[i+1]-x[i]; ya:=y[i+1]-y[i]; for k:=i to j-1 do begin if ki then begin a[1]:=x[k]; a[2]:=xa; a[3]:=0; a[4]:=0; b[1]:=y[k]; b[2]:=ygrada; b[3]:=3*ya-2*ygrada-ygradb; b[4]:=ygrada+ygradb-2*ya; end; if k=i then begin a[1]:=x[k]; a[2]:=xa; a[3]:=0; a[4]:=0; b[1]:=y[k]; b[2]:=2*ya-ygradb; b[3]:=ygradb-ya; b[4]:=0; end; end; if k=j-1 then begin a[1]:=x[k]; a[2]:=xa; a[3]:=0; a[4]:=0; b[1]:=y[k]; b[2]:=ygrada; b[3]:=ya-ygrada; b[4]:=0; end; spline(a,b,coordsystem); xa:=xb; ya:=yb; ygrada:=q*ygradb; end; end; end; procedure OCURVE(x,y:coord_array; i,j,coordsystem:integer); var a,b:array5; xa,ya,sa,xb,yb,sb,q,xgradb,ygradb,xgrada,ygrada:real; k,kk:integer; begin movto(coordsystem,x[i],y[i]); if (j-i)=0 then errormessage; if (j-i)=1 then linto(coordsystem,x[j],y[j]); if (j-i)>1 then begin xa:=x[i+1]-x[i]; ya:=y[i+1]-y[i]; sa:=sqrt(xa*xa+ya*ya); for k:=i to j+1 do begin kk:=k; if ki then begin a[1]:=x[kk]; a[2]:=xgrada; a[3]:=3*xa-2*xgrada-xgradb; a[4]:=xgrada+xgradb-2*xa; b[1]:=y[kk]; b[2]:=ygrada; b[3]:=3*ya-2*ygrada-ygradb; b[4]:=ygrada+ygradb-2*ya; spline(a,b,coordsystem); end; xa:=xb; ya:=yb; sa:=sb; xgrada:=q*xgradb; ygrada:=q*ygradb; end; {for k=i to j+1} end; {if j-i>1} end; {ocurve} begin nsteps:=50; end.