{$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 (xx<x2)) then
    begin
      inc(nintercept,1);
      yy:=y1+(xx-x1)*(y2-y1)/(x2-x1);
    end;
    x1:=x2; y1:=y2;
  end;
end;


Procedure  INTERCEPT_CURVE(x,y:coord_array; i,j:integer;xx:real;
var nintercept:integer; var yy:real);
var a,b:array5;
xa,ya,sa,xb,yb,sb,q,xgradb,ygradb,xgrada,ygrada:real;
k:integer;
begin
  nintercept:=0;
  if (j-i)<>0 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 k<j-1 then
      begin
        xb:=x[k+2]-x[k+1];
        yb:=y[k+2]-y[k+1];
        q:=xb/xa;
        ygradb:=(yb/q+q*ya)/(1.0+q);
        if k>i 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 k<j-1 then
      begin
        xb:=x[k+2]-x[k+1];
        yb:=y[k+2]-y[k+1];
      end;
      if k=j-1 then   
      begin
        xb:=x[i]-x[j];
        yb:=y[i]-y[j];
      end;
      if k=j then
      begin
        xb:=x[i+1]-x[i];
        yb:=y[i+1]-y[i];
      end;
      if k=j+1 then
      begin
        xb:=x[i+2]-x[i+1];
        yb:=y[i+2]-y[i+1];
        kk:=i;
      end;
      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);
      if k<>i 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.