* given unit-radius. Scale those points to fit the desired radius
 *)
procedure defineCircleCpts (rad : ScaledPts; centx, centy : ScaledPts;
				var CircleCpt : ControlPoints;
				var numpts : integer);
const UnitRadius = 16777216; (* TWO24 scaledpts *)
var ratio : real;
begin
  if (rad = 0) then
    begin
    complain (ERRBAD);
    writeln(logfile,'Error in fig#',pgfigurenum:0,' on page ',currpagenum:0);
    writeln(logfile,'Zero length radius for circle! Setting to 1 sp');
    rad := 1;
    end;
  ratio := float(rad) / float(UnitRadius);
  numpts := 16;
  CircleCpt[1,1] := round (ratio * 16777216.00000) + centx;
  CircleCpt[1,2] := 0 + centy; {round (ratio *      0.00000)}
  CircleCpt[2,1] := round (ratio * 15500126.47492) + centx;
  CircleCpt[2,2] := round (ratio * 6420362.60441) + centy;
  CircleCpt[3,1] := round (ratio * 11863283.20303) + centx;
  CircleCpt[3,2] := round (ratio * 11863283.20303) + centy;
  CircleCpt[4,1] := round (ratio * 6420362.60441) + centx;
  CircleCpt[4,2] := round (ratio * 15500126.47492) + centy;
  CircleCpt[5,1] := 0 + centx; {round (ratio *     -0.00000) }
  CircleCpt[5,2] := round (ratio * 16777216.00000) + centy;
  CircleCpt[6,1] := round (ratio * -6420362.60441) + centx;
  CircleCpt[6,2] := round (ratio * 15500126.47492) + centy;
  CircleCpt[7,1] := round (ratio * -11863283.20303) + centx;
  CircleCpt[7,2] := round (ratio * 11863283.20303) + centy;
  CircleCpt[8,1] := round (ratio * -15500126.47492) + centx;
  CircleCpt[8,2] := round (ratio * 6420362.60441) + centy;
  CircleCpt[9,1] := round (ratio * -16777216.00000) + centx;
  CircleCpt[9,2] := 0 + centy; {round (ratio *     -0.00000)}
  CircleCpt[10,1] := round (ratio * -15500126.47492) + centx;
  CircleCpt[10,2] := round (ratio * -6420362.60441) + centy;
  CircleCpt[11,1] := round (ratio * -11863283.20303) + centx;
  CircleCpt[11,2] := round (ratio * -11863283.20303) + centy;
  CircleCpt[12,1] := round (ratio * -6420362.60441) + centx;
  CircleCpt[12,2] := round (ratio * -15500126.47492) + centy;
  CircleCpt[13,1] := 0 + centx; {round (ratio *      0.00000) }
  CircleCpt[13,2] := round (ratio * -16777216.00000) + centy;
  CircleCpt[14,1] := round (ratio * 6420362.60441) + centx;
  CircleCpt[14,2] := round (ratio * -15500126.47492) + centy;
  CircleCpt[15,1] := round (ratio * 11863283.20303) + centx;
  CircleCpt[15,2] := round (ratio * -11863283.20303) + centy;
  CircleCpt[16,1] := round (ratio * 15500126.47492) + centx;
  CircleCpt[16,2] := round (ratio * -6420362.60441) + centy;
 (*   create the pre-list phantom *)
  CircleCpt[0,1] := CircleCpt[16,1];
  CircleCpt[0,2] := CircleCpt[16,2];  
end;


{---------------------------------------------------------------}
(* compute control points for an arc going from startangle to 
 * stopangle, centered at (centx, centy)
 *)
procedure definearcpts (rad : ScaledPts; centx, centy : ScaledPts;
			startang, stopang : integer;
			var cpts : ControlPoints;
			var nknots : integer);
var n : integer;
    a, b, curr, delta: real;
    i : integer;
begin
  a := startang * DEGTORAD;
  b := stopang * DEGTORAD;
  n := 16;

  if (a > b) then
   begin
    a := a - (2 * PI);
   end;

  delta := abs(b - a) / n;
 
  if (a = b) then
   begin
   complain (ERRNOTBAD);
   writeln(logfile,'Error in compute arc points:: should be a circle');
   end;
 curr := a;
 i := 1;
 while ((curr <= b)) do
   begin	 (* make arc about (centx,centy) *)
   cpts[i,1] := round (rad * cos (curr)) + centx;
   cpts[i,2] := round (rad * sin (curr)) + centy;
   i := i + 1;
   curr := curr + delta;
   end;  (* while *)

(* go one point beyond --
 *  around the arc so that we can have good smoothness
 *  for this phantom point 
 *)

 cpts[i,1] := round (rad * cos (b + delta)) + centx;
 cpts[i,2] := round (rad * sin (b + delta)) + centy;

(* and one phantom point before the list *)
 cpts[0,1] := round (rad * cos (a - delta)) + centx;
 cpts[0,2] := round (rad * sin (a - delta)) + centy;


 nknots := i-1;
end; 
			  
  

(* &&Module spline.p *)
(*
 Procedures below may make free use of the global variables
        arrayXY   [list of control points]
        pointmatrix [list of spline segments]
        knot    [list of spline knots]
        catrommtx  [matrix for Catmull-Rom splines]
        bsplmtx   [matrix for B-splines]
        lastPoint, intervals
*)


{-----------------------------------------------------}
function max (a, b: integer):integer;
begin
  if (a > b) then
    max := a
  else
    max := b;
end;

{-----------------------------------------------------}
function min (a, b: integer):integer;
begin
  if (a < b) then
    min := a
  else
    min := b;
end;

{---------------------------------------------------------------------}
(* initialize the Catmull-Rom basis matrix *)

procedure initcrmatrix;
begin
  catrommtx[1,1] := -0.5; catrommtx[1,2] := 1.5;
  catrommtx[1,3] := -1.5; catrommtx[1,4] := 0.5;
  catrommtx[2,1] := 1.0;  catrommtx[2,2] := -2.5;
  catrommtx[2,3] := 2.0;  catrommtx[2,4] := -0.5;
  catrommtx[3,1] := -0.5; catrommtx[3,2] := 0.0;
  catrommtx[3,3] := 0.5;  catrommtx[3,4] := 0.0;
  catrommtx[4,1] := 0.0;  catrommtx[4,2] := 1.0;
  catrommtx[4,3] := 0.0;  catrommtx[4,4] := 0.0;
end;

{-----------------------------------------------------}
procedure initbsplmatrix;
begin
  bsplmtx[1,1] := -1.0/6.0;     bsplmtx[1,2] := 0.5;
  bsplmtx[1,3] := -0.5;         bsplmtx[1,4] := 1.0/6.0;
  bsplmtx[2,1] := 0.5;          bsplmtx[2,2] := -1.0;
  bsplmtx[2,3] := 0.5;          bsplmtx[2,4] := 0.0;
  bsplmtx[3,1] := -0.5;         bsplmtx[3,2] := 0.0;
  bsplmtx[3,3] := 0.5;          bsplmtx[3,4] := 0.0;
  bsplmtx[4,1] := 1.0/6.0;      bsplmtx[4,2] := 2.0/3.0;
  bsplmtx[4,3] := 1.0/6.0;      bsplmtx[4,4] := 0.0;
end;

{--------------------------------------------------------}    
(* init the Cardinal Spline Matrix *)
procedure initcardmatrix;
begin
  cardmtx[1,1] := -1.0; cardmtx[1,2] := 1.0;
  cardmtx[1,3] := -1.0; cardmtx[1,4] := 1.0;
  cardmtx[2,1] := 2.0;  cardmtx[2,2] := -2.0;
  cardmtx[2,3] := 1.0;  cardmtx[2,4] := -1.0;
  cardmtx[3,1] := -1.0; cardmtx[3,2] := 0.0;
  cardmtx[3,3] := 1.0;  cardmtx[3,4] := 0.0;
  cardmtx[4,1] := 0.0;  cardmtx[4,2] := 1.0;
  cardmtx[4,3] := 0.0;  cardmtx[4,4] := 0.0;
end;

{--------------------------------------------------------}    
procedure initallspline;
  begin
  initcrmatrix;
  initbsplmatrix;
  initcardmatrix;
  end;


{-----------------------------------------------------}
procedure matXvector (var m: Fourby4Matrix; (* IN *)
			var v: Oneby4Vector; (* IN *)
                        var result: Oneby4Vector); (* OUT *)
var t: Oneby4Vector;
begin
  t[1] := v[1]*m[1,1] + v[2]*m[1,2] + v[3]*m[1,3] + v[4]*m[1,4];
  t[2] := v[1]*m[2,1] + v[2]*m[2,2] + v[3]*m[2,3] + v[4]*m[2,4];
  t[3] := v[1]*m[3,1] + v[2]*m[3,2] + v[3]*m[3,3] + v[4]*m[3,4];
  t[4] := v[1]*m[4,1] + v[2]*m[4,2] + v[3]*m[4,3] + v[4]*m[4,4];
  result[1] := t[1]; result[2] := t[2];
  result[3] := t[3]; result[4] := t[4];
end;

{-----------------------------------------------------}
(* actually the dot-product *)
function vecXvec (var v1, v2: Oneby4Vector) : real;
begin
  vecXvec := v1[1]*v2[1] + v1[2]*v2[2] + v1[3]*v2[3] + v1[4]*v2[4];
end;


{------------------------------------------------------}
(* basXctl is the pre-computed BasisMatrix times the control-point vector *)

function splinePosition (var basXctl : Oneby4Vector; (* IN *)
			t : real ) : real;
var tvect : Oneby4Vector;	{ vector of t values for spline matrix}
begin
  tvect[4] := 1.0;
  tvect[3] := t;
  tvect[2] := t * t;
  if (tvect[2] <= MINREAL) then
    begin			(* avoid underflow problems *)
    tvect[2] := 0.0;
    end;
  tvect[1] := t * tvect[2];  (* t^3 *)
  splinePosition := vecXvec (tvect, basXctl);  
end;  
			
{-------------------------------------------------}
function TwoToThe (n : integer) : integer;
label 78;
var i : integer;
    tmp : integer;
begin
tmp := 1;
if (n <= 0) then
  goto 78;
if (n < 6) then
  begin
    case n of
      1 : tmp := 2;
      2 : tmp := 4;
      3 : tmp := 8;
      4 : tmp := 16;
      5 : tmp := 32;
    end; (* case *)
  end  (* if *)
else
  begin
  tmp := 32;
  for i := 6 to n do
   tmp := tmp * 2;
  end;
78:
  TwoToThe := tmp;
end;  

{------------------------------------------------------}
function distance (x0, y0, x1, y1 : real) : real;
var res : real;
begin
  res := sqrt ( (x1 - x0)*(x1 - x0) + (y1 - y0)*(y1 - y0));
  distance := res;
end;  


{------------------------------------------------------}
(* compute the number of subdivisions for this span.
   We do this by a quadrature method and a simple linear-distance
   metric. This is not optimal in the number of subdivisions actually
   required, but is computationally efficient and accurate to the 
   nearest power of 2 .
   *)
function numsubdivisions (var XtimesBas, YtimesBas : Oneby4Vector;
			  resolution : ScaledPts): integer;
var n : integer;
    d : integer;  
    t : real;
    x0, y0, xt, yt : real;
begin
  x0 := splinePosition (XtimesBas, 0.0);
  y0 := splinePosition (YtimesBas, 0.0);

  t := 1.0;
  n := 0;
  xt := splinePosition (XtimesBas, t);
  yt := splinePosition (YtimesBas, t);  

  while ((round (distance (x0, y0, xt, yt)) > resolution) or
  	 (n < 1)) do
    begin
    t := t / 2.0; (* perform the quadrature *)
    n := n + 1;
    xt := splinePosition (XtimesBas, t);
    yt := splinePosition (YtimesBas, t);  
    end;  (* while *)
  numsubdivisions := TwoToThe (n);  
end;  

{------------------------------------------------------------------------}
(*  compute new control vertices such that the resulting spline
 * will interpolate through the old control points.
 * This will work as long as the actual arc length
 * between consecutive nodes does not vary from span to span.
 * The method is noted in Wu and Abel's CACM 20(10) Oct 77 paper 
 * but the actual working method is from
 *	Barsky and Greenberg's paper in
 *	CG&IP 14(3) Nov 1980 pp.203-226
 *)
procedure invertsplvertices (numpts : integer; 
				isclosed : boolean;
				var xys : ControlPoints); (* INOUT *)
var i : integer;
    beta, Xrprime, Yrprime : array[0..MAXCTLPTS] of real;
    tempxys : ControlPoints;
begin
	(* compute the values of beta *)
  beta[1] := 0.25;
  for i := 2 to numpts + 1 do
    beta[i] := 1.0 / (4.0 - beta[i - 1]);

	(* and the r primes from the original vertices *)
  Xrprime[1] := beta[1] * xys[1,1] * 5.0;
  Yrprime[1] := beta[1] * xys[1,2] * 5.0;
  for i := 2 to numpts -1 do
    begin
    Xrprime[i] := beta[i] * (6.0 * xys[i,1] - Xrprime[i - 1]);
    Yrprime[i] := beta[i] * (6.0 * xys[i,2] - Yrprime[i - 1]);
    end;  (* for *)
  Xrprime[numpts] := beta[numpts] * (5.0 * xys[numpts,1] - Xrprime[numpts - 1]);
  Yrprime[numpts] := beta[numpts] * (5.0 * xys[numpts,2] - Yrprime[numpts - 1]);

(* Now perform the back-substitution from the bottom up *)
  tempxys[numpts,1] := round (Xrprime[numpts]);
  tempxys[numpts,2] := round (Yrprime[numpts]);
  for i := numpts - 1 downto 1 do
    begin
    tempxys[i,1] := round (Xrprime[i] - beta[i] * tempxys[i + 1, 1]);
    tempxys[i,2] := round (Yrprime[i] - beta[i] * tempxys[i + 1, 2]);
    end;

if (isclosed) then
  begin
 (* at this point, we've probably been through one control-point
  *  adjustment, so let's not muck it up 
  *)
  tempxys[numpts+1,1] := tempxys[1,1];
  tempxys[numpts+1,2] := tempxys[1,2];
  tempxys[numpts+2,1] := tempxys[2,1];
  tempxys[numpts+2,2] := tempxys[2,2];
  tempxys[0,1] := tempxys[numpts,1];
  tempxys[0,2] := tempxys[numpts,2];
	  (* copy them back *)
  for i := 0 to (numpts+2) do
    begin
    xys[i,1] := tempxys[i,1];
    xys[i,2] := tempxys[i,2];
    end;  
  end  (* closed *)
else
  begin
  (* copy back *)
  for i := 2 to numpts -1 do
   begin
    xys[i,1] := tempxys[i,1];
    xys[i,2] := tempxys[i,2];
   end;
  end;  (* open*)
end; 
			      

{-----------------------------------------------------}
(*  adjust the list of control points so that we can use
 *   it for  B-spline interpolation.  
 *  Add any "phantom" vertices necessary so that the end
 *   conditions will be correct for interpolation
 *)
procedure Bctlptadjust (isclosed : boolean; isarc : boolean;
			 var n: integer; (* INOUT *)
                         var xys: ControlPoints; (* INOUT *)
                         var thx: ThickAryType); (* INOUT *)
var j : integer;
    tmp : ControlPoints;
    tmpthx : ThickAryType;
begin   (* ctlpt adjust*)

if (isclosed) then
  begin
(* here, we have to supply the last 'real' point for the user,
   and add three phantoms-- one before, and two after *)

  if (n = 2) then
    begin
    complain (ERRBAD);
    writeln(logfile,'A closed spline requires more than 2 control points ');
    writeln(logfile,'making a temporary fix in order to continue...');
    xys[3,1] := xys[1,1];
    xys[3,2] := xys[1,2];
    end;  

  for j := 1 to (n) do
    begin
    tmp[j, 1] := xys[j, 1];
    tmp[j, 2] := xys[j, 2];
    tmpthx[j] := thx[j];
    end;
		(* Now take care of the 'phantom' vertices *)    
  tmp[n+1, 1] := xys[1, 1];
  tmp[n+1, 2] := xys[1, 2];
  tmpthx[n+1] := thx[1];
  tmp[n+2, 1] := xys[2, 1];
  tmp[n+2, 2] := xys[2, 2];
  tmpthx[n+2] := thx[2];
  tmp[n+3, 1] := xys[3, 1]; 
  tmp[n+3, 2] := xys[3, 2];
  tmpthx[n+3] := thx[3];

  if (not isarc) then
    begin
    tmp[0,1] := xys[n, 1]; (* wrap around to the real last point *)
    tmp[0,2] := xys[n, 2];
    tmpthx[0] := thx[n];
    end
  else
    begin
    tmp[0,1] := xys[0,1];
    tmp[0,2] := xys[0,2];
    tmpthx[0] := thx[0];
    end;

  n := n + 1; 	(* we supplied the 'last' point for the user *)

  for j := 0 to n+2 do
    begin
    xys[j,1] := tmp[j,1];
    xys[j,2] := tmp[j,2];
    thx[j] := tmpthx[j];
    end;  (* for *)
  end  (* if closed *)
else 
  begin		 (* OPEN SPLINE *)
  if (not isarc) then
    begin
    tmp[0,1] := 2 * xys[1, 1] - xys[2,1];
    tmp[0,2] := 2 * xys[1, 2] - xys[2,2];
    end
  else
    begin
    tmp[0,1] := xys[0,1];
    tmp[0,2] := xys[0,2];
    end;
  tmpthx[0] := thx[1];

  for j := 1 to (n) do
    begin
    tmp[j, 1] := xys[j, 1];
    tmp[j, 2] := xys[j, 2];
    tmpthx[j] := thx[j];
    end;
  
  tmp[n+1, 1] := 2 * xys[n, 1] - xys[n-1,1];
  tmp[n+1, 2] := 2 * xys[n, 2] - xys[n-1,2];
  tmpthx[n+1] := thx[n];

  tmp[n+2, 1] := tmp[n+1, 1];
  tmp[n+2, 2] := tmp[n+1, 2];
  tmpthx[n+2] := thx[n];

  for j := 0 to n+2 do
    begin
    xys[j,1] := tmp[j,1];
    xys[j,2] := tmp[j,2];
    thx[j] := tmpthx[j];
    end;  (* for *)
  end; (*  if open *)
  
end;



{-----------------------------------------------------}
(*  adjust the list of control points so that we can use
 *       it for simple Catmull-Rom spline interpolation.  
 *  Add any "phantom" vertices necessary so that the end
 *   conditions will be correct for interpolation
 *)
procedure CRctlptadjust (isclosed : boolean; isarc : boolean;
			 var n: integer; (* INOUT *)
                         var xys: ControlPoints; (* INOUT *)
                         var thx: ThickAryType); (* INOUT *)
var j : integer;
    tmp : ControlPoints;
    tmpthx : ThickAryType;
begin   (* ctlpt adjust*)
if (isclosed) then
  begin
(* here, we have to supply the last 'real' point for the user,
   and add three phantoms-- one before, and two after *)

  if (n = 2) then
    begin
      complain (ERRBAD);
      writeln(logfile,'A closed spline requires more than 2 control points ');
      writeln(logfile,'making a temporary fix in order to continue...');
      xys[3,1] := xys[1,1];
      xys[3,2] := xys[1,2];
    end;  


  for j := 1 to (n) do
    begin
    tmp[j, 1] := xys[j, 1];
    tmp[j, 2] := xys[j, 2];
    tmpthx[j] := thx[j];
    end;
			(* the phantom vertices *)    
    tmp[n+1, 1] := xys[1, 1];
    tmp[n+1, 2] := xys[1, 2];
    tmpthx[n+1] := thx[1];
    tmp[n+2, 1] := xys[2, 1];
    tmp[n+2, 2] := xys[2, 2];
    tmpthx[n+2] := thx[2];
    tmp[n+3, 1] := xys[3, 1];
    tmp[n+3, 2] := xys[3, 2];
    tmpthx[n+3] := thx[3];
  
    if (not isarc) then
      begin
      tmp[0,1] := xys[n, 1]; (* wrap around to the real last point *)
      tmp[0,2] := xys[n, 2];
      tmpthx[0] := thx[n];
      end
    else
      begin
      tmp[0,1] := xys[0,1];
      tmp[0,2] := xys[0,2];
      tmpthx[0] := thx[0];
      end;
    n := n + 1; (* we supplied the 'last' point for the user *)
  
    for j := 0 to n+2 do
      begin
      xys[j,1] := tmp[j,1];
      xys[j,2] := tmp[j,2];
      thx[j] := tmpthx[j];
      end;  (* for *)
  end  (* if closed *)
else
  begin (* OPEN SPLINE *)
  if (not isarc) then
    begin
    tmp[0,1] := xys[1, 1]; (* double the first point *)
    tmp[0,2] := xys[1, 2];
    end
  else
    begin
    tmp[0,1] := xys[0,1];
    tmp[0,2] := xys[0,2];
    end;  
  tmpthx[0] := thx[1];

  for j := 1 to (n) do
    begin
    tmp[j, 1] := xys[j, 1];
    tmp[j, 2] := xys[j, 2];
    tmpthx[j] := thx[j];
    end;
    
  tmp[n+1, 1] := xys[n, 1]; (* and triple the last *)
  tmp[n+1, 2] := xys[n, 2];
  tmpthx[n+1] := thx[n];
  tmp[n+2, 1] := xys[n, 1];
  tmp[n+2, 2] := xys[n, 2];
  tmpthx[n+2] := thx[n];

  for j := 0 to n+2 do
    begin
    xys[j,1] := tmp[j,1];
    xys[j,2] := tmp[j,2];
    thx[j] := tmpthx[j];
    end;  (* for *)
  end; (* if open *)
end;    (* ctlpt adjust *)

     

{----------------------------------------------------------}

procedure interpsplines (splinetype: SplineKind;
			 isclosed: boolean;
			 isanArc: boolean;
			 linepatt : LineStyle;
                         var basismatrix : Fourby4Matrix; (* IN *)
                         numctls: integer; 
                         var arrayXY: ControlPoints; (* IN *)
                         var pointmatrix: SplineSegments; (* OUT *)
                         varythicks: boolean;
                         var thickmatrix: ThickAryType; (* IN *)
                         var TTmatrix: ThickAryType); (* OUT *)
label 32;
var xctl, yctl,		{ vectors of x, y posits of control points}
    wctl : Oneby4Vector; {vector of thicknesses at each ctl pt}
    t, incr: real;
    Pi: integer;	{ P sub i }
    i, currpt : integer;    
    theresolution : ScaledPts;

begin (* interp splines*)
  if ((not isclosed) and (isanArc)) then
    numctls := numctls + 1; (* lie a little *)

   case (splinetype) of

     BSPL: Bctlptadjust (isclosed, isanArc, numctls, arrayXY, thickmatrix);
     
     CARD,
     CATROM:  CRctlptadjust (isclosed, isanArc, numctls, arrayXY, thickmatrix);
    
     INTBSPL: begin
     		if (isclosed) then
		  begin
		  Bctlptadjust (true, isanArc, numctls, arrayXY, thickmatrix);
		  invertsplvertices (numctls, true, arrayXY);
		  end
		else 
		  begin
		  invertsplvertices (numctls, false, arrayXY);
		  Bctlptadjust (false, isanArc, numctls, arrayXY, thickmatrix);
		  end;  (* else *)
     	      end; (* Interpolating Bsplines *)
   end;

  if ((not isclosed) and (isanArc)) then
    numctls := numctls - 1; (* UN-lie a little *)


(* this is the scheme:
 *    val :=  t-vector   *  Basis matrix     * point matrix
 *        [t^3  t^2 t 1] *      [[Ms]]       * [Pi-1 Pi Pi+1 Pi+2]
 *    where "Pi-1" is "P sub (i-1)", etc.
 *
 *  But we do this in a round about way:
 *        Point matrix * basis
 *   then   * t-vector   will yield the single value
 *   
 *   there are certainly faster ways to do this, 
 *   but this is the easiest to understand
 *)

  currpt := 1;
  case linepatt of
     solid : theresolution := MAXVECLENsp;
     dotted,
     dashed,
     dotdash : theresolution := 3 * MAXVECLENsp; {###}
   end;

  for Pi := 1 to (numctls - 1) do
    begin
    xctl[1] := float(arrayXY[Pi-1, 1]);
    xctl[2] := float(arrayXY[Pi,   1]);
    xctl[3] := float(arrayXY[Pi+1, 1]);
    xctl[4] := float(arrayXY[Pi+2, 1]);
    yctl[1] := float(arrayXY[Pi-1, 2]);
    yctl[2] := float(arrayXY[Pi,   2]);
    yctl[3] := float(arrayXY[Pi+1, 2]);
    yctl[4] := float(arrayXY[Pi+2, 2]);
    matXvector (basismatrix, xctl, xctl);
    matXvector (basismatrix, yctl, yctl);

	(* compute the delta-t increment for this segment
		based on a metric for subdivision *)
    intervals := numsubdivisions (xctl, yctl, theresolution);
    if ((linepatt = solid) and (intervals <= 2)) then
      intervals := intervals * 2;
    incr := 1.0 / intervals;

	(* avoid over-flowing the "pointmatrix" *)
    if ((currpt + intervals - 1) >= MAXSPLINESEGS) then
       begin
       complain (ERRREALBAD);
       writeln (logfile,'error: Too many spline segments required.');
       writeln (logfile,' Reducing the number of control points to get output.');
       goto 32;
       end;
  
    t := 0.0;
    while (t < 0.999999999) do
      begin
	pointmatrix[currpt, 1] := round (splinePosition (xctl, t));
	pointmatrix[currpt, 2] := round (splinePosition (yctl, t));

	if (varythicks) then
	  begin
	    wctl[1] := float(thickmatrix[Pi-1]);
	    wctl[2] := float(thickmatrix[Pi  ]);
	    wctl[3] := float(thickmatrix[Pi+1]);
	    wctl[4] := float(thickmatrix[Pi+2]);
	    matXvector (catrommtx, wctl, wctl);  (* requires using Catmull-Rom *)
	    TTmatrix[currpt] := round (splinePosition (wctl, t));
	  end;
    
        t := t + incr;
        currpt := currpt + 1;
      end; (* while loop *)


    end; (* for loop *)

32:
	(* the END-condtion *)
    pointmatrix[currpt, 1] := round (splinePosition (xctl, 1.0));
    pointmatrix[currpt, 2] := round (splinePosition (yctl, 1.0));    
    if (varythicks) then
      begin
	wctl[1] := thickmatrix[numctls-2];
	wctl[2] := thickmatrix[numctls-1];
	wctl[3] := thickmatrix[numctls];
	wctl[4] := thickmatrix[numctls+1];
	matXvector (catrommtx, wctl, wctl);  (* requires using Catmull-Rom *)
	TTmatrix[currpt] := round (splinePosition (wctl, 1.0));
      end;

    lastPoint := currpt;

end; (* interpsplines *)


{----------------------------------------------------------------}
procedure drawSpline (splinetype : SplineKind;
		     isclosed: boolean;
		     isanArc: boolean;
		     patt : LineStyle;
                     numctls: integer;
                     var arrayXY: ControlPoints; (* IN *)
                     var pointmatrix: SplineSegments; (* OUT *)
                     varythicks: boolean;
                     var thickmatrix: ThickAryType; (* IN *)
                     var TTmatrix: ThickAryType); (* OUT *)
begin
  lastPoint := 0;


  case (splinetype) of
    CATROM : interpsplines (splinetype, isclosed, isanArc, patt, catrommtx,
			   numctls, arrayXY, pointmatrix,
              		   varythicks, thickmatrix, TTmatrix);

    CARD : interpsplines (splinetype, isclosed, isanArc, patt, cardmtx, 
			   numctls, arrayXY, pointmatrix, 
	                   varythicks, thickmatrix, TTmatrix);

    BSPL : interpsplines (splinetype, isclosed, isanArc, patt, bsplmtx, 
			   numctls, arrayXY, pointmatrix, 
	                   varythicks, thickmatrix, TTmatrix);

    INTBSPL : interpsplines (splinetype, isclosed, isanArc, patt, bsplmtx,
			   numctls, arrayXY, pointmatrix, 
	                   varythicks, thickmatrix, TTmatrix);
  end; (*Case *)	      	     
end;


(* &&module TeXtyl *)
{----------------------------------------------------------------}
(* rotate a (x,y) point about mx, my *)
procedure ptrotate (var x, y : integer;
                        mx, my: integer;
                        angle : real);
var tmpx, tmpy : integer;
    cosa, sina : real;
begin
  tmpx := x - mx;       
  tmpy := y - my;
  cosa := cos(angle * DEGTORAD); 
  sina := sin(angle * DEGTORAD);
  x := round(tmpx * cosa - tmpy * sina) + mx;
  y := round(tmpx * sina + tmpy * cosa) + my;
end;

{----------------------------------------------------------------}
(* transform two line points: scale, rotate and translate 
*)
procedure xfmlinepts (var x1, y1, x2, y2 : ScaledPts;
                        offh, offv : ScaledPts;
                        midx, midy : ScaledPts;
                        scalefact : real;
                        theta : real;
                        dx, dy : ScaledPts;
                        sx, sy : real);
begin
  if ((sx = 0.0) or (sy = 0.0)) then
    begin
      complain (ERRBAD);
      writeln(logfile,'?? Some scale factor is Zero... continuing anyway');
    end;
        (* scale about center of item*)
  if ((sx <> 1.0) or (sy <> 1.0)) then
   begin
   x1 := round((x1 - midx) * sx) + midx;
   x2 := round((x2 - midx) * sx) + midx;
   y1 := round((y1 - midy) * sy) + midy;     
   y2 := round((y2 - midy) * sy) + midy;
   end;
      (* rotate if necessary *)
   if (theta <> 0.0) then
     begin  (* rotate about the midpoint *)
     ptrotate(x1, y1, midx, midy, theta);
     ptrotate(x2, y2, midx, midy, theta);
     end;
      (* translate *)
   x1 := (x1 + round(dx * scalefact) + offh);
   x2 := (x2 + round(dx * scalefact) + offh);
   y1 := (y1 + round(dy * scalefact) + offv);
   y2 := (y2 + round(dy * scalefact) + offv);
end;  (* xfmlinepts *)

{----------------------------------------------------------------}
procedure xfmcontpts (var xpts : ControlPoints; xknots : integer;
                        offh, offv : ScaledPts; midx, midy : ScaledPts;
                        scalefact : real;
                        theta : real; dx, dy : ScaledPts; sx, sy : real);
var i : integer;
begin
    (* scale about center of item *)
 if ((sx <> 1.0) or (sy <> 1.0)) then
  for i := 0 to xknots do
     begin
     xpts[i,1] := round((xpts[i,1] - midx) * sx) + midx;
     xpts[i,2] := round((xpts[i,2] - midy) * sy) + midy;
     end;

  if (theta <> 0.0) then
    begin (* rotate about center *)
    for i := 0 to xknots do
      begin
      ptrotate (xpts[i,1], xpts[i,2], midx, midy, theta);
      end;
    end;
    (* translate *)
  for i := 0 to xknots do
    begin
    xpts[i,1] := (xpts[i,1] + round(dx * scalefact) + offh);
    xpts[i,2] := (xpts[i,2] + round(dy * scalefact) + offv);
    end;
end;  (* xfmcontpts *)


{----------------------------------------------------------------}
(* convert into DVI space and offset by H & V *)
procedure dvilinepts (var x1, y1, x2, y2 : ScaledPts;
			offh, offv : ScaledPts);
begin
   x1 := (x1  + offh);
   x2 := (x2  + offh);
   y1 := (y1 * (-1) + offv);
   y2 := (y2 * (-1) + offv);
end;

{----------------------------------------------------------------}
(* convert into DVI space and offset by H & V *)
procedure dvicontpts (var xpts : ControlPoints; xknots : integer;
                        offh, offv : ScaledPts);
var i : integer;
begin
  for i := 0 to xknots do
    begin
    xpts[i,1] := (xpts[i,1]  + offh);
    xpts[i,2] := (xpts[i,2] * (-1) + offv);
    end;
end;

{----------------------------------------------------------------}
(*	transform all the figure's elements according to the 
	top-level tranformation requirements in 1st Quadrant space.
	then reset the toplevel's xfms.
*)
procedure toplevelxfm (toplev, curfig : pItem; recurlevel : integer);
var pi : pItem;
    null1, null2 : ScaledPts;
    old1, old2 : ScaledPts;
    midx, midy : ScaledPts;
begin
  with toplev^ do
    begin
    midy := (BBty - BBby) div 2;
    midx := (BBrx - BBlx) div 2;
    end;
  pi := curfig^.body^.things;  { if recur==0, this is same as toplev }
  while (pi <> nil) do
    begin
    with pi^ do
      begin
      case (kind) of
	Aline : begin
		xfmlinepts (lx1, ly1, lx2, ly2, 0, 0, midx, midy, 1.0, 
		      toplev^.figtheta, toplev^.fdx, toplev^.fdy,
		      toplev^.fsx, toplev^.fsy);
		end;
	Aspline : begin
		  xfmcontpts (spts, nsplknots, 0, 0, midx, midy, 1.0,
		      toplev^.figtheta, toplev^.fdx, toplev^.fdy,
		      toplev^.fsx, toplev^.fsy);
		  end;
	Attspline : begin
		  xfmcontpts (ttpts, nttknots, 0, 0, midx, midy, 1.0,
		      toplev^.figtheta, toplev^.fdx, toplev^.fdy,
		      toplev^.fsx, toplev^.fsy);
		    end;
	Aarc : begin
	       null1 := 0; null2 := 0;
	       old1 := acentx; old2 := acenty;
	       xfmlinepts (acentx, acenty, null1, null2, 0,0, midx, midy, 1.0,
		    toplev^.figtheta, toplev^.fdx, toplev^.fdy,
		    toplev^.fsx, toplev^.fsy);		
			      
	       xfmcontpts (arcpts, narcknots + 1, 0, 0, old1, old2, 1.0,
		      toplev^.figtheta, 
		      toplev^.fdx + (acentx - old1), 
		      toplev^.fdy + (acenty - old2),
		      toplev^.fsx, toplev^.fsy);
	       end;		  	
	Alabel : begin
		 null1 := 0; null2 := 0;
		 xfmlinepts (labx, laby, null1, null2, 0, 0, midx, midy, 1.0,
		      toplev^.figtheta, toplev^.fdx, toplev^.fdy,
		      toplev^.fsx, toplev^.fsy);		
		 end;
	Abeam : ;   (* not transformable *)

	Atieslur: ; (* not transformable *)

	Afigure : begin
		    toplevelxfm (toplev, pi, recurlevel + 1);
		  end;
      end; (* case *)
    end; (* with *)
    pi := pi^.nextitem;
    end;  (* while *)
  if (recurlevel = 0) then
    begin (* reset the toplevel's xfms *)
    with toplev^ do
      begin
      figtheta := 0.0;
      fsx := 1.0; fsy := 1.0;
      fdx := 0;   fdy := 0;
      end;    
    end;
end;


{----------------------------------------------------------------}
function scalefitfactor (actualwid, actualht, 
			 goalwid, goalht: ScaledPts): real;
var sx, sy : real;
begin
  sx := goalwid/actualwid;
  sy := goalht/actualht;
  if (sx < sy) then
    scalefitfactor := sx
  else
    scalefitfactor := sy;
end;  



(* ---- The handlers for each primitive ---- 
 *   The result of calling each handler is either immediate
 *       output to the buffer of the commands to produce the
 *       primitive, OR the primitive gets pushed onto a stack/list
 *       that defines a current 'figure' (set of prims) for
 *       output at a later time
 *
 *  Look at linehandle for a basic idea of how the handlers
 *  work. the others follow pretty closely.
 *)


{------------------------------------------------------------}
procedure linehandle (figdepth : integer; scalefact: real; 
                     x1, y1, x2, y2 : ScaledPts;
                     dvih, dviv : ScaledPts; (* possible dvi-offsets *)
                     thk : VThickness; vk : VectKind;
		     patt : LineStyle;
		     minx, maxx, miny, maxy : ScaledPts;
                     tx, ty: ScaledPts; sx, sy, r : real);
var midx, midy : ScaledPts;                  
    lineitem : pItem;
begin
   midx := (minx + maxx) div 2;
   midy := (miny + maxy) div 2;

	(* do local primitive -level transformations *)
   xfmlinepts (x1, y1, x2, y2, dvih, dviv,
                midx, midy, scalefact, r, tx, ty, sx, sy);