program ppbspln (input,output);


(**********************************************************************)
(*                                                                    *)
(*  record type definition selection:                                 *)
(*                                                                    *)
(*  comments:  some users may wish to reformulate the following       *)
(*             record definitions in terms of linked lists.           *)
(*             for instance, the vector components, knot sequences,   *)
(*             break points, b-spline coefficients and weights could  *)
(*             be expressed in that way.                              *)
(*                                                                    *)
(**********************************************************************)
const
  maxdeg = 5;

type 
  tuple = array [1..4] of DOUBLE

(**********************************************************************)
(*                                                                    *)
(*  record type vector                                                *)
(*                                                                    *)
(*  comments:  the basic building blcok for coefficients              *)
(*                                                                    *)
(**********************************************************************)

  vector = record
     dim : integer;    (* the dimension fo the vector *)
component : tuple;     (* the array of components *)
  end;                 (* the definition of record type vector *)

real44 = array [0..44] of DOUBLE;
vec44 = array [0.44] of vector;
integer44 = array [0..44] of integer;
real4444 = array [0..44,0..44] of DOUBLE;

real22 = array [0..22] of DOUBLE;
vec22 = array [0.22] of vector;
integer22 = array [0..22] of integer;
real2222 = array [0..22,0..22] of DOUBLE;
vec2222 = array [0..22,0..22] fo vector;

vecmd = array [0..maxdeg] of vector;
realmd = array [0..maxdeg] of DOUBLE;
vecmdmd = array [0..maxdeg,0..maxdeg] of vector;
realmdmd = array [0..maxdeg,0..maxdeg] of DOUBLE;

(**********************************************************************)
(*                                                                    *)
(*  record type; poly                                                 *)
(*                                                                    *)
(*  comments:  vector-valued polynomial (power basis form)            *)
(*             the basic building block for the pp and                *)
(*             rational_poly record types (see below)                 *)
(*                                                                    *)
(**********************************************************************)
poly = record
degree : integer;    (* degree *)
      dim : integer; (* dimension of the coefficients *)
    coeff : vecmd;   (* coeff(.i.) is the coefficient of t**i *)
end;
(* definition of record type poly *)


(**********************************************************************)
(*  record type:  pp              *)
(*                                                                    *)
(*     comments:  piecewise polynomial.                               *)
(*                                                                    *)
(**********************************************************************)
pp = record       
   degree   : integer;        (* degree *)
   num_segs : integer;        (* number of segments *)
   dim      : integer;        (* dimension of coefficients *)
   break_point : real22;      (* break points *)
   segment : array(.1..22.) of poly;         (* the polynomials in pp *)
end;                (* definition of record type pp *)


(**********************************************************************)
(*  record type:  bspline         *)
(*                                *)
(*     comments:  bspline curve   *)
(**********************************************************************)
bspline = record  
   degree : integer;                      (* degree *)
   dim : integer;  (* dimension of the coefficients *)
   number_of_knots : integer;    (* number of knots *)
   knot : real44;                  (* knot sequence *)
   coeff : vec44;       (* the bspline coefficients *)
end;           (* definition of record type bspline *)


(**********************************************************************)
(*              *)
(*  record type:  rational_poly   *)
(*                                *)
(*     comments:  rational vector-valued polynomial (power basis form)*)
(*                the basic building block for record type            *)
(*                rational_pp (see below)           *)
(**********************************************************************)
rational_poly = record              
   degree : integer;    (* degree *)
   dim : integer;  (* dimension of the coefficients *)
   top_coeff : vecmd;     (* top_coeff(i) is the coeff of t**i in num *)
   bot_coeff : realmd;  (* bot_coeff(i) is the coeff of t**i in denom *)
   numerator : poly;         (* redundant info stored for convenience *)
   denominator : poly;       (* redundant info stored for convenience *)
end;     (* definition of record type rational_poly *)


(**********************************************************************)
(*  record type:  rational_pp     *)
(*                                *)
(*     comments:  rational piecewise polynomial.    *)
(*                generalization of iges parametric spline curve*)
(*                entity in the following ways:     *)
(*                (1)  degree is arbitrary          *)
(*                (2)  dimension is arbitrary       *)
(*                (3)  form is rational             *)
(**********************************************************************)
rational_pp = record                
      degree : integer; (* degree *)
      num_segs : integer;   (* number of rational polynomial segments *)
      dim : integer;                 (* dimension of the coefficients *)
      break_point : real22;         (* break points *)
      rpp_segment : array(.1..22.) of rational_poly;(* the components *)
end;       (* definition of record type rational_pp *)
                  
                  
(**********************************************************************)

(*  record type:  rational_bspline                  *)

(*     comments:  rational bspline curve.           *)
(*                generalization of iges rational b-spline curve      *)
(*                entity, in that dimension of coefficients is        *)
(*                arbitrary.      *)
2
(**********************************************************************)
rational_bspline = record           
   degree : integer;    (* degree *)
   dim : integer;  (* dimension of the coefficients *)
   number_of_knots : integer;     (* # knots, counting multiplicities *)
   knot : real44;            (* set of knots including multiplicities *)
   coeff : vec44;      (* the b-spline coefficients *)
   weight : real44;                  (* the weights *)
   numerator : bspline;      (* redundant info stored for convenience *)
   denominator : bspline;    (* redundant info stored for convenience *)
end;  (* definition of record type rational_bspline *)
                  
                  
(**********************************************************************)

(*  record type:  bsurf           *)

(*     comments:  b-spline surface.                 *)

(**********************************************************************)
bsurf = record    
   degree1 : integer;        (* degree in 1st param *)
   degree2 : integer;        (* degree in 2nd param *)
   dim : integer;      (* dimension of coefficients *)
   num_of_knots1 : integer;  (* #knots in 1st param *)
   num_of_knots2 : integer;  (* #knots in 2nd param *)
   knot1 : real44;    (* knot sequence in 1st param *)
   knot2 : real44;    (* knot sequence in 2nd param *)
   coeff : vec2222;    (* the b-spline coefficients *)
end;             (* definition of record type bsurf *)
                  
                  
(**********************************************************************)

(*  record type:  rbsurf          *)

(*     comments:  rational b-spline surface.        *)
(*                a generalization of the iges rational b-spline      *)
(*                surface entity in that the dimension of the         *)
(*                coefficients is arbitrary.        *)

(**********************************************************************)
rbsurf = record   
   degree1 : integer;        (* degree in 1st param *)
   degree2 : integer;        (* degree in 2nd param *)
   dim : integer;      (* dimension of coefficients *)
   num_of_knots1 : integer; (* # knots in 1st param *)
   num_of_knots2 : integer; (* # knots in 2nd param *)
   knot1 : real44;    (* knot sequence in 1st param *)
   knot2 : real44;    (* knot sequence in 2nd param *)
   coeff : vec2222;    (* the b-spline coefficients *)
   weight : real2222;                (* the weights *)
   numerator : bsurf;        (* redundant info stored for convenience *)
   denominator : bsurf;      (* redundant info stored for convenience *)
end;            (* definition of record type rbsurf *)
                  
                  
(**********************************************************************)

(*  record type:  psurf           *)

(*     comments:  polynomial surface entity (power basis)             *)
(*                the basic building block for record types           *)
(*                ppsurf and rpsurf (see below).    *)

(**********************************************************************)
psurf = record    
   degree1 : integer;        (* degree in 1st param *)
   degree2 : integer;        (* degree in 2nd param *)
   dim : integer;  (* dimension of the coefficients *)
   coeff :vecmdmd;        (* coeff(i,j) is the coeff of (s**i)*(t**j) *)
end;             (* definition of record type psurf *)
                  
                  
(**********************************************************************)

(*  record type:  ppsurf          *)

(*     comments:  piecewise polynomial surface      *)

(**********************************************************************)
ppsurf = record   
   num_segs1 : integer;      (* number of patches along 1st parameter *)
   num_segs2 : integer;      (* number of patches along 2nd parameter *)
   degree1 : integer;    (* degree in 1st parameter *)
   degree2 : integer;    (* degree in 2nd parameter *)
   dim : integer;      (* dimension of coefficients *)
   break_point1 : real22;             (* breakpoints in 1st parameter *)
   break_point2 : real22;             (* breakpoints in 2nd parameter *)
   ppatch : array(.1..4,1..4.) of psurf;      (* the array of patches *)
end;            (* definition of record type ppsurf *)
                  
                  
(**********************************************************************)

(*  record type:  rpsurf          *)

(*     comments:  rational polynomial surface.  the basic building    *)
(*                block for record type rppsurf (see below).          *)

(**********************************************************************)
rpsurf = record   
   degree1 : integer;    (* degree in 1st parameter *)
   degree2 : integer;    (* degree in 2nd parameter *)
   dim : integer;          (* dimension of the numerator coefficients *)
   top_coeff :vecmdmd;(* top_coeff(i,j) = coe of (s**i)*(t**j) in num *)
   bot_coeff : realmdmd;(* bot_coeff(i,j) " " "    "      "  in denom *)
   numerator : psurf;        (* redundant info stored for convenience *)
   denominator : psurf;      (* redundant info stored for convenience *)
end;            (* definition of record type rpsurf *)
                  
                  
(**********************************************************************)

(*  record type:  rppsurf         *)

(*     comments:  rational piecewise polynomial surface.              *)
(*                a generalization of the iges parametric spline      *)
(*                surface entity in the following ways:               *)
(*                (1)  degree is arbitrary          *)
(*                (2)  dimension is arbitrary       *)
(*                (3)  form is rational             *)

(**********************************************************************)
rppsurf = record  
   num_segs1 : integer;      (* number of patches along 1st parameter *)
   num_segs2 : integer;      (* number of patches along 2nd parameter *)
   degree1 : integer;    (* degree in 1st parameter *)
   degree2 : integer;    (* degree in 2nd parameter *)
   dim : integer;  (* dimension of the coefficients *)
   break_point1 : real22;            (* break points in 1st parameter *)
   break_point2 : real22;            (* break points in 2nd parameter *)
   rppatch : array(.1..4,1..4.) of rpsurf;    (* the array of patches *)
end;           (* definition of record type rppsurf *)
                  
(* global variables *)              
var               
    n_knots_to_add: integer;        
    n_knots_to_remove: integer;     
    return_code : integer;          
    machine_epsilon : real;         
    idebug : integer;               
    debug : text; 
    report : text;                  
    i : integer;  
    f,f0,f1 : rational_pp;          
    g,g0 : rational_bspline;        
    b : rbsurf;   
    p : rppsurf;  
                  
(**********************************************************************)
(**********************************************************************)
(**********************************************************************)

(*      name:  writevec           *)

(*  category:  vector operations  *)

(*   purpose:  writes the value of a vector to a file                 *)

(*     input:  outfil - the name of the file to write to              *)
(*             a - the name of the vector to write  *)

(*    output:  no output parameters.                *)
(*             the value of the vector is written to the file.        *)
(*             the number of decimal places written depends upon the  *)
(*             dimension of the vector.             *)

(**********************************************************************)
                  
procedure writevec(var outfil : text;                 
 var a : vector );                  
   var            
      i : integer;                  
    f,d : integer;                  
begin (* function writevec *)       
   with a do      
   begin (* with a *)               
      d := 15;    
      if (dim > 2) then d := 6;     
      f := d + 4; 
      for i := 1 to dim do          
      begin       
         write(outfil,component(.i.):f:d);            
         write(outfil,' ');         
      end;        
   end;  (* with a *)               
end;  (* function writevec *)       
                  
(**********************************************************************)
%page             
(**********************************************************************)

(*      name:  sum                *)

(*  category:  vector operations  *)

(*   purpose:  finds the sum of two vectors         *)

(*     input:  a - a vector       *)
(*             b - a vector       *)

(*    output:  sum(a,b) = a + b   *)
(*             if the input vectors are not of the same dimension,    *)
(*             an error message is written          *)

(**********************************************************************)
                  
function sum(a:vector; b:vector):vector;              
var               
   i : integer;   
begin (* function sum *)            
   if not(a.dim = b.dim)            
   then           
begin (* error report *)            
   writeln(' in function sum dims not equal ');       
   writeln(' a =');writevec(output,a);                
   writeln(' a.dim =',a.dim:1);     
   writeln(' b =');writevec(output,b);                
   writeln(' b.dim =',b.dim:1);     
end   (* error report *)            
   else           
   begin (* normal case *)          
      sum.dim := a.dim;             
      for i := 1 to a.dim do        
      sum.component(.i.) :=         
      a.component(.i.) + b.component(.i.);            
   end;  (* normal case *)          
end;  (* function sum *)            
                  
(**********************************************************************)
(**********************************************************************)

(*      name:  difference         *)

(*  category:  vector operations  *)

(*   purpose:  finds the difference of two vectors  *)

(*     input:  a - a vector       *)
(*             b - a vector       *)

(*    output:  difference(a,b) = a - b              *)
(*             if the input vectors are not of the same dimension,    *)
(*             an error message is written          *)

(**********************************************************************)
                  
function difference(a:vector; b:vector):vector;       
var               
   i : integer;   
begin (* function difference *)     
   if not(a.dim = b.dim)            
   then           
begin (* error report *)            
   writeln(' in function difference dims not equal ');                  
   writeln(' a =');writevec(output,a);                
   writeln(' a.dim =',a.dim:1);     
   writeln(' b =');writevec(output,b);                
   writeln(' b.dim =',b.dim:1);     
end   (* error report *)            
   else           
   begin (* normal case *)          
      difference.dim := a.dim;      
      for i := 1 to a.dim do        
      difference.component(.i.) :=  
      a.component(.i.) - b.component(.i.);            
   end;  (* normal case *)          
end;  (* function difference *)     
                  
(**********************************************************************)
(**********************************************************************)

(*      name:  norm               *)

(*  category:  vector operations  *)

(*   purpose:  find the euclidean norm of a vector  *)

(*     input:  a - a vector       *)

(*    output:  norm(a) = !! a !!  *)

(**********************************************************************)
                  
function norm(a:vector) : real;     
var               
    i : integer;  
    norm_squared : real;            
   begin (* function norm *)        
      norm_squared := 0.0;          
      with a do   
      for i := 1 to dim do          
      norm_squared := norm_squared + sqr(component(.i.));               
                  
      norm := sqrt(norm_squared);   
                  
   end;  (* function norm *)        
                  
(**********************************************************************)
(**********************************************************************)

(*      name:  dist               *)

(*  category:  vector operations  *)

(*   purpose:  compute the euclidean distance between two vectors     *)

(*     input:  a - a vector       *)
(*             b - a vector       *)

(*    output:  dist(a,b) = !! (a - b) !!            *)

(**********************************************************************)
                  
function dist(a : vector ; b : vector) : real;        
begin (* function dist *)           
   dist := norm(difference(a,b));   
end;  (* function dist *)           
                  
(**********************************************************************)
(**********************************************************************)

(*      name:  product            *)

(*  category:  vector operations  *)

(*   purpose:  compute the vector which is the product of a scalar    *)
(*             and a vector       *)

(*     input:  s - a scalar       *)
(*             v - a vector       *)

(*    output:  product(s,v) = s*v *)

(**********************************************************************)
                  
function product(s : real; v : vector) : vector;      
var               
   i : integer;   
begin (* function product *)        
   product.dim := v.dim;            
   for i := 1 to v.dim do           
   product.component(.i.) :=        
   s*v.component(.i.);              
end;  (* function product *)        
                  
(**********************************************************************)
(**********************************************************************)

(*      name:  quotient           *)

(*  category:  vector operations  *)

(*   purpose:  compute the vector which is the quotient of a vector   *)
(*             and a scalar       *)

(*     input:  v - a vector       *)
(*             s - a scalar       *)

(*    output:  quotient(v,s) = v/s                  *)

(**********************************************************************)
                  
function quotient(v : vector; s : real) : vector;     
var               
   i : integer;   
begin (* function quotient *)       
   quotient.dim := v.dim;           
   for i := 1 to v.dim do           
   quotient.component(.i.) :=       
   v.component(.i.)/s;              
end;  (* function quotient *)       
(**********************************************************************)
(**********************************************************************)

(*      name:  zerovec            *)

(*  category:  vector operations  *)

(*   purpose:  creates the zero vector of a given dimension           *)

(*     input:  dimension - the desired dimension    *)

(*    output:  zerovec(dimension) = (0,...,0)       *)

(**********************************************************************)
                  
function zerovec(dimension : integer):vector;         
var               
   i : integer;   
begin (* function zerovec *)        
   zerovec.dim := dimension;        
   for i := 1 to dimension do       
   zerovec.component(.i.) := 0.0;   
end;  (* function zerovec *)        
                  
(**********************************************************************)
(**********************************************************************)

(*      name:  project            *)

(*  category:  vector operations  *)

(*   purpose:  project a vector onto a vector having one less         *)
(*             dimension.         *)

(*     input:  v - a vector. v = (v1,v2,...vn-1,vn) *)

(*    output:  project(v) = the vector (v1,v2,...vn-1)                *)

(**********************************************************************)
function project(v:vector):vector;  
   var            
      i : integer;                  
      temp : vector;                
begin (* function project *)        
   temp.dim := v.dim - 1;           
   for i := 1 to temp.dim do        
   temp.component(.i.) := v.component(.i.);           
   project := temp;                 
end;  (* function project *)        
(**********************************************************************)
(**********************************************************************)

(*      name:  append             *)

(*  category:  vector operations  *)

(*   purpose:  given a vector and a real number, create a new vector  *)
(*             of one dimension larger by using the real number as    *)
(*             the new last component.              *)

(*     input:  v -  a vector.  v = (v1,...,vn)      *)
(*             r - a real number  *)

(*    output:  append(v,r) = (v1,...,vn,r)          *)

(**********************************************************************)
function append(v:vector; r:real):vector;             
   var            
      i : integer;                  
      temp : vector;                
begin (* function append *)         
   temp.dim := v.dim + 1;           
   for i := 1 to v.dim do           
   temp.component(.i.) := v.component(.i.);           
   temp.component(.temp.dim.) := r;
   append := temp;                  
end;  (* function append *)         
(**********************************************************************)
(**********************************************************************)

(*      name:  findeps            *)

(*  category:  mathematical utilities               *)

(*   purpose:  finds an approximation to the machine epsilon          *)

(*     input:  no  input parameters                 *)

(*    output:  eps - approximate machine epsilon    *)


(**********************************************************************)
procedure findeps(var eps:real);    
var               
   epsp1 : real;
begin (* procedure findeps *)       
   eps := 1.0;    
   epsp1 := eps + 1.0;              
   while epsp1 > 1 do               
   begin (* while loop *)           
      eps := 0.5*eps;               
      epsp1 := eps + 1.0;           
   end;  (* while loop *)           
   writeln(report);writeln(report,  
   '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%');          
   writeln(report,' machine epsilon is approximately ',2.0*eps);        
   writeln(report,                  
   '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%');          
   writeln(report);                 
end;  (* procedure findeps *)       
(**********************************************************************)
(**********************************************************************)

(*      name:  symmetric          *)

(*  category:  mathematical utilities               *)

(*   purpose:  given a degree n and a parameter x   *)
(*             values of the symmetric functions of degree n          *)
(*             are computed.  these functions are used by the         *)
(*             deboor_fix alogorithm for computing b-spline           *)
(*             coefficients.      *)

(*     input:  n - the degree desired for the symmetric functions     *)
(*             x - parameter value at which the symmetric functions   *)
(*                 are to be evaluated              *)

(*    output:  s - s(j)(x) is the j-th symmetric function of          *)
(*                 degree n evaluated at x.         *)

(**********************************************************************)
                  
procedure symmetric(n : integer;    
  x : realmd;     
               var  s : realmd);    
var               
   i,j : integer;
stemp : realmdmd;
                  
begin (* procedure symmetric *)     
                  
   for i := 0 to n do stemp(.i,0.) := 1.0;            
   for i := 0 to n do stemp(.i,i+1.) := 0.0;          
                  
   for i:= 1 to n do                
   for j:= 1 to i do                
   stemp(.i,j.) := x(.i.)*stemp(.i-1,j-1.) + stemp(.i-1 ,j.);           
                  
   for j := 0 to n do               
   s(.j.) := stemp(.n,j.);          
                  
end;  (* procedure symmetric *)      
  
(**********************************************************************)
(**********************************************************************)

(*      name:  factorial           *)

(*  category:  mathematical utilities               *)

(*   purpose:  find the factorial of a given integer                  *)

(*     input:  m - an integer      *)

(*    output:  factorial(m) = m]   *)

(**********************************************************************)
  
function factorial(m : integer): integer;             
var                
i,f : integer;     
begin (* function factorial *)       
   if(m=0) or (m=1) then factorial := 1               
   else            
   begin(* iteration*)               
      f:=1;        
      for i:=2 to m do f:=f*i;       
      factorial := f;                
   end; (* iteration *)              
end;  (* function factorial *)       
  
(**********************************************************************)
(**********************************************************************)

(*      name:  power               *)

(*  category:  mathematical utilities               *)

(*   purpose:  since this version of pascal does not have an          *)
(*             exponentiation operator, we have had to write          *)
(*             this function.      *)

(*     input:  i - an integer      *)
(*             j - a non-negative integer           *)

(*    output:  power(i,j) = i**j   *)

(**********************************************************************)
  
function power(i:integer; j:integer):integer;         
var                
   k,p:integer;    
begin (* function power *)           
   if j=0 then power:=1              
   else            
   begin (* iteration *)             
      p := 1;      
      for k:=1 to j do p := p*i;     
      power := p;  
   end;  (* iteration *)             
end;  (* function power *)           
(**********************************************************************)
(**********************************************************************)

(*      name:  rp_fraction         *)

(*  category:  create fractional parts              *)

(*   purpose:  given a rational polynomial, create the numerator      *)
(*             and denominator as separate polynomial entities.       *)
(*             this takes up more space, but is convenient for        *)
(*             evaluation.         *)

(*     input:  rp - a rational polynomial (presumably with the        *)
(*             numerator and denominator uninitialized)               *)

(*    output:  rp - the rational polynomial with the numerator        *)
(*             and denominator calculated.          *)

(**********************************************************************)
  
procedure rp_fraction(var rp : rational_poly);        
var                
i : integer;       
begin (* procedure rp_fraction *)    
  
   rp.numerator.degree := rp.degree; 
   rp.numerator.dim := rp.dim;       
   rp.numerator.coeff := rp.top_coeff;                
  
   rp.denominator.degree := rp.degree;                
   rp.denominator.dim := 1;          
      for i := 0 to rp.degree do     
      begin (* finding denominator.coeff(.i.) *)      
         rp.denominator.coeff(.i.).dim := 1;          
         rp.denominator.coeff(.i.).component(.1.) := rp.bot_coeff(.i.); 
      end;  (* finding denominator.coeff(.i.) *)      
  
end;  (* procedure rp_fraction *)    
  
(**********************************************************************)
(**********************************************************************)

(*      name:  rpsurf_fraction     *)

(*  category:  create fractional parts              *)

(*   purpose:  given a rational polynomial surface, create the        *)
(*             numerator and denominator as separate polynomial       *)
(*             surface entities. this takes up more space, but is     *)
(*             convenient for evaluation.           *)

(*     input:  rp - a rational polynomial surface  (presumably        *)
(*             with the numerator and denominator uninitialized)      *)

(*    output:  rp - the rational polynomial with the numerator        *)
(*             and denominator calculated.          *)

(**********************************************************************)
  
procedure rpsurf_fraction(var rp : rpsurf);           
var             
   i,j : integer;
begin (* procedure rpsurf_fraction *)              
  
   rp.numerator.degree1 := rp.degree1;             
   rp.numerator.degree2 := rp.degree2;             
   rp.numerator.dim := rp.dim;       
   rp.numerator.coeff := rp.top_coeff;             
  
   rp.denominator.degree1 := rp.degree1;           
   rp.denominator.degree2 := rp.degree2;           
   rp.denominator.dim := 1;          
    for i := 0 to rp.degree1 do      
    for j := 0 to rp.degree2 do      
    begin (* finding denominator.coeff(.i,j.) *)   
     rp.denominator.coeff(.i,j.).dim := 1;         
     rp.denominator.coeff(.i,j.).component(.1.) := rp.bot_coeff(.i,j.); 
    end;  (* finding denominator.coeff(.i,j.) *)   
  
end;  (* procedure rpsurf_fraction *)              
          
(**********************************************************************)
(**********************************************************************)
        
(*      name:  rbs_fraction        *)        
        
(*  category:  create fractional parts                                *)
        
(*   purpose:  given a rational b-spline curve, create the            *)
(*             numerator and denominator as separate b-spline         *)
(*             curve entities. this takes up more space, but is       *)
(*             convenient for evaluation.                             *)
        
(*     input:  rbs - a rational b-spline curve  (presumably           *)
(*             with the numerator and denominator uninitialized)      *)
        
(*    output:  rbs - the rational b-spline with the numerator         *)
(*             and denominator calculated.                            *)
        
(**********************************************************************)
          
procedure rbs_fraction(var rbs : rational_bspline);
var       
i : integer;            
begin (* procedure rbs_fraction *)           
          
   rbs.numerator.degree := rbs.degree;       
   rbs.numerator.dim := rbs.dim;             
   rbs.numerator.number_of_knots := rbs.number_of_knots;                
   rbs.numerator.knot := rbs.knot;           
   for i:= 1 to rbs.number_of_knots - rbs.degree - 1 do                 
   rbs.numerator.coeff(.i.) := product(rbs.weight(.i.),rbs.coeff(.i.)); 
          
   rbs.denominator.degree := rbs.degree;     
   rbs.denominator.dim := 1;                 
   rbs.denominator.number_of_knots := rbs.number_of_knots;              
   rbs.denominator.knot := rbs.knot;         
   for i := 1 to rbs.number_of_knots - rbs.degree - 1 do                
   begin (* finding coeff(.i.)*)             
      rbs.denominator.coeff(.i.).dim := 1;   
      rbs.denominator.coeff(.i.).component(.1.) := rbs.weight(.i.);     
   end;  (* finding coeff(.i.) *)            
          
end;  (* procedure rbs_fraction *)           
          
(**********************************************************************)
(**********************************************************************)
        
(*      name:  rbsurf_fraction     *)        
        
(*  category:  create fractional parts                                *)
        
(*   purpose:  given a rational b-spline surface, create the          *)
(*             numerator and denominator as separate b-spline         *)
(*             surface entities. this takes up more space, but is     *)
(*             convenient for evaluation.                             *)
        
(*     input:  rbs - a rational b-spline surface  (presumably         *)
(*             with the numerator and denominator uninitialized)      *)
        
(*    output:  rbs - the rational b-spline surface with the numerator *)
(*             and denominator calculated.                            *)
        
(**********************************************************************)
          
procedure rbsurf_fraction(var b : rbsurf);   
var       
   i,j : integer;       
begin (* procedure rbsurf_fraction *)        
          
   b.numerator.degree1 := b.degree1;         
   b.numerator.degree2 := b.degree2;         
   b.numerator.dim := b.dim;                 
   b.numerator.num_of_knots1 := b.num_of_knots1;   
   b.numerator.num_of_knots2 := b.num_of_knots2;   
   b.numerator.knot1 := b.knot1;             
   b.numerator.knot2 := b.knot2;             
   for i:= 1 to b.num_of_knots1 - b.degree1 - 1 do 
   for j:= 1 to b.num_of_knots2 - b.degree2 - 1 do 
   b.numerator.coeff(.i,j.) := product(b.weight(.i,j.),b.coeff(.i,j.)); 
          
   b.denominator.degree1 := b.degree1;       
   b.denominator.degree2 := b.degree2;       
   b.denominator.dim := 1;
   b.denominator.num_of_knots1 := b.num_of_knots1; 
   b.denominator.num_of_knots2 := b.num_of_knots2; 
   b.denominator.knot1 := b.knot1;           
   b.denominator.knot2 := b.knot2;           
   for i := 1 to b.num_of_knots1 - b.degree1 - 1 do
   for j := 1 to b.num_of_knots2 - b.degree2 - 1 do
   begin (* finding coeff(.i,j.)*)           
      b.denominator.coeff(.i,j.).dim := 1;   
      b.denominator.coeff(.i,j.).component(.1.) := b.weight(.i,j.);     
   end;  (* finding coeff(.i,j.) *)          
          
end;  (* procedure rbsurf_fraction *)        
          
(**********************************************************************)
(**********************************************************************)
        
(*      name:  read_knots          *)        
        
(*  category:  input/output        *)        
        
(*   purpose:  prompt the user and read in a knot sequence.  this     *)
(*             procedure ensures that the input knot sequence is      *)
(*             nondecreasing.      *)        
        
(*     input:  degree - the degree of the b-spline curve              *)
        
(*    output:  number_of_knots - the number of knots, including       *)
(*                               multiplicity                         *)
(*             knot - the knot sequence.  start and end values        *)
(*                    are given multiplicity = degree + 1             *)
        
(**********************************************************************)
   procedure read_knots(degree : integer;    
                   var  number_of_knots : integer; 
                   var  knot : real44);      
var       
   i : integer;         
   temp : real44;       
   nondecreasing : boolean;                  
          
begin (* procedure read_knots *)             
   writeln(' enter number of knots (counting multiplicities) ');        
   repeat (* until number_of_knots >= 2*(degree + 1) *)                 
   writeln(' the number of knots must be at least ',
   2*(degree + 1):1);   
   readln(number_of_knots);     
   until number_of_knots >= 2*(degree + 1);        
            
   repeat (* until nondecreasing and knot(1) < knot(number_of_knots) *) 
   writeln(' enter the knot values.');             
   writeln(' the first and last values you enter will');                
   writeln(' automatically be given multiplicity ',degree+1:1);         
   writeln(' therefore, enter ',number_of_knots-2*degree:1,' values');  
            
   for i := 1 to number_of_knots - 2*degree do read(temp(.i.));         
            
   for i := 1 to degree + 1 do 
   begin (* creating multiple end knots*)          
   knot(.i.) := temp(.1.);     
   knot(.number_of_knots - i + 1.):=temp(.number_of_knots-2*degree.);   
   end;  (* creating multiple end knots *)         
            
   if number_of_knots > 2*(degree + 1) then        
   for i := 1 to number_of_knots - 2*(degree + 1) do
   knot(.i + degree + 1.):= temp(.i + 1.);         
            
   nondecreasing := true;      
   for i := 1 to number_of_knots - 1 do            
   nondecreasing := nondecreasing and              
   (knot(.i.) <= knot(.i + 1.));
            
   if nondecreasing = false then
   writeln(' knot sequence must be nondecreasing. ');
            
   if knot(.1.) >= knot(.number_of_knots.) then    
   writeln(' last knot value must be greater than first knot value.');  
            
   if (nondecreasing = false) or
   (knot(.1.) >= knot(.number_of_knots.)) then     
   begin    
      writeln(' try again] '); 
      writeln;                 
   end;     
            
   until nondecreasing and (knot(.1.)<knot(.number_of_knots.));         
            
end;  (* procedure read_knots *)
(**********************************************************************)
%page       
(**********************************************************************)

(*      name:  read_rational_bspline                                  *)

(*  category:  input/output                                           *)

(*   purpose:  prompt the user and read in a rational b-spline        *)
(*             curve.                                                 *)

(*     input:  no input parameters.  input is from the user.          *)

(*    output:  g - a rational b-spline curve                          *)

(**********************************************************************)
procedure read_rational_bspline(var g : rational_bspline);              
var         
   temp : real44;              
   nondecreasing : boolean;    
   i,j : integer;              
   order : integer;            
begin (* procedure read_rational_bspline *)        
   writeln; writeln;           
   writeln(' enter degree desired for rational bspline :');             
   readln(g.degree);           
   if g.degree > 0 then        
   with g do
   begin (* reading in data for rational bspline *)
            
      read_knots(degree,g.number_of_knots,g.knot); 
            
      order := g.degree + 1;   
      writeln(' enter dimension for b-spline coefficients :');          
      readln(g.dim);           
            
      writeln(' enter ',g.number_of_knots - order:1,' coefficients.');  
      for i := 1 to g.number_of_knots - order do   
      begin 
         for j:= 1 to dim do read(coeff(.i.).component(.j.));           
         coeff(.i.).dim := g.dim;                  
      end;  
            
            
      writeln(' enter ',g.number_of_knots - order:1,' weights.');       
      for i := 1 to g.number_of_knots - order do   
      read(weight(.i.));       
            
   end;  (* reading in data for rational bspline *)
            
rbs_fraction(g);               
            
end;  (* procedure read_rational_bspline *)        
            
(**********************************************************************)
(**********************************************************************)

(*      name:  read_rbsurf                                            *)

(*  category:  input/output                                           *)

(*   purpose:  prompt the user and read in a rational b-spline        *)
(*             surface                                                *)

(*     input:  no input parameters.  input is from the user.          *)

(*    output:  g - a rational b-spline surface                        *)

(**********************************************************************)
procedure read_rbsurf(var b : rbsurf);             
var         
   temp : real44;              
   nondecreasing : boolean;    
   i,j,k : integer;            
   order1,order2 : integer;    
begin (* procedure read_rbsurf *)                  
   writeln; writeln;           
   writeln(' enter degree in 1st parameter :');    
   readln(b.degree1);          
   writeln(' enter degree in 2nd parameter :');    
   readln(b.degree2);          
   if (b.degree1 > 0) and (b.degree2 > 0) then     
   with b do
   begin (* reading in data for rational bspline surface *)             
            
   writeln(' define first knot sequence:  ');      
   read_knots(degree1,num_of_knots1,knot1);        
            
   writeln(' define second knot sequence:  ');     
   read_knots(degree2,num_of_knots2,knot2);        
            
   order1 := degree1 + 1;      
   order2 := degree2 + 1;      
            
   writeln(' enter dimension for b-spline coefficients :');             
   readln(dim);                
            
   i := num_of_knots1 - order1;
   j := num_of_knots2 - order2;
   writeln(' enter the ',i:1,'x',j:1,' coefficients as follows:');      
   for i := 1 to num_of_knots1 - order1 do         
   for j := 1 to num_of_knots2 - order2 do         
   begin(* reading in the i,j-th coefficient *)    
      writeln('enter coefficient(',i:1,',',j:1,'):  ');                 
      for k:= 1 to dim do read(coeff(.i,j.).component(.k.));            
      coeff(.i,j.).dim := dim; 
   end;(* reading in the i,j-th coefficient *)     
            
   i := num_of_knots1 - order1; 
   j := num_of_knots2 - order2; 
             
   repeat (* until k = 0 or k = 1 *)                
   writeln(' enter 0 for rational case, 1 otherwise');
   writeln(' if you enter 1, all weights will be set to 1');             
   readln(k);
   if k = 0 then                
   begin (* reading in weights *)
      writeln(' enter the ',i:1,'x',j:1,' weights as follows:');         
      for i := 1 to num_of_knots1 - order1 do       
      for j := 1 to num_of_knots2 - order2 do       
      begin(* reading in the i,j-th weight *)       
         writeln('enter weight(',i:1,',',j:1,'):  ');
         readln(weight(.i,j.)); 
      end;(* reading in the i,j-th weight *)        
   end   (* reading in weights *)
   else      
   if k = 1 then                
   for i := 1 to num_of_knots1 - order1 do          
   for j := 1 to num_of_knots2 - order2 do          
   weight(.i,j.) := 1.0;        
   until (k = 0) or (k = 1);    
             
  rbsurf_fraction(b);           
             
             
   end;  (* reading in data for rational bspline surface *)              
end;  (* procedure read_rbsurf *)
             
(**********************************************************************) 
(**********************************************************************) 
 
(*      name:  write_rbsurf                                           *) 
 
(*  category:  input/output                                           *) 
 
(*   purpose:  write the defining data for a rational b-spline        *) 
(*             surface to a designated output file                    *) 
 
(*     input:  outfil - the file to which the data is to be written   *) 
(*             b - the rational b-spline surface                      *) 
 
(*    output:  no output parameters.  the defining data for b         *) 
(*             is written to outfil.                                  *) 
 
(**********************************************************************) 
procedure write_rbsurf(var outfil : text; b : rbsurf);
var          
   i,j,k : integer;             
   order1,order2 : integer;     
             
begin (* procedure write_rbsurf *)                  
   with b do 
   begin (* writing values of b to outfil *)        
   writeln(outfil,' ');         
   writeln(outfil,' degree in 1st parameter =',degree1:1);               
   writeln(outfil,' degree in 2nd parameter =',degree2:1);               
   writeln(outfil,' ');         
   writeln(outfil,' no of knots in 1st parameter = ',num_of_knots1:1);   
   for i := 1 to num_of_knots1 do
   writeln(outfil,' knot # ',i:2,' = ',knot1(.i.):20:15);                
   writeln(outfil,' ');         
   writeln(outfil,' no of knots in 2nd parameter = ',num_of_knots2:1);   
   for i := 1 to num_of_knots2 do
   writeln(outfil,' knot # ',i:2,' = ',knot2(.i.):20:15);                
             
   order1 := degree1 + 1;       
   order2 := degree2 + 1;       
   writeln(outfil);             
   for i := 1 to num_of_knots1 - order1 do          
   for j := 1 to num_of_knots2 - order2 do          
   begin (* writing (i-j)th coeff and wgt *)        
      write(outfil,'coeff(',i:1,',',j:1,') = ');    
      writevec(outfil,coeff(.i,j.));                
      write(outfil,'    wgt(',i:1,',',j:1,') = ');  
      write(outfil,weight(.i,j.):10:6);             
      writeln(outfil);          
   end;  (* writing (i-j)th coeff and wgt *)        
   end;  (* writing values of b to outfil *)        
end;  (* procedure write_rbsurf *)                  
             
(**********************************************************************) 
(**********************************************************************) 
 
(*      name:  write_rbs                                              *) 
 
(*  category:  input/output                                           *) 
 
(*   purpose:  write the defining data for a rational b-spline        *) 
(*             curve to a designated output file                      *) 
 
(*     input:  outfil - the file to which the data is to be written   *) 
(*             b - the rational b-spline curve                        *) 
 
(*    output:  no output parameters.  the defining data for g         *) 
(*             is written to outfil.                                  *) 
 
(**********************************************************************) 
procedure write_rbs(var outfil : text; g : rational_bspline);            
var          
   i : integer;                 
   order : integer;             
             
begin (* procedure write_rbs *) 
   with g do 
   begin (* writing values of g to outfil*)         
   order := degree + 1;         
   writeln(outfil,' ');         
   writeln(outfil,' degree = ',degree:1);           
   writeln(outfil,' number of knots =',number_of_knots:1);               
   writeln(outfil,' number of coeffs = ',number_of_knots-order:1);       
   writeln(outfil,' ');         
             
   for i:= 1 to number_of_knots do                  
   writeln(outfil,' knot # ',i:1,' = ',knot(.i.):20:15);                 
             
   writeln(outfil,' ');         
   writeln(outfil,' ');         
   for i:=1 to number_of_knots - order do           
   begin     
      write(outfil,' coeff(',i:2,')=');             
      writevec(outfil,coeff(.i.));                  
      writeln(outfil);          
   end;      
   writeln(outfil,' ');         
   for i:=1 to number_of_knots - order do           
   begin     
      write(outfil,' weight(',i:2,')=');            
      write(outfil,weight(.i.):15:10);              
      writeln(outfil);          
   end;      
   end; (* writing values of g to outfil *)         
end;  (* procedure write_rbs *) 
             
(**********************************************************************) 
(**********************************************************************) 
 
(*      name:  read_rational_pp                                       *) 
 
(*  category:  input/output                                           *) 
 
(*   purpose:  prompt the user and read in a rational piecewise       *) 
(*             polynomial                                             *) 
 
(*     input:  no input parameters.  input is from the user.          *) 
 
(*    output:  g - a rational piecewise polynomial                    *) 
 
(**********************************************************************) 
procedure read_rational_pp(var f : rational_pp);    
var          
   d,i,j,k : integer;           
   increasing : boolean;        
   order : integer;             
begin (* procedure read_rational_pp *)              
   writeln; writeln;            
   writeln(' enter degree desired for piecewise polynomials ');          
   readln(f.degree);   d := f.degree;               
   if f.degree >= 1 then        
   begin (* normal case *)      
   with f do 
   begin (* reading in data for rational_pp *)      
             
   order := degree + 1;         
   writeln(' enter number of segments:   ');        
   readln(num_segs);            
             
   repeat (* until increasing *)
   writeln(' enter ',num_segs+1:1,' breakpoints:  ');
   for i := 1 to num_segs+1 do read(break_point(.i.));
   increasing := true;          
             
   for i:= 1 to num_segs do     
   increasing := increasing and 
   (break_point(.i.) < break_point(.i + 1.));       
             
   if increasing=false then     
   writeln(' breakpoints were not strictly increasing.  try again');     
   until increasing;            
             
   writeln(' enter dimension for coefficients: ');  
   readln(dim);                 
             
   cms('clrscrn',return_code);  
   writeln(' each segment is expressed in the form:');
   writeln('        top_coeff(.0.)+...+top_coeff(.',d:1,'.)*s**',d:1);   
   writeln(' f(s) = ----------------------------------------------');    
   writeln('        bot_coeff(.0.)+...+bot_coeff(.',d:1,'.)*s**',d:1);   
   writeln;  
   writeln(' where s is the distance from the left end point');          
   writeln(' ');                
for i := 1 to num_segs do       
   begin(* reading in data for segment i *)         
      rpp_segment(.i.).degree := degree;            
writeln(' enter top_coeff(.0.),...,top_coeff(.',d:1,'.) for seg #',i:2); 
      for j := 0 to degree do   
      begin  
     for k := 1 to dim do       
      read(rpp_segment(.i.).top_coeff(.j.).component(.k.));              
      rpp_segment(.i.).top_coeff(.j.).dim := dim;   
      end;   
writeln(' enter bot_coeff(.0.),...bot_coeff(.',d:1,'.) for seg #',i:2);  
      for j := 0 to degree do read(rpp_segment(.i.).bot_coeff(.j.));     
             
   rp_fraction(rpp_segment(.i.));
             
   end;(* reading in data for segment i *)          
end; (* reading in data for rational_pp *)          
end; (* normal case *)          
end;  (* procedure read_rational_pp *)              
             
(**********************************************************************) 
%page        
(**********************************************************************)

(*      name:  read_rppsurf                                           *)

(*  category:  input/output                                           *)

(*   purpose:  prompt the user and read in a rational piecewise       *)
(*             polynomial surface                                     *)

(*     input:  no input parameters.  input is from the user.          *)

(*    output:  g - a rational piecewise polynomial surface            *)

(**********************************************************************)
procedure read_rppsurf(var f : rppsurf);           
var         
   d1,d2,order1,order2,i,j,k,m,n : integer;        
   u,v : real;                 
   del1,del2 : real;           
   start1,start2,end1,end2 : real;                 
   rppvalue : vector;          
   increasing : boolean;       
begin (* procedure read_rppsurf *)                 
   writeln; writeln;           
   writeln(' enter degree desired in first parameter ');                
   readln(f.degree1);   d1 := f.degree1;           
   writeln(' enter degree desired in second parameter ');               
   readln(f.degree2);   d2 := f.degree2;           
   if(f.degree1 >= 1) and (f.degree2 >= 1) then    
   begin (* normal case *)     
   with f do
   begin (* reading in data for rppsurf *)         
            
   order1 := degree1 + 1;      
   order2 := degree2 + 1;      
            
   writeln(' enter number of segments in first direction:   ');         
   readln(num_segs1);          
   writeln(' enter number of segments in second direction:   ');        
   readln(num_segs2);          
            
   repeat (* until increasing *)
   writeln(' enter ',num_segs1+1:1,' breakpoints in 1st dir:  ');       
   for i := 1 to num_segs1+1 do read(break_point1(.i.));                
   increasing := true;         
            
   for i:= 1 to num_segs1 do   
   increasing := increasing and
   (break_point1(.i.) < break_point1(.i + 1.));    
            
   if increasing=false then    
   writeln(' breakpoints were not strictly increasing.  try again');    
   until increasing;           
            
   repeat (* until increasing *)
   writeln(' enter ',num_segs2+1:1,' breakpoints in 2nd dir:  ');       
   for i := 1 to num_segs2+1 do read(break_point2(.i.));                
   increasing := true;         
            
   for i:= 1 to num_segs2 do   
   increasing := increasing and
   (break_point2(.i.) < break_point2(.i + 1.));    
            
   if increasing=false then    
   writeln(' breakpoints were not strictly increasing.  try again');    
   until increasing;           
            
   writeln(' enter dimension for coefficients: '); 
   readln(dim);                
            
   cms('clrscrn',return_code); 
   writeln(' each patch is expressed in the form:');
   writeln('           sum{top_coeff(m,n)*(s**m)*(t**n)');              
   writeln(' f(s,t) = ----------------------------------------------'); 
   writeln('           sum{bot_coeff(m,n)*(s**m)*(t**n)');              
   writeln;                
   writeln(' where s,t are the distances from the left end points');    
   writeln(' ');           
for i := 1 to num_segs1 do 
for j := 1 to num_segs2 do 
   begin(* reading in data for patch (i,j) *)   
      writeln(' begin to enter data for patch (',i:1,',',j:1,'):');     
      rppatch(.i,j.).degree1 := degree1;        
      rppatch(.i,j.).degree2 := degree2;        
      for m := 0 to degree1 do                  
      for n := 0 to degree2 do     
      begin (* reading in (m,n)-th top coeff *)         
         writeln(' enter top_coeff(',m:1,',',n:1,')');  
         for k := 1 to dim do      
         read(rppatch(.i,j.).top_coeff(.m,n.).component(.k.));                  
         rppatch(.i,j.).top_coeff(.m,n.).dim :=dim;     
      end;  (* reading in (m,n)-th top coeff *)         
               
      for m := 0 to degree1 do     
      for n := 0 to degree2 do     
      begin (* reading in (m,n)-th bot coeff *)         
         writeln(' enter bot_coeff(',m:1,',',n:1,')');  
         read(rppatch(.i,j.).bot_coeff(.m,n.));         
      end;  (* reading in (m,n)-th bot coeff *)         
      rpsurf_fraction(rppatch(.i,j.));                  
   cms('clrscrn',return_code);     
   end;(* reading in data for patch (i,j) *)            
end; (* reading in data for rppsurf *)                  
end; (* normal case *)             
end;  (* procedure read_rppsurf *) 
               
(**********************************************************************)        
%page          
(**********************************************************************)        
        
(*      name:  write_rppsurf                                          *)        
        
(*  category:  input/output                                           *)        
        
(*   purpose:  write the defining data for a rational piecewise       *)        
(*             polynomial surface to a designated output file         *)        
        
(*     input:  outfil - the file to which the data is to be written   *)        
(*             f - the rational piecewise polynomial surface          *)        
        
(*    output:  no output parameters.  the defining data for f         *)        
(*             is written to outfil.                                  *)        
        
(**********************************************************************)        
procedure write_rppsurf(var outfil : text; f : rppsurf);
var            
   i,j,m,n : integer;              
               
begin (* procedure write_rppsurf *)
with f do      
begin (* writing values of f to outfil *)               
   writeln(outfil,' ');            
   writeln(outfil,' degree1 = ',degree1:1);             
   writeln(outfil,' degree2 = ',degree2:1);             
   writeln(outfil,' num_segs1 = ',num_segs1:1);         
   writeln(outfil,' num_segs2 = ',num_segs2:1);         
   writeln(outfil,' ');            
   for i := 1 to num_segs1 + 1 do writeln(outfil,       
   ' break_point1(',i:1,')=',break_point1(.i.):20:15);  
   writeln(outfil,' ');            
   for i := 1 to num_segs2 + 1 do writeln(outfil,       
   ' break_point2(',i:1,')=',break_point2(.i.):20:15);  
   writeln(outfil,' ');            
   for i:= 1 to num_segs1 do       
   for j:= 1 to num_segs2 do       
   begin (* writing top coeffs for (i,j)-th patch *)    
      writeln(outfil,' ');         
      for m:= 0 to degree1 do      
      for n:= 0 to degree2 do      
      begin (* writing (m,n)-th top coefficient *)      
         write(outfil,             
         ' patch # (',i:1,',',j:1,')   top_coeff(',m:1,',',n:1,') = ');         
         writevec(outfil,f.rppatch(.i,j.).top_coeff(.m,n.));
         writeln(outfil);          
      end;  (* writing (m,n)-th top coefficient *)      
   end;  (* writing top coefficients for (i,j)-th patch *) 
   writeln(outfil,' ');            
   for i:= 1 to num_segs1 do       
   for j:= 1 to num_segs2 do       
   begin (* writing bottom coefficients for (i,j)-th patch *)
      writeln(outfil,' ');         
      for m:= 0 to degree1 do      
      for n:= 0 to degree2 do      
      begin(* writing (m,n)-th bottom coefficient *)    
      write(outfil,                
      ' patch # (',i:1,',',j:1,')   bot_coeff(',m:1,',',n:1,') = ');            
      write(outfil,f.rppatch(.i,j.).bot_coeff(.m,n.):20:15);
      writeln(outfil);             
      end; (* writing (m,n)-th bottom coefficient *)    
   end;  (* writing bottom coefficients for (i,j)-th patch *)
end;  (* writing values of f to outfil *)               
               
end;  (* procedure write_rppsurf *)
(**********************************************************************)        
%page          
(**********************************************************************)        
        
(*      name:  write_rpp                                              *)        
        
(*  category:  input/output                                           *)        
        
(*   purpose:  write the defining data for a rational piecewise       *)        
(*             polynomial curve to a designated output file           *)        
        
(*     input:  outfil - the file to which the data is to be written   *)        
(*             f - the rational piecewise polynomial curve            *)        
        
(*    output:  no output parameters.  the defining data for f         *)        
(*             is written to outfil.                                  *)        
        
(**********************************************************************)        
procedure write_rpp(var outfil : text; f : rational_pp);
var            
   i,j : integer;                  
               
begin (* procedure write_rpp *)    
   with f do   
   begin (* writing values of f to outfil *)            
   writeln(outfil,' ');            
   writeln(outfil,' degree = ',degree:2);               
   writeln(outfil,' num_segs = ',num_segs:2);           
   writeln(outfil,' ');            
   for i := 1 to num_segs + 1 do writeln(outfil,        
   ' break_point(.',i:1,'.)=',break_point(.i.):20:15);  
   writeln(outfil,' ');            
   for i:= 1 to num_segs do        
   begin (* writing top coeffs for i-th segment *)      
      writeln(outfil,' ');         
      for j:= 0 to degree do       
      begin (* writing j-th top coefficient *)          
         write(outfil,' seg # = ',i:1,'   top_coeff(.',j:1,'.) = ');            
         writevec(outfil,f.rpp_segment(.i.).top_coeff(.j.));
         writeln(outfil);          
      end;  (* writing j-th top coefficient *)          
   end;  (* writing top coefficients for i-th segment *)
   writeln(outfil,' ');            
   for i:= 1 to num_segs do        
   begin (* writing bottom coefficients for i-th segment *)
      writeln(outfil,' ');         
      for j:= 0 to degree do       
      begin(* writing j-th bottom coefficient *)        
         write(outfil,' seg # = ',i:1,'   bot_coeff(.',j:1,'.) = ');            
         write(outfil,f.rpp_segment(.i.).bot_coeff(.j.):20:15);                 
         writeln(outfil);          
      end;  (* writing j-th bottom coefficient *)       
   end;  (* writing bottom coefficients for i-th segment *)
   end;  (* writing values of f to outfil *)            
               
end;  (* procedure write_rpp *)    
(**********************************************************************)        
   %page       
(**********************************************************************)        
        
(*      name:  bracket                                                *)        
        
(*  category:  interval finding                                       *)        
        
(*   purpose:  given a parameter and a set of indexed break_points,   *)        
(*             finds an index so that the corresponding break_points  *)        
(*             bracket the parameter                                  *)        
        
(*     input:  num_segs - the number of segments determined by        *)        
(*                        the break points                            *)        
(*             break_point - the array of break_points                *)        
(*             u - the parameter to be bracketed                      *)        
        
(*    output:  bracket(num_segs,break_point,u) = i, where i           *)        
(*             is the largest integer such that                       *)        
(*             break_point(.i.) <= u <= break_point(.i+1.)            *)        
        
(**********************************************************************)        
function bracket(num_segs : integer;
                 break_point : real22;                  
                 u : real) : integer;
var            
   i:integer;  
begin (* function bracket *)       
        for i:= num_segs downto 1 do
        if(break_point(.i.) <= u) and (u <= break_point(.i + 1.))               
        then   
        begin  
           bracket := i;           
           leave;                  
        end;   
 end; (* function bracket *)       
               
(**********************************************************************)        
%page          
(**********************************************************************)        
        
(*      name:  knot_interval                                          *)        
        
(*  category:  interval finding                                       *)        
        
(*   purpose:  given a parameter and a set of indexed knots           *)        
(*             finds an index so that the corresponding knots         *)        
(*             bracket the parameter                                  *)        
        
(*     input:  number_of_knots - the number of knots                  *)        
(*             t -  the array of knots                                *)        
(*             u - the parameter to be bracketed                      *)        
        
(*    output:  bracket(number_of_knots,t,u)    = i, where i           *)        
(*             is the largest integer such that                       *)        
(*             t(.i.) <= u <= t(.i+1.)                                *)        
        
(**********************************************************************)        
               
function knot_interval(number_of_knots:integer;         
                                     t:real44;          
                                     u:real):integer;   
var            
i:integer;     
begin (* function knot_interval *) 
        for i:= number_of_knots - 1 downto 1 do         
        if(t(.i.) < t(.i+1.)) and  
          (t(.i.) <= u) and (u <= t(.i + 1.)) then      
          begin
             knot_interval:=i;     
             leave;                
          end; 
 end; (* function knot_interval *) 
(**********************************************************************)        
%page          
(**********************************************************************)        
        
(*      name:  leeval                                                 *)        
        
(*  category:  evaluation                                             *)        
        
(*   purpose:  evaluates the value and derivatives of a linear        *)        
(*             combination of b-spline basis functions over one span  *)        
(*             the method used is that described in the paper         *)        
(*             " a simplified b-spline computation routine "          *)        
(*             by e. t. y. lee which appeared in computing, vol 29,   *)        
(*             pp. 365-371 (1982).                                    *)        
        
(*     input:  m - degree of b-spline basis functions                 *)        
(*             t - t(.0.),...,t(.2*m+1.) is the knot sequence         *)        
(*             c - c(.0.),...,c(.m.) are the b-spline coefficients    *)        
(*             x - parameter to evaluate at: t(.m.) <= x <= t(.m+1.)  *)        
(*             d - highest derivative desired:  0 <= d <= m           *)        
        
(*    output:  q - q(.j.) = j-th deriv of function at x               *)        
        
(**********************************************************************)        
procedure leeval(m:integer;        
                 t:real44;         
                 c:vec44;          
                 x:real;           
                 d:integer;        
             var q:vecmd);         
var            
   a,b : real44;
   i,j:integer;
   fromright:boolean;              
   vec1,vec2,vec3 : vector;        
               
begin (* procedure leeval *)       
/*             
if idebug = 1 then                 
begin (* writing to debug file *)  
writeln(debug);
writeln(debug,' entering leeval with: ');               
writeln(debug,' m = ',m:3);        
writeln(debug,' t= ');for i:=0 to 2*m+1                 
do write(debug,t(.i.):3:0);        
writeln(debug,' ');                
writeln(debug,' c= ');             
for i := 0 to m do                 
begin          
writevec(debug,c(.i.));            
writeln(debug,' ');                
end;           
writeln(debug,' d= ',d);           
writeln(debug,' x= ',x:4:2);       
end; (* writing to debug file *)   
*/             
   if m=0 then q(.0.):=c(.0.)      
else           
begin (* non-trivial case *)       
   for i:=1 to m do                
   begin       
      a(.i.) := x - t(.i.);        
      b(.i.) := t(.m+i.) - x;      
   end;        
               
fromright := (b(.1.) > a(.m.));    
               
if fromright then                  
begin (* fromright *)              
   for i:= 0 to m do               
   q(.i.) := c(.i.);               
               
(*------------------------------------------------------------------*)          
   for j:=1 to m do                
   for i:=0 to m-j do              
   begin       
   (* =============================================================             
   q(.i.):=(a(.i+j.)*q(.i+1.)+b(.i+1.)*q(.i.))/(a(.i+j.)+b(.i+1.));             
   ==============================================================*)             
   vec1 := product(a(.i+j.),q(.i+1.));                  
   vec2 := product(b(.i+1.),q(.i.));
   vec3:= sum(vec1,vec2);          
   q(.i.) := quotient(vec3,a(.i+j.)+b(.i+1.));          
   end;        
(*------------------------------------------------------------------*)          
               
(*------------------------------------------------------------------*)          
   for j:= 1 to d do               
   for i:= d downto j do           
   begin       
   (* =============================================================             
   q(.i.):=(q(.i.)-q(.i-1.))/(b(.i-j+1.)/(m-j+1));      
   ==============================================================*)             
   vec1 := difference(q(.i.),q(.i-1.));                 
   q(.i.) := quotient(vec1,b(.i-j+1.)/(m-j+1))          
   end;        
(*------------------------------------------------------------------*)          
               
   end (* fromright *)             
   else        
   begin (* not fromright *)       
   for i:= 0 to m do               
   q(.i.) := c(.m-i.);             
               
(*------------------------------------------------------------------*)          
   for j:= 1 to m do               
   for i:= 0 to m-j do             
   begin       
   (* =============================================================             
q(.i.):=(b(.m-j+1-i.)*q(.i+1.)+a(.m-i.)*q(.i.))/(a(.m-i.)+b(.m-j+1-i.));        
   ==============================================================*)             
   vec1 := product(b(.m-j+1-i.),q(.i+1.));              
   vec2 := product(a(.m-i.),q(.i.));
   vec3 := sum(vec1,vec2);         
   q(.i.) := quotient(vec3,a(.m-i.)+b(.m-j+1-i.));      
   end;        
(*------------------------------------------------------------------*)          
               
(*------------------------------------------------------------------*)          
   for j:=1 to d do                
   for i := d downto j do          
   begin       
   (* =============================================================             
   q(.i.):=(-1.0)*(q(.i.)-q(.i-1.))/(a(.m-i+j.)/(m-j+1));  
   ==============================================================*)             
   vec1 := difference(q(.i-1.),q(.i.));                 
   q(.i.) := quotient(vec1,a(.m-i+j.)/(m-j+1));         
   end;        
(*------------------------------------------------------------------*)          
               
   end; (* not fromright *)        
   end; (* non-trivial case *)     
/*             
   if idebug = 1 then              
   begin (* writing to debug *)    
   writeln(debug);                 
   for i := 0 to d do              
   begin       
   write(debug,' leaving leeval with q(',i:1,')= ');    
   writevec(debug,q(.i.));         
   writeln(debug);                 
   end;        
   writeln(debug);                 
   end;  (* writing to debug *)    
   */          
   end; (* procedure leeval *)     
               
(**********************************************************************)        
%page          
(**********************************************************************)        
        
(*      name:  horner                                                 *)        
        
(*  category:  evaluation                                             *)        
        
(*   purpose:  uses horner's method to evaluate a polynomial          *)        
(*             with vector coefficients                               *)        
        
(*     input:  p - a polynomial                                       *)        
(*             x - a parameter value                                  *)        
        
(*    output:  horner(p,x) = p(x) evaluated by horner's method        *)        
        
(**********************************************************************)        
               
function horner(p:poly; x:real): vector;                
var            
   h:vector;   
   i:integer;  
begin (* function horner *)        
   with p do   
   begin (* horner's method *)     
      h := coeff(.degree.);        
      for i:=degree-1 downto 0 do  
      (*=======================================================                 
      h := h*x + coeff(.i.);       
      =======================================================*)                 
      h := sum(product(x,h),coeff(.i.));                
               
      horner:=h;
   end;  (* horner's method *)     
end;   (* function horner *)       
               
(**********************************************************************)        
%page          
(**********************************************************************)        
        
(*      name:  formal_derivative                                      *)        
        
(*  category:  evaluation                                             *)        
        
(*   purpose:  finds the polynomial which is the derivative of        *)        
(*             a given polynomial.  the derivative is not             *)        
(*             evaluated at any parameter value in this               *)        
(*             procedure                                              *)        
        
(*     input:  p - a polynomial with vector coefficients              *)        
        
(*    output:  q - the polynomial which is the derivative of p        *)        
        
(**********************************************************************)        
               
procedure formal_derivative(p:poly; (* the input polynomial *)                  
                        var q:poly  (* the derivative of p *)  );               
var            
   i:integer;  
begin(* procedure formal_derivative *)                  
   if p.degree = 0 then            
   begin (* zero degree case *)    
      q.degree := 0;               
      q.dim := p.dim;              
      q.coeff(.0.) := zerovec(p.dim);
   end   (* zero degree case *)    
   else        
   begin (* positive degree case *)
      q.degree := p.degree - 1;    
      for i := 0 to p.degree - 1 do
      (*====================================================
      q.coeff(.i.) := (i+1)*p.coeff(.i+1.);             
      ====================================================*)
      q.coeff(.i.) := product(i+1,p.coeff(.i+1.));      
               
   end   (* positive degree case *)
end; (* procedure formal derivative *)                  
               
(**********************************************************************)        
%page          
(**********************************************************************)        
        
(*      name:  eval_poly                                              *)        
        
(*  category:  evaluation                                             *)        
        
(*   purpose:  evaluates a vector-valued polynomial and its           *)        
(*             derivatives at a given parameter value                 *)        
        
(*     input:  p - a vector valued polynomial                         *)        
(*             x - a parameter value                                  *)        
(*             d - maximum derivative desired                         *)        
        
(*    output:  pvalue(i) = i-th derivative of p at x (i = 0,...d )    *)        
        
(**********************************************************************)        
               
procedure eval_poly(p : poly;      
                    x : real;      
                    d : integer;   
          var pvalue  : vecmd);    
var            
   i : integer;
   q : array(.0..20.) of poly;     
begin (* procedure eval_poly *)    
   q(.0.) := p;
   pvalue(.0.) := horner(q(.0.),x);
   for i := 1 to d do              
   begin (* evaluating i-th derivative *)               
      formal_derivative(q(.i-1.),q(.i.));               
      pvalue(.i.) := horner(q(.i.),x);                  
   end;  (* evaluating i-th derivative *)               
end;  (* procedure eval_poly *)    
(**********************************************************************)        
%page          
(**********************************************************************)        
        
(*      name:  eval_psurf                                             *)        
        
(*  category:  evaluation                                             *)        
        
(*   purpose:  evaluates a vector-valued polynomial of two            *)        
(*             variables at a given pair of parameter values          *)        
        
(*     input:  p - the vector-valued polynomial of two variables      *)        
(*             d1 - maximum derivative desired in first parameter     *)        
(*             d2 - maximum derivative desired in second parameter    *)        
(*             s - first parameter                                    *)        
(*             t - second parameter                                   *)        
        
(*    output:  svalue(i,j) = (i,j)-th mixed partial derivative        *)        
(*             of p at (s,t)                                          *)        
        
(**********************************************************************)        
               
procedure eval_psurf(p : psurf;    
                    d1 : integer;  
                    d2 : integer;  
                     s : real;     
                     t : real;     
            var svalue :vecmdmd);  
var            
   i,j,k : integer;                
   f,g : poly; 
   pvalue : vecmd;                 
   q : vecmdmd;
begin (* procedure eval_psurf *)   
   g.degree := p.degree2;          
   g.dim := p.dim;                 
      f.degree := p.degree1;       
      f.dim := p.dim;              
               
   for j := 0 to p.degree2 do      
   begin (* j loop *)              
               
      for i := 0 to p.degree1 do   
      f.coeff(.i.) := p.coeff(.i,j.);
               
      eval_poly(f,s,d1,pvalue);    
      for k := 0 to d1 do          
      q(.j,k.) := pvalue(.k.);     
               
   end;  (* j loop *)              
               
   for i := 0 to d1 do             
   begin (* i loop *)              
               
      for k := 0 to p.degree2 do   
      g.coeff(.k.) := q(.k,i.);    
               
      eval_poly(g,t,d2,pvalue);    
               
      for j := 0 to d2 do          
      svalue(.i,j.) := pvalue(.j.);
               
   end;  (*i loop *)               
               
end;  (* procedure eval_psurf *)   
(**********************************************************************)        
%page          
(**********************************************************************)        
        
(*      name:  eval_ppsurf                                            *)        
        
(*  category:  evaluation                                             *)        
        
(*   purpose:  evaluates a vector-valued piecewise polynomial of two  *)        
(*             variables at a given pair of parameter values          *)        
        
(*     input:  p - the vector-valued piecewise polynomial             *)        
(*                 of two variables                                   *)        
(*             d1 - maximum derivative desired in first parameter     *)        
(*             d2 - maximum derivative desired in second parameter    *)        
(*             u - first parameter                                    *)        
(*             v - second parameter                                   *)        
        
(*    output:  svalue(i,j) = (i,j)-th mixed partial derivative        *)        
(*             of p at (u,v)                                          *)        
        
(**********************************************************************)        
               
procedure eval_ppsurf(p : ppsurf;  
                     d1 : integer; 
                     d2 : integer; 
                      u : real;    
                      v : real;    
             var svalue :vecmdmd); 
var            
   i,j : integer;                  
   s,t : real; 
begin (* procedure eval_ppsurf *)  
   with p do   
   begin (* with p *)              
      i :=  bracket(num_segs1,break_point1,u);          
      j :=  bracket(num_segs2,break_point2,v);          
               
      s := u - break_point1(.i.);  
      t := v - break_point2(.j.);  
               
      eval_psurf(ppatch(.i,j.),d1,d2,s,t,svalue);       
               
   end;  (* with p *)              
               
end;  (* procedure eval_ppsurf *)  
(**********************************************************************)        
%page          
(**********************************************************************)        
        
(*      name:  eval_rppsurf                                           *)        
        
(*  category:  evaluation                                             *)        
        
(*   purpose:  evaluates a rational piecewise polynomial              *)        
(*             at a given pair of parameter values.                   *)        
(*             derivatives are not evaluated; only the function.      *)        
        
(*     input:  p : a rational piecewise vector-valued polynomial      *)        
(*             u - the first parameter                                *)        
(*             v - the second parameter                               *)        
        
(*    output:  eval_rppsurf(p,u,v) = p(u,v)                           *)        
        
(**********************************************************************)        
function eval_rppsurf(p : rppsurf; 
                      u : real;    
                      v : real): vector;                
var            
   i,j : integer;                  
   top,temp : vector;              
   bot : real; 
   s,t : real; 
   svalue : vecmdmd;               
begin (* function eval_rppsurf *)  
with p do      
begin (* with p *)                 
               
      i :=  bracket(num_segs1,break_point1,u);          
      j :=  bracket(num_segs2,break_point2,v);          
               
      s := u - break_point1(.i.);  
      t := v - break_point2(.j.);  
               
      eval_psurf(rppatch(.i,j.).numerator,0,0,s,t,svalue); 
      top := svalue(.0,0.);        
      eval_psurf(rppatch(.i,j.).denominator,0,0,s,t,svalue);
      temp := svalue(.0,0.);       
      bot := temp.component(.1.);  
      eval_rppsurf := quotient(top,bot);                
               
end;  (* with p *)                 
end;  (* function eval_rppsurf *)  
(**********************************************************************)        
%page          
(**********************************************************************)        
        
(*      name:  eval_pp                                                *)        
        
(*  category:  evaluation                                             *)        
        
(*   purpose:  evaluates a piecewise vector-valued polynomial         *)        
(*             and its derivatives at a given parameter value         *)        
        
(*     input:  f - a piecewise vector-valued polynomial               *)        
(*             x - a parameter value                                  *)        
(*             d - the maximum derivative desired                     *)        
        
(*    output:  ppvalue(.i.) = i-th deriv of f at x                    *)        
(*             for i = 0 to d                                         *)        
        
(**********************************************************************)        
               
procedure eval_pp(f : pp;          
                  x : real;        
                  d : integer;     
        var ppvalue : vecmd);      
var            
   i : integer;
   s : real;   
begin (* procedure eval_pp *)      
   i := bracket(f.num_segs,f.break_point,x);            
   with f do   
   begin (* with f *)              
      s := x - break_point(.i.);   
      eval_poly(segment(.i.),s,d,ppvalue);              
   end;  (* with f *)              
end;  (* procedure eval_pp *)      
               
(**********************************************************************)        
%page          
(**********************************************************************)        
        
(*      name:  eval_rpp                                               *)        
        
(*  category:  evaluation                                             *)        
        
(*   purpose:  evaluates a vector-valued rational piecewise           *)        
(*             polynomial at a given parameter value                  *)        
        
(*     input:  rpp - a vector-valued rational piecewise polynomial    *)        
(*             x - a parameter value                                  *)        
        
(*    output:  eval_rpp(rpp,x) = rpp(x)                               *)        
        
(**********************************************************************)        
               
function eval_rpp(rpp : rational_pp;
                    x : real):vector;
var            
   i,j : integer;                  
   ppvalue : vecmd;                
       top : vector;               
       bot : real;                 
         s : real;                 
               
begin (* function eval_rpp *)      
with rpp do    
begin (* with rpp *)               
   i := bracket(num_segs,break_point,x);                
   s := x - break_point(.i.);      
               
   with rpp_segment(.i.) do        
   begin (* with rpp_segment(.i.) *)
      eval_poly(numerator,s,0,ppvalue);                 
      top := ppvalue(.0.);         
      eval_poly(denominator,s,0,ppvalue);               
      bot := ppvalue(.0.).component(.1.);               
   end;  (* with rpp_segment(.i.) *)
               
   eval_rpp := quotient(top,bot);  
               
end;  (* with rpp *)               
end;  (* function eval_rpp *)      
               
(**********************************************************************)        
%page          
(**********************************************************************)        
        
(*      name:  knoteval                                               *)        
        
(*  category:  evaluation                                             *)        
        
(*   purpose:  finds the left and right derivatives of prescribed     *)        
(*             orders of a piecewise polynomial.  derivatives are     *)        
(*             taken at a prescribed breakpoint.  this is done in     *)        
(*             order to find the degree of continuity at that         *)        
(*             breakpoint.  knowing the degree of continuity there    *)        
(*             will enable us to assign the appropriate knot          *)        
(*             multiplicity                                           *)        
        
(*     input:  f - a piecewise polynomial                             *)        
(*             i - the index of the breakpoint                        *)        
(*             n - the highest order derivative desired               *)        
        
(*    output:  right_deriv(j) -  the j-th right derivative of f       *)        
(*                               at the breakpoint    (j = 0 to n)    *)        
(*             left_deriv(j)  -  the j-th left derivative of f        *)        
(*                               at the breakpoint    (j = 0 to n)    *)        
        
(**********************************************************************)        
               
procedure knoteval (f : pp;        
                    i : integer;   
                    n : integer;   
      var right_deriv : vecmd;     
      var  left_deriv : vecmd);    
var            
   s : real;   
   j : integer;
begin (* procedure knoteval *)     
   with f do   
   begin (* evaluating derivatives *)
      s := break_point(.i.) - break_point(.i-1.);       
      eval_poly(segment(.i.),0.0,n,right_deriv);        
      eval_poly(segment(.i-1.),s,n,left_deriv);         
   end;  (* evaluating derivatives *)
end;  (* procedure knoteval *)     
(**********************************************************************)        
%page          
(**********************************************************************)        
        
(*      name:  eval_bspline                                           *)        
        
(*  category:  evaluation                                             *)        
        
(*   purpose:  evaluate a bspline curve and its derivatives at        *)        
(*             a given parameter.  this routine finds the span        *)        
(*             and passes the appropriate portion of the spline's     *)        
(*             defining data to leeval, which performs the actual     *)        
(*             evaluation.                                            *)        
        
(*     input:  g - a bspline curve                                    *)        
(*             x - a parameter value                                  *)        
(*             d - the maximum derivative desired                     *)        
        
(*    output:  e_lee(i) = the i-th derivative of g at x               *)        
        
(**********************************************************************)        
               
procedure eval_bspline(g : bspline;
                       x : real;   
                       d : integer;
              var e_lee  : vecmd); 
var            
   i,k : integer;                  
   knot       : real44;            
   coeff      : vec44;             
               
begin (* procedure eval_bspline *) 
      i := knot_interval(g.number_of_knots,g.knot,x);   
/*             
      if idebug = 1 then           
      writeln(debug,'in eval_bspline knot_interval = ',i:1);
*/             
      for k := 0 to 2*g.degree + 1 do
      begin (* aligining the sequences *)               
        knot(.k.) := g.knot(.i - g.degree + k.);        
        coeff(.k.) := g. coeff(.i - g.degree + k.);     
      end;  (* aligning the sequences *)                
               
   leeval(g.degree,                
          knot,
          coeff,
          x,   
          d,   
          e_lee);                  
               
               
end;  (* procedure eval_bspline *) 
(**********************************************************************)        
%page          
(**********************************************************************)        
        
(*      name:  eval_rbs                                               *)        
        
(*  category:  evaluation                                             *)        
        
(*   purpose:  evaluates the value (not derivatives) of a             *)        
(*             rational bspline function at a given parameter value.  *)        
        
(*     input:  rbs - a rational bspline function                      *)        
(*               x - the parameter value                              *)        
        
(*    output:  eval_rbs(rbs,x) = rbs(x)                               *)        
        
(**********************************************************************)        
               
function eval_rbs(rbs : rational_bspline;               
                    x : real):vector;
var            
e_lee : vecmd; 
   top : vector;
   bot : real; 
begin (* function eval_rbs *)      
   with rbs do 
   begin (* with rbs *)            
      eval_bspline(numerator,x,0,e_lee);                
      top := e_lee(.0.);           
      eval_bspline(denominator,x,0,e_lee);              
      bot := e_lee(.0.).component(.1.);                 
   end;  (* with rbs *)            
               
   eval_rbs := quotient(top,bot);  
               
end;  (* function eval_rbs *)      
(**********************************************************************)        
%page          
(**********************************************************************)        
        
(*      name:  eval_bsurf                                             *)        
        
(*  category:  evaluation                                             *)        
        
(*   purpose:  evaluate a tensor product b-spline surface             *)        
        
(*     input:  b - a tensor product b-spline surface                  *)        
(*            d1 - maximum derivative desired in 1st parameter        *)        
(*            d2 - maximum derivative desired in 2nd parameter        *)        
(*             s - 1st parameter                                      *)        
(*             t - 2nd parameter                                      *)        
        
(*    output:  svalue(i,j) = the (i,j)-th derivative of b at (s,t)    *)        
        
(**********************************************************************)        
procedure eval_bsurf(b : bsurf;    
                    d1 : integer;  
                    d2 : integer;  
                     s : real;     
                     t : real;     
            var svalue : vecmdmd); 
var            
   i,j,k : integer;                
   f,g : bspline;                  
   cvalue : vecmd;                 
     q : vecmdmd;                  
begin (* procedure eval_bsurf *)   
               
   g.degree := b.degree2;          
   g.dim := b.dim;                 
   g.number_of_knots := b.num_of_knots2;                
   g.knot := b.knot2;              
               
      f.degree := b.degree1;       
      f.dim := b.dim;              
      f.number_of_knots := b.num_of_knots1;             
      f.knot := b.knot1;           
               
   for j := 1 to b.num_of_knots2 - b.degree2 - 1 do     
   begin (* j loop *)              
      for i:= 1 to b.num_of_knots1 - b.degree1 - 1 do   
      f.coeff(.i.) := b.coeff(.i,j.);
               
      eval_bspline(f,s,d1,cvalue); 
               
      for k := 0 to d1 do          

      q(.j,k.) := cvalue(.k.);     
   end;  (* j loop *)              
               
   for i := 0 to d1 do             
   begin (* i loop *)              
      for k := 1 to b.num_of_knots2 - b.degree2 - 1 do  
      g.coeff(.k.) := q(.k,i.);    
               
      eval_bspline(g,t,d2,cvalue); 
               
      for j := 0 to d2 do          
      svalue(.i,j.) := cvalue(.j.);
               
   end;  (* i loop *)              
               
end;  (* procedure eval_bsurf *)   
(**********************************************************************)        
%page          
(**********************************************************************)        
        
(*      name:  eval_rbsurf                                            *)        
        
(*  category:  evaluation                                             *)        
        
(*   purpose:  evaluate a tensor product rational b-spline surface    *)        
(*             function value only.  no derivatives                   *)        
        
(*     input:  rb - a tensor product rational b-spline surface        *)        
(*              s - the 1st parameter                                 *)        
(*              t - the 2nd parameter                                 *)        
        
(*    output:   eval_rbsurf(rb,s,t) = rb(s,t)                         *)        
        
(**********************************************************************)        
function eval_rbsurf(rb : rbsurf;  
                      s : real;    
                      t : real): vector;                
var            
   svalue : vecmdmd;               
   top,botvec : vector;            
   bot : real; 
begin (* function eval_rbsurf *)   
with rb do     
   begin (* with rb *)             
      eval_bsurf(numerator,0,0,s,t,svalue);             
      top := svalue(.0,0.);        
               
      eval_bsurf(denominator,0,0,s,t,svalue);           
      botvec := svalue(.0,0.);     
      bot := botvec.component(.1.);
               
      eval_rbsurf := quotient(top,bot);                 
   end;  (* with rb *)             
end;  (* function eval_rbsurf *)   
(**********************************************************************)        
%page          
(**********************************************************************)        
        
(*      name:  determine_continuity                                   *)        
        
(*  category:  knot finding                                           *)        
        
(*   purpose:  given a piecewise polynomial and a tolerance,          *)        
(*             find the level of derivative continuity (to within     *)        
(*             tolerance) at each breakpoint.  this information       *)        
(*             will be used to determine the knot multiplicities.     *)        
        
(*     input:  f - a piecewise polynomial                             *)        
(*             tol - tolerance used to determine derivative           *)        
(*                   continuity                                       *)        
        
(*    output:  nu  - an array defined as follows.  let c(i) =         *)        
(*                   the level of derivative continuity at the i-th   *)        
(*                   breakpoint.  then nu(i) = c(i) + 1.              *)        
(*                   for instance, if nu(i) = 0, the function is      *)        
(*                   not even continuous at the i-th breakpoint.      *)        
(*                   if nu(i) = 2, the function and its first         *)        
(*                   derivative are continuous there.                 *)        
        
(**********************************************************************)        
               
procedure determine_continuity(f:pp;
                               tol:real;                
                           var nu:integer44);           
var            
   left_deriv,right_deriv:vecmd;   
   i,j,maxdeg:integer;             
   gap:real;   
   x : real;   
   m,n : integer;                  
   temp : real;
   factor : real;                  
   continuous : boolean;           
begin (* procedure determine_continuity *)              
   with f do   
   begin (* with f *)              
      maxdeg := 0;                 
               
 
      for i:= 1 to num_segs do     
      maxdeg := max(maxdeg,segment(.i.).degree);        
                
for i:=2 to num_segs do            
      begin (* finding nu at the i-th breakpoint *)     
         knoteval(f,i,maxdeg,right_deriv,left_deriv);   
         nu(.i.):=0;                    
         j:=-1;                         
         repeat (* until (j=maxdeg) or (continuous = false) *)                  
         j := j+1;                      
         gap := dist(right_deriv(.j.),left_deriv(.j.));    
                                        
(* the following is a heuristic method for determining continuity *)            
      temp := max(norm(right_deriv(.j.)),norm(left_deriv(.j.)));                
      if temp <> 0.0 then               
      factor := gap/(temp*machine_epsilon)                 
      else factor := 0.0;               
                                        
      if (gap < tol) or (factor < power(10,j))             
      then                              
      continuous := true                
      else                              
      continuous := false;              
                                        
if idebug = 1 then                      
begin                                   
writeln(report);                        
writeln(report,' breakpoint = ',i:1,' deriv = ',j:1,       
' gap = ',gap,'   factor =',factor:5:2);
writeln(report,' continuous =',continuous);                
end;                                    
                                        
         if continuous then nu(.i.) := nu(.i.) + 1;        
         until (j=maxdeg) or (continuous = false);         
end;  (* finding nu at the i-th breakpoint *)              
         nu(.1.):=0;                    
         nu(.num_segs + 1.) := 0;       
   end;  (* with f *)                   
end; (* procedure determine_continuity *)                  
(**********************************************************************)        
%page                                   
(**********************************************************************)        
        
(*      name:  find_knots                                             *)        
        
(*  category:  knot_finding                                           *)        
        
(*   purpose:  given a piecewise polynomial, find the knot values     *)        
(*             of the corresponding b-spline curve, with each knot    *)        
(*             value given the appropriate multiplicity.              *)        
        
(*     input:  f - a piecewise polynomial                             *)        
(*             order - degree + 1                                     *)        
(*             tol - tolerance used to determine level of derivative  *)        
(*                   continuity.                                      *)        
        
(*    output:  number_of_knots - the total number of knots, with each *)        
(*                               knot counted as many times as its    *)        
(*                               multiplicity.                        *)        
(*             t - the knot sequence.  the first and last knots are   *)        
(*                 always given multiplicity = degree + 1             *)        
        
(**********************************************************************)        
                                        
procedure find_knots(f : pp;            
                     order: integer;    
                     tol : real;        
                     var number_of_knots : integer;        
                     var t : real44);   
   var                                  
      i,j,k : integer;                  
      nu : integer44;                   
      knot_multiplicity : array(.1..44.) of integer;       
                                        
begin (* procedure find_knots *)        
   determine_continuity(f,tol,nu);      
   with f do                            
   begin (* determining knots *)        
      for i:= 1 to num_segs+1 do        
      begin (* assigning knot_multiplicity(.i.)*)          
      knot_multiplicity(.i.):=order - nu(.i.);             
      end;  (* assigning knot_multiplicity(.i.) *)         
                                        
      k := 0;                           
                                        
      for i := 1 to num_segs+1 do       
      for j := 1 to knot_multiplicity(.i.) do              
      begin                             
         k := k + 1;                    
         t(.k.) := break_point(.i.);    
      end;                              
                                        
   number_of_knots := k;                
                                        
    end; (* determining knots *)        
                                        
end; (* procdure find_knots *)          
                                        
(**********************************************************************)        
%page                                   
(**********************************************************************)        
        
(*      name:  find_knots2                                            *)        
        
(*  category:  knot finding                                           *)        
        
(*   purpose:  given a piecewise polynomial of two variables,         *)        
(*             knot sequences are found in each parameter which       *)        
(*             are for the corresponding piecewise bezier             *)        
(*             surface represented in b-spline form.  knots that      *)        
(*             are not really needed (due to derivative continuity)   *)        
(*             can be removed by the procedure remove_all_surf_knots. *)        
        
(*     input:  p - a piecewise polynomial surface                     *)        
        
(*    output:  n1 - number of knots in 1st parameter                  *)        
(*             n2 - number of knots in 2nd parameter                  *)        
(*             t1 - 1st knot sequence.  each breakpoint in 1st param  *)        
(*                  is given multiplicity degree1 + 1                 *)        
(*             t2 - 2nd knot sequence.  each breakpoint in 2nd param  *)        
(*                  is given multiplicity degree2 + 1                 *)        
        
(**********************************************************************)        
procedure find_knots2(p : ppsurf;       
                 var n1 : integer;      
                 var n2 : integer;      
                 var t1 : real44;       
                 var t2 : real44);      
var                                     
   i,j,k : integer;                     
begin (* procedure find_knots2 *)       
with p do                               
begin (* with p *)                      
                                        
   n1 := (num_segs1 + 1)*(degree1 + 1); 
   n2 := (num_segs2 + 1)*(degree2 + 1); 
                                        
   k := 0;                              
   for i := 1 to num_segs1 + 1 do       
   for j := 1 to degree1 + 1 do         
   begin                                
      k := k + 1;                       
      t1(.k.) := break_point1(.i.);     
   end;                                 
                                        
   k := 0;                              
   for i := 1 to num_segs2 + 1 do       
   for j := 1 to degree2 + 1 do         
   begin                                
      k := k + 1;                       
      t2(.k.) := break_point2(.i.);     
   end;                                 
                                        
end;  (* with p *)                      
end;  (* procedure find_knots2 *)       
(**********************************************************************)        
%page                                   
(**********************************************************************)        
        
(*      name:  find_distinct_knots                                    *)        
        
(*  category:  knot finding                                           *)        
        
(*   purpose:  given a knot sequence, find the distinct knots.        *)        
(*             these will be the breakpoints of the corresponding     *)        
(*             piecewise polynomial.                                  *)        
        
(*     input:  number_of_knots - the number of knots (not surprising) *)        
(*             knot - the knot sequence                               *)        
        
(*    output:  n_distinct_knots - the number of distinct knots        *)        
(*               distinct_knot  - the sequence of distinct knots      *)        
        
(**********************************************************************)        
                                        
procedure find_distinct_knots(number_of_knots : integer;   
                                         knot : real44;    
                         var n_distinct_knots : integer;   
                            var distinct_knot : real44);   
var                                     
   i : integer;                         
begin (* procedure find_distinct_knots *)                  
   n_distinct_knots := 1;               
   distinct_knot(.1.) := knot(.1.);     
   for i := 1 to number_of_knots - 1 do 
   if knot(.i.) < knot(.i+1.) then      
   begin (* updating *)                 
      n_distinct_knots := n_distinct_knots + 1;            
      distinct_knot(.n_distinct_knots.) := knot(.i+1.);    
   end;  (* updating *)                 
end;  (* procedure find_distinct_knots *)                  
                                        
(**********************************************************************)        
%page                                   
(**********************************************************************)        
        
(*      name:  add_knot                                               *)        

*)        
(*  category:  knot manipulation                                      *)        

(*   purpose:  given a b-spline curve and a parameter value           *)        
(*             make that parameter value another knot.  adding this   *)        
(*             knot does not change the parameterization of the       *)        
(*             curve.  however, some new b-spline coefficients        *)        
(*             have to be calculated to reflect the presence of the   *)        
(*             new knot.   this procedure follows the algorithm       *)        
(*             given in the paper by w. bohm, "inserting new knots    *)        
(*             into b-spline curves", computer-aided design, vol. 12, *)        
(*             pp. 199-201                                            *)        
        
(*     input:  f - a b-spline curve                                   *)        
(*             x - the new knot to be added                           *)        

(*    output:  g - the same b-spline curve expressed in terms of      *)        
(*                 the new knot sequence                              *)        
(**********************************************************************)        
procedure add_knot(f : bspline;             
 x : real;                
               var g : bspline);            
var     
   i,j,m : integer;       
   alpha,beta : real;     
   q1,q2 : vector;        
        
begin (* procedure add_knot *)              
        
   (* check to make sure that the knot to be added is within the *)             
   (* original parameter range.  if it is not, issue an error message *)        
   (* and set the output bspline equal to the input bspline *)                  
   with f do              
   if (x < knot(.1.)) or (x > knot(.number_of_knots.)) then   
   begin (* abnormal case *)                
      writeln(' the knot to be added is outside the parameter range ');         
      writeln(' therefore, knot not added.');                 
      g := f;             
   end   (* abnormal case *)                
   else
   begin (* normal case *)                  
   g.degree := f.degree; 
   g.dim := f.dim;        
   g.number_of_knots := f.number_of_knots + 1;                
        
   m := f.degree;         
   i := knot_interval(f.number_of_knots,f.knot,x);            
        
   for j := 1 to i do     
   g.knot(.j.) := f.knot(.j.);              
        
   g.knot(.i+1.) := x;    
        
   for j := i + 2 to g.number_of_knots do   
   g.knot(.j.) := f.knot(.j - 1.);          
        
   for j := 1 to i - m do
   g.coeff(.j.) := f.coeff(.j.);            
        
   for j := i - m + 1 to i do               
   begin (* defining new coefficients in terms of old *)      
      alpha := (f.knot(.j+m.)-x)/(f.knot(.j+m.)-f.knot(.j.)); 
      beta := 1.0 - alpha;                  
      q1 := product(alpha,f.coeff(.j-1.));  
      q2 := product(beta,f.coeff(.j.));     
      g.coeff(.j.) := sum(q1,q2);           
   end;  (* defining new coefficients in terms of old *)      
        
   for j := i + 1 to g.number_of_knots - m - 1 do             
   g.coeff(.j.) := f.coeff(.j - 1.);        
   end;  (* normal case *)                  
        
end;  (* procedure add_knot *)              
(**********************************************************************)        
(**********************************************************************)        
        
(*      name:  add_surf_knot      *)        
        
(*  category:  knot manipulation  *)        
        
(*   purpose:  given a b-spline surface, add a knot in the first      *)        
(*             or second parameter.  the corresponding curve          *)        
(*             algorithm add_knot is repeatedly invoked.              *)        
        
(*     input:  b - a b-spline surface               *)        
(*             param - 1 = add knot in 1st param    *)        
(*   2 = add knot in 2nd param    *)        
(*             x - knot value to be added           *)        
        
(*    output:  a - the resulting b-spline surface with the new        *)        
(*                 knot added.    *)        
        
(**********************************************************************)        
procedure add_surf_knot(b : bsurf;          
  param : integer;        
      x : real;           
  var a : bsurf);         
var     
   i,j,c1,c2 : integer;   
   d,e : bspline;         
begin (* procedure add_surf_knot *)         
        
   c1 := b.num_of_knots1 - b.degree1 - 1;   
   c2 := b.num_of_knots2 - b.degree2 - 1;   
        
   a := b; (* to start with, set everything equal *)          
        
   case param of          
   1 : 
   begin (* adding knot in 1st parameter *)
      d.degree := b.degree1;                
      d.dim :=  b.dim;    
      for j := 1 to c2 do
      begin (* j loop *) 
         d.number_of_knots := b.num_of_knots1;                
         d.knot := b.knot1;                 
         for i := 1 to c1 do d.coeff(.i.) := b.coeff(.i,j.);  
         add_knot(d,x,e);
         for i := 1 to c1 + 1 do a.coeff(.i,j.) := e.coeff(.i.);                
      end;  (* j loop *) 
      a.num_of_knots1 := e.number_of_knots;
      a.knot1 := e.knot; 
   end;  (* adding knot in 1st parameter *)
   2 : 
   begin (* adding knot in 2nd parameter *)
      d.degree := b.degree2;                
      d.dim :=  b.dim;    
      for i := 1 to c1 do
      begin (* i loop *) 
         d.number_of_knots := b.num_of_knots2;                
         d.knot := b.knot2;                 
         for j := 1 to c2 do d.coeff(.j.) := b.coeff(.i,j.);  
         add_knot(d,x,e);
         for j := 1 to c2 + 1 do a.coeff(.i,j.) := e.coeff(.j.);                
      end;  (* i loop *) 
      a.num_of_knots2 := e.number_of_knots;
      a.knot2 := e.knot; 
   end;  (* adding knot in 2nd parameter *)
   otherwise              
   begin (* error processing *)             
      writeln(' incorrect input param =',param:1,'in add_surf_knot.');          
      writeln(' output surface set equal to input surface.'); 
      a := b;             
   end;  (* error processing *)             
   end;  (*case statment *)                 
end;  (* procedure add_surf_knot *)         
        
(**********************************************************************)        
(**********************************************************************)        
        
(*      name:  remove_knot        *)        
        
(*  category:  knot manipulation  *)        
        
(*   purpose:  remove a knot from a b-spline curve if it can be       *)        
(*             done.  bohm's algorithm (see procedure add_knot)       *)        
(*             is reversed.  at a certain step in this reversed       *)        
(*             algorithm one of the new coefficients can be           *)        
(*             expressed in two different ways.  if these two         *)        
(*             representations are equal (to within tolerance)        *)        
(*             we can remove the knot.              *)        
        
(*     input:  g - a b-spline curve                 *)        
(*             i - index of knot to remove          *)        
(*           tol - tolerance used to decide whether knot can          *)        
(*                 be removed.    *)        
        
(*    output:  removed = true if knot was removed, false otherwise.   *)        
(*             f - the b-spline curve with a knot removed             *)        
(*                 (use f only if removed = true)   *)        
        
(**********************************************************************)        
procedure remove_knot(g : bspline;          
    i : integer;          
  tol : real;             
            var removed : boolean;          
                  var f : bspline);         
var     
   j,l,d : integer;       
   alpha,beta : real;     
   x : real;              
   v1,v2 : vector;        
   p,q : vec44;           
   t : real44;            
   k : integer;           
   max_coeff_norm : real;
   distance : real;       
factor : real;            
begin (* procedure remove_knot *)           
        
   d := g.degree;         
(* check if the index of knot to remove is a legal one *)     
if (i <= d + 1) or (i >= g.number_of_knots - d) then          
begin (* error processing *)                
   removed := false;      
   writeln(output,i:2,' is not a valid index for knot removal');                
   writeln(' index must be > ',d+1:1,' and < ',               
   g.number_of_knots - d:1);                
end   (* error processing *)                
else    
begin (* normal processing *)               
(* check if this is a multiple knot *)      
(* if it is, reset i to be the lowest index for the multiple knot *)            
   repeat (* until g.knot(.i-1.) < g.knot(.i.) *)             
   if g.knot(.i-1.) = g.knot(.i.) then i := i - 1;            
   until g.knot(.i-1.) < g.knot(.i.);       
        
   for j := 1 to i-1 do   
   t(.j.) := g.knot(.j.);
        
   x := g.knot(.i.);      
        
   for j := i to g.number_of_knots - 1 do   
   t(.j.) := g.knot(.j+1.);                 
        
   for j := 1 to g.number_of_knots - d - 1 do                 
   q(.j.) := g.coeff(.j.);                  
        
   l := i - 1;            
   for j := 1 to l - d do
   p(.j.) := q(.j.);      
        
   for j := l - d + 1 to l do               
   begin (* backing out p's in terms of q's *)                
      alpha := (t(.j+d.) - x)/(t(.j+d.)-t(.j.));              
      beta := 1.0 - alpha;                  
      v1 := product(alpha,p(.j-1.));        
      v2 := difference(q(.j.),v1);          
      p(.j.) := quotient(v2,beta);          
   end;  (* backing out p's in terms of q's *)                
        
if idebug = 1 then        
   begin                  
   max_coeff_norm := 0.0;
   for k := 1 to l do     
   max_coeff_norm := max(max_coeff_norm,norm(q(.k.)));        
   distance := dist(p(.l.),q(.l+1.));       
   if max_coeff_norm > 0.0 then             
   factor := distance/(max_coeff_norm*machine_epsilon*d)      
   else factor := 0.0;    
   writeln(output,' in remove_knot dist(p(.l.),q(.l+1.) =',   
          distance);      
   writeln(output,' factor =',factor:6:2);  
   writeln(debug,' in remove_knot dist(p(.l.),q(.l+1.) =',    
          distance);      
   writeln(debug,' factor =',factor:6:2);   
   end;
        
   if dist(p(.l.),q(.l+1.)) < tol then      
   begin (* completing the process of removing a knot *)      
      for j := l to g.number_of_knots - d - 2 do              
      p(.j.) := q(.j+1.);
        
      for j := l downto l - (d div 2) + 1 do                  
      begin (* backing out half the coeffs from the top *)    
         alpha := (t(.j+d.) - x)/(t(.j+d.)-t(.j.));           
         if alpha = 0.0 then leave;         
         beta := 1.0 - alpha;               
         v1 := product(beta,p(.j.));        
         v2 := difference(q(.j.),v1);       
         p(.j-1.) := quotient(v2,alpha);    
      end;  (* backing out half the coeffs from the top *)    
        
      f.degree := g.degree;                 
      f.dim := g.dim;     
      f.number_of_knots := g.number_of_knots - 1;             
        
      for j := 1 to f.number_of_knots do    
      f.knot(.j.) := t(.j.);                
        
      for j := 1 to f.number_of_knots - f.degree - 1 do       
      f.coeff(.j.) := p(.j.);               
        
      removed := true;    
   end   (* completing the process of removing a knot *)      
   else
      removed := false;   
end;  (* normal processing *)               
end;  (* procedure remove_knot *)           
(**********************************************************************)        
(**********************************************************************)        
        
(*      name:  remove_surf_knot   *)        
        
(*  category:  knot manipulation  *)        
        
(*   purpose:  remove a knot from one of the parameters of a b-spline *)        
(*             surface.  if the knot can be removed from all the      *)        
(*             corresponding curves, it is removed from the surface.  *)        
        
(*     input:  b - a b-spline surface               *)        
(*             param = 1 (remove knot in 1st param) *)        
(*   2 (remove knot in 2nd param) *)        
(*             r - index of knot to be removed      *)        
(*             tol - tolerance to be used to determine whether knot   *)        
(* can be removed                 *)        
        
(*    output:   removed = true if knot can be removed                 *)        
(*    = false otherwise           *)        
(*  a - the surface with the knot removed, use      *)        
(*      only if remove = true     *)        
        
(**********************************************************************)        
procedure remove_surf_knot(b : bsurf;       
     param : integer;     
         r : integer;     
       tol : real;        
                 var removed : boolean;     
     var a : bsurf);      
var     
   i,j,c1,c2 : integer;   
   f,g : bspline;         
begin (* procedure remove_surf_knot *)      
        
   c1 := b.num_of_knots1 - b.degree1 - 1;   
   c2 := b.num_of_knots2 - b.degree2 - 1;   
        
   a := b; (* to start with, set everything equal *)          
        
   case param of          
   1 : 
   begin (* removing knot in 1st parameter *)                 
      f.degree := b.degree1;                
      f.dim :=  b.dim;    
      for j := 1 to c2 do
      begin (* j loop *) 
         f.number_of_knots := b.num_of_knots1;                
         f.knot := b.knot1;                 
         for i := 1 to c1 do f.coeff(.i.) := b.coeff(.i,j.);  
         remove_knot(f,r,tol,removed,g);    
         if removed then 
         for i := 1 to c1 - 1 do a.coeff(.i,j.) := g.coeff(.i.)                 
         else             
         leave;           
      end;  (* j loop *) 
      a.num_of_knots1 := g.number_of_knots;
      a.knot1 := g.knot; 
   end;  (* removing knot in 1st parameter *)                 
   2 : 
   begin (* removing knot in 2nd parameter *)                 
      f.degree := b.degree2;                
      f.dim :=  b.dim;    
      for i := 1 to c1 do
      begin (* i loop *) 
         f.number_of_knots := b.num_of_knots2;                
         f.knot := b.knot2;                 
         for j := 1 to c2 do f.coeff(.j.) := b.coeff(.i,j.);  
         remove_knot(f,r,tol,removed,g);    
         if removed then 
         for j := 1 to c2 - 1 do a.coeff(.i,j.) := g.coeff(.j.)                 
         else             
         leave;           
      end;  (* i loop *) 
      a.num_of_knots2 := g.number_of_knots;
      a.knot2 := g.knot; 
   end;  (* removing knot in 2nd parameter *)                 
   end;  (*case statment *)                 
end;  (* procedure remove_surf_knot *)      
        
(**********************************************************************)        
(**********************************************************************)        
        
(*      name:  remove_all_surf_knots                *)        
        
(*  category:  knot manipulation  *)        
        
(*   purpose:  remove all knots that can be removed (to within        *)        
(*             tolerance) from a given surface in both parameters.    *)        
        
(*     input:  b - a b-spline surface               *)        
(*             tol - the tolerance used to determine whether a knot   *)        
(* can be removed.                *)        
        
(*    output:  a - the b-spline surface with all possible knots       *)        
(*                 removed.       *)        
        
(**********************************************************************)        
procedure remove_all_surf_knots(b : bsurf;  
            tol : real;   
          var a : bsurf);
var     
   removed : boolean;     
   temp : bsurf;          
   current_knot_index : integer;            
   number_left : integer;
begin (* procedure remove_all_surf_knots *)
   if idebug = 1 then     
   begin                  
   writeln(output,'entering remove_all_surf_knots');          
   writeln(debug,'entering remove_all_surf_knots');           
   end;
        
(*--------------------------------------------------------------------*)        
   number_left := b.num_of_knots1 - 2*(b.degree1 + 1);        
   current_knot_index := b.degree1 + 2;     
   while number_left > 0 do                 
   begin (* trying to remove knots in 1st parameter *)        
      remove_surf_knot(b,1,current_knot_index,tol,removed,temp);                
if idebug = 1 then        
begin (* debug output *) 
  writeln(output,'removing knots in 1st parameter: ');        
  writeln(output,'current_knot_index =',current_knot_index);  
  writeln(output,'number_left =',number_left);                
  writeln(output,'removed =',removed);      
  writeln(debug,'removing knots in 1st parameter: ');         
  writeln(debug,'current_knot_index =',current_knot_index);   
  writeln(debug,'number_left =',number_left);                 
  writeln(debug,'removed =',removed);       
end;  (* debug output *) 
      if removed then     
      begin (* successful knot removal update *)              
         a := temp;       
         b := a;          
         number_left := number_left - 1     
         (* current_knot_index stays the same *)              
      end   (* successful knot removal update *)              
      else                
      begin(* unsuccessful knot removal update *)             
         i := 0;          
         repeat (* until b.knot1(.current_knot_index + i.)    
     > b.knot1(.current_knot_index.) *)     
         i := i + 1;      
         number_left := number_left - 1;    
         until b.knot1(.current_knot_index + i.)              
             > b.knot1(.current_knot_index.);                 
         current_knot_index := current_knot_index + i;        
      end; (* unsuccessful knot removal update *)             
   end;  (* trying to remove knots in 1st parameter *)        
(*--------------------------------------------------------------------*)        
   number_left := b.num_of_knots2 - 2*(b.degree2 + 1);        
   current_knot_index := b.degree2 + 2;     
   while number_left > 0 do                 
   begin (* trying to remove knots in 2nd parameter *)        
      remove_surf_knot(b,2,current_knot_index,tol,removed,temp);                
if idebug = 1 then        
begin (* debug output *) 
  writeln(output,'removing knots in 2nd parameter: ');        
  writeln(output,'current_knot_index =',current_knot_index);  
  writeln(output,'number_left =',number_left);                
  writeln(output,'removed =',removed);      
  writeln(debug,'removing knots in 2nd parameter: ');         
  writeln(debug,'current_knot_index =',current_knot_index);   
  writeln(debug,'number_left =',number_left);                 
  writeln(debug,'removed =',removed);       
end;  (* debug output *) 
      if removed then     
      begin (* successful knot removal update *)              
         a := temp;       
         b := a;          
         number_left := number_left - 1     
         (* current_knot_index stays the same *)              
      end   (* successful knot removal update *)              
      else                
      begin(* unsuccessful knot removal update *)             
         i := 0;          
         repeat (* until b.knot2(.current_knot_index + i.)    
     > b.knot2(.current_knot_index.) *)     
         i := i + 1;      
         number_left := number_left - 1;    
         until b.knot2(.current_knot_index + i.)              
             > b.knot2(.current_knot_index.);                 
         current_knot_index := current_knot_index + i;        
      end; (* unsuccessful knot removal update *)             
   end;  (* trying to remove knots in 2nd parameter *)        
   if idebug = 1 then     
   begin                  
   writeln(output,'leaving remove_all_surf_knots.');          
   writeln(debug,'leaving remove_all_surf_knots.');           
   end;
end;  (* procedure remove_all_surf_knots *)
(**********************************************************************)        

(**********************************************************************)        
        
(*      name:  poly_to_rpoly      *)        
        
(*  category:  conversion between like entities.    *)        
        
(*   purpose:  given a vector-valued parametric polynomial curve,     *)        
(*             create the corresponding vector-valued rational        *)        
(*             parametric polynomial curve by using the last          *)        
(*             component of the input coefficients as the             *)        
(*             coefficient of the denominator.      *)        
        
(*     input:  p - a vector valued parametric polynomial curve.       *)        
        
(*    output:  rp - the corresponding vector valued (one less         *)        
(*                  dimension) rational parametric polynomial         *)        
(*                  curve.        *)        
        
(**********************************************************************)        
procedure poly_to_rpoly(p : poly;           
  var rp : rational_poly);                  
var     
i : integer;              
begin (* procedure poly_to_rpoly *)         
   rp.degree := p.degree;
   rp.dim := p.dim - 1;   
        
   for i := 0 to p.degree do                
   begin                  
      rp.top_coeff(.i.) := project(p.coeff(.i.));             
      rp.bot_coeff(.i.) := p.coeff(.i.).component(.p.dim.);   
   end;
end;  (* procedure poly_to_rpoly *)         
(**********************************************************************)        

(**********************************************************************)        
        
(*      name:  rpoly_to_poly      *)        
        
(*  category:  conversion between like entities.    *)        
        
(*   purpose:  given a vector-valued rational parametric polynomial   *)        
(*             curve, create the corresponding vector-valued          *)        
(*             parametric polynomial curve by setting the last        *)        
(*             component of the output coefficients equal to the      *)        
(*             coefficient of the input denominator *)        
        
(*     input:  rp - a vector-valued rational parametric polynomial    *)        
(*                  curve.        *)        
        
(*    output:   p - the corresponding vector valued (one more         *)        
(*                  dimension) parametric polynomial                  *)        
(*                  curve.        *)        
        
(**********************************************************************)        
procedure rpoly_to_poly(rp : rational_poly;
   var p : poly);         
var     
i : integer;              
begin (* procedure rpoly_to_poly *)         
   p.degree := rp.degree;
   p.dim := rp.dim + 1;   
        
   for i:= 0 to p.degree do                 
   p.coeff(.i.) := append(rp.top_coeff(.i.),rp.bot_coeff(.i.));                 
end;  (* procedure rpoly_to_poly *)         
(**********************************************************************)        

(**********************************************************************)        
        
(*      name:  pp_to_rpp          *)        
        
(*  category:  conversion between like entities.    *)        
        
(*   purpose:  given a vector-valued piecewise polynomial curve,      *)        
(*             create the corresponding vector-valued rational        *)        
(*             piecewise polynomial curve by using the last           *)        
(*             component of the input coefficients as the             *)        
(*             coefficient of the denominator.      *)        
        
(*     input:  f - a vector valued piecewise polynomial curve.        *)        
        
(*    output:  rf - the corresponding vector valued (one less         *)        
(*                  dimension) rational piecewise polynomial          *)        
(*                  curve.        *)        
        
(**********************************************************************)        
procedure pp_to_rpp(f : pp;                 
              var  rf : rational_pp);       
var     
i : integer;              
begin (* procedure pp_to_rpp *)             
   rf.degree := f.degree;
   rf.num_segs := f.num_segs;               
   rf.dim := f.dim - 1;   
   rf.break_point := f.break_point;         
        
   for i := 1 to f.num_segs do              
   poly_to_rpoly(f.segment(.i.),rf.rpp_segment(.i.));         
        
end;  (* procedure pp_to_rpp *)             
(**********************************************************************)        

(**********************************************************************)        
        
(*      name:  rpp_to_pp          *)        
        
(*  category:  conversion between like entities.    *)        
        
(*   purpose:  given a vector-valued rational piecewise polynomial    *)        
(*             curve, create the corresponding vector-valued          *)        
(*             piecewise polynomial curve by setting the last         *)        
(*             component of the output coefficients equal to the      *)        
(*             coefficient of the output denominator                  *)        
        
(*     input:  rf - a vector-valued rational piecewise polynomial     *)        
(*                  curve.        *)        
        
(*    output:   f - the corresponding vector valued (one more         *)        
(*                  dimension) piecewise polynomial *)        
(*                  curve.        *)        
        
(**********************************************************************)        
procedure rpp_to_pp(rf : rational_pp;       
               var   f : pp);               
var     
i : integer;              
begin (* procedure rpp_to_pp *)             
  f.degree := rf.degree; 
  f.num_segs := rf.num_segs;                
  f.dim := rf.dim + 1;    
  f.break_point := rf.break_point;          
        
for i := 1 to f.num_segs do                 
   rpoly_to_poly(rf.rpp_segment(.i.),f.segment(.i.));         
end;  (* procedure rpp_to_pp *)             
(**********************************************************************)        

(**********************************************************************)        
        
(*      name:  bspline_to_rbspline                  *)        
        
(*  category:  conversion between like entities.    *)        
        
(*   purpose:  given a vector-valued bspline curve, *)        
(*             create the corresponding vector-valued rational        *)        
(*             bspline curve by interpreting the input                *)        
(*             curve as being in homogeneous form.  *)        
        
(*     input:  g - a bspline curve.                 *)        
        
(*    output:  rg - the corresponding vector valued (one less         *)        
(*                  dimension) rational bspline curve.                *)        
        
(**********************************************************************)        
procedure bspline_to_rbspline(g : bspline;  
       var rg : rational_bspline);          
var     
   i : integer;           
   temp : vector;         
begin (* procedure bspline_to_rbspline *)   
   rg.degree := g.degree;
   rg.dim := g.dim - 1;   
   rg.number_of_knots := g.number_of_knots;
   rg.knot := g.knot;     
        
for i := 1 to g.number_of_knots - g.degree - 1 do             
begin
   rg.weight(.i.) := g.coeff(.i.).component(.g.dim.);         
   temp := project(g.coeff(.i.));           
   rg.coeff(.i.) := quotient(temp,rg.weight(.i.));            
end;    
end;  (* procedure bspline_to_rbspline *)   
(**********************************************************************)        

(**********************************************************************)        
        
(*      name:  rbspline_to_bspline                  *)        
        
(*  category:  conversion between like entities.    *)        
        
(*   purpose:  given a vector-valued rational b-spline                *)        
(*             curve, create the corresponding vector-valued          *)        
(*             b-spline curve by expressing the input                 *)        
(*             rational b-spline in homogeneous form                  *)        
        
(*     input:  rg - a vector-valued rational b-spline curve           *)        
        
(*    output:   g - the corresponding vector valued (one more         *)        
(*                  dimension) parametric b-spline  *)        
(*                  curve.        *)        
        
(**********************************************************************)        
procedure rbspline_to_bspline(rg : rational_bspline;          
         var g : bspline);                  
var     
   i : integer;           
   temp : vector;         
begin (* procedure rbspline_to_bspline *)   
   g.degree := rg.degree;
   g.dim := rg.dim + 1;   
   g.number_of_knots := rg.number_of_knots;
   g.knot := rg.knot;     
        
   for i := 1 to g.number_of_knots - g.degree- 1 do           
   begin                  
      temp := product(rg.weight(.i.),rg.coeff(.i.));          
      g.coeff(.i.) := append(temp,rg.weight(.i.));            
   end;
        
end;  (* procedure rbspline_to_bspline *)   
(**********************************************************************)        

(**********************************************************************)        
        
(*      name:  psurf_to_rpsurf    *)        
        
(*  category:  conversion between like entities.    *)        
        
(*   purpose:  given a vector-valued parametric polynomial surface    *)        
(*             create the corresponding vector-valued rational        *)        
(*             parametric polynomial surface by using the last*)                
(*             component of the input coefficients as the             *)        
(*             coefficient of the denominator.      *)        
        
(*     input:  p - a vector valued parametric polynomial surface      *)        
        
(*    output:  rp - the corresponding vector valued (one less         *)        
(*                  dimension) rational parametric polynomial         *)        
(*                  surface       *)        
        
(**********************************************************************)        
procedure psurf_to_rpsurf(p : psurf;        
  var rp : rpsurf);    
var                                     
   i,j : integer;                       
begin (* procedure psurf_to_rpsurf *)   
   rp.degree1 := p.degree1;             
   rp.degree2 := p.degree2;             
   rp.dim := p.dim - 1;                 
                                        
   for i := 0 to p.degree1 do           
   for j := 0 to p.degree2 do           
   begin                                
      rp.top_coeff(.i,j.) := project(p.coeff(.i,j.));      
      rp.bot_coeff(.i,j.) := p.coeff(.i,j.).component(.p.dim.);                 
   end;                                 
end;  (* procedure psurf_to_rpsurf *)   
(**********************************************************************)        
                                   
(**********************************************************************)        
        
(*      name:  rpsurf_to_psurf                                        *)        
        
(*  category:  conversion between like entities.                      *)        
        
(*   purpose:  given a vector-valued rational parametric polynomial   *)        
(*             surface, create the corresponding vector-valued        *)        
(*             parametric polynomial surface by setting the last      *)        
(*             component of the output coefficients equal to the      *)        
(*             coefficient of the input denominator                   *)        
        
(*     input:  rp - a vector-valued rational parametric polynomial    *)        
(*                  surface.                                          *)        
        
(*    output:   p - the corresponding vector valued (one more         *)        
(*                  dimension) parametric polynomial                  *)        
(*                  surface.                                          *)        
        
(**********************************************************************)        
procedure rpsurf_to_psurf(rp : rpsurf;  
                     var p : psurf);    
var                                     
   i,j : integer;                       
begin (* procedure rpsurf_to_psurf *)   
   p.degree1 := rp.degree1;             
   p.degree2 := rp.degree2;             
   p.dim := rp.dim + 1;                 
                                        
   for i:= 0 to p.degree1 do            
   for j:= 0 to p.degree2 do            
   p.coeff(.i,j.) := append(rp.top_coeff(.i,j.),rp.bot_coeff(.i,j.));           
end;  (* procedure rpsurf_to_psurf *)   
(**********************************************************************)        
                                   
(**********************************************************************)        
        
(*      name:  ppsurf_to_rppsurf                                      *)        
        
(*  category:  conversion between like entities.                      *)        
        
(*   purpose:  given a vector-valued piecewise polynomial surface     *)        
(*             create the corresponding vector-valued rational        *)        
(*             piecewise polynomial surface by using the last         *)        
(*             component of the input coefficients as the             *)        
(*             coefficient of the denominator.                        *)        
        
(*     input:  f - a vector valued piecewise polynomial surface       *)        
        
(*    output:  rf - the corresponding vector valued (one less         *)        
(*                  dimension) rational piecewise polynomial          *)        
(*                  surface                                           *)        
        
(**********************************************************************)        
procedure ppsurf_to_rppsurf(f : ppsurf; 
                      var  rf : rppsurf);                  
var                                     
   i,j : integer;                       
begin (* procedure ppsurf_to_rppsurf *) 
  rf.degree1 := f.degree1;              
  rf.degree2 := f.degree2;              
  rf.num_segs1 := f.num_segs1;          
  rf.num_segs2 := f.num_segs2;          
  rf.dim := f.dim - 1;                  
  rf.break_point1 := f.break_point1;    
  rf.break_point2 := f.break_point2;    
                                        
   for i := 1 to f.num_segs1 do         
   for j := 1 to f.num_segs2 do         
   psurf_to_rpsurf(f.ppatch(.i,j.),rf.rppatch(.i,j.));     
                                        
end;  (* procedure ppsurf_to_rppsurf *) 
(**********************************************************************)        
                                   
(**********************************************************************)        
        
(*      name:  rppsurf_to_ppsurf                                      *)        
        
(*  category:  conversion between like entities.                      *)        
        
(*   purpose:  given a vector-valued rational piecewise polynomial    *)        
(*             surface, create the corresponding vector-valued        *)        
(*             piecewise polynomial surface by setting the last       *)        
(*             component of the output coefficients equal to the      *)        
(*             coefficient of the input denominator                   *)        
        
(*     input:  rf - a vector-valued rational piecewise polynomial     *)        
(*                  surface.                                          *)        
        
(*    output:   f - the corresponding vector valued (one more         *)        
(*                  dimension) piecewise polynomial                   *)        
(*                  surface.                                          *)        
        
(**********************************************************************)        
procedure rppsurf_to_ppsurf(rf : rppsurf;                  
                       var   f : ppsurf);                  
var                                     
   i,j : integer;                       
begin (* procedure rppsurf_to_ppsurf *) 
  f.degree1 := rf.degree1;              
  f.degree2 := rf.degree2;              
  f.num_segs1 := rf.num_segs1;          
  f.num_segs2 := rf.num_segs2;          
  f.dim := rf.dim + 1;                  
  f.break_point1 := rf.break_point1;    
  f.break_point2 := rf.break_point2;    
                                        
   for i := 1 to f.num_segs1 do         
   for j := 1 to f.num_segs2 do         
   rpsurf_to_psurf(rf.rppatch(.i,j.),f.ppatch(.i,j.));     
end;  (* procedure rppsurf_to_ppsurf *) 
(**********************************************************************)        
                                   
(**********************************************************************)        
        
(*      name:  bsurf_to_rbsurf                                        *)        
        
(*  category:  conversion between like entities.                      *)        
        
(*   purpose:  given a vector-valued bspline surface,                 *)        
(*             create the corresponding vector-valued rational        *)        
(*             bspline surface by interpreting the input b-spline     *)        
(*             surface as being in homogeneous form.                  *)        
        
(*     input:  g - a bspline surface.                                 *)        
        
(*    output:  rg - the corresponding vector valued (one less         *)        
(*                  dimension) rational bspline surface.              *)        
        
(**********************************************************************)        
procedure bsurf_to_rbsurf(g : bsurf;    
                     var rg : rbsurf);  
var                                     
   i,j : integer;                       
   temp : vector;                       
begin (* procedure bsurf_to_rbsurf *)   
   rg.degree1 := g.degree1;             
   rg.degree2 := g.degree2;             
   rg.dim := g.dim - 1;                 
   rg.num_of_knots1 := g.num_of_knots1; 
   rg.num_of_knots2 := g.num_of_knots2; 
   rg.knot1 := g.knot1;                 
   rg.knot2 := g.knot2;                 
                                        
   for i := 1 to g.num_of_knots1 - g.degree1 - 1 do        
   for j := 1 to g.num_of_knots2 - g.degree2 - 1 do        
   begin                                
      rg.weight(.i,j.) := g.coeff(.i,j.).component(.g.dim.); 
      temp := project(g.coeff(.i,j.));  
      rg.coeff(.i,j.) := quotient(temp,rg.weight(.i,j.));  
   end;                                 
end;  (* procedure bsurf_to_rbsurf *)   
(**********************************************************************)        
                                   
(**********************************************************************)        
        
(*      name:  rbsurf_to_bsurf                                        *)        
        
(*  category:  conversion between like entities.                      *)        
        
(*   purpose:  given a vector-valued rational b-spline                *)        
(*             surface, create the corresponding vector-valued        *)        
(*             b-spline surface by expressing the input               *)        
(*             surface in homogeneous form.                           *)        
        
(*     input:  rg - a vector-valued rational b-spline surface         *)        
        
(*    output:   g - the corresponding vector valued (one more         *)        
(*                  dimension) parametric b-spline                    *)        
(*                  surface.                                          *)        
        
(**********************************************************************)        
procedure rbsurf_to_bsurf(rg : rbsurf;  
                       var g : bsurf);  
var                                     
   i,j : integer;                       
   temp : vector;                       
begin (* procedure rbsurf_to_bsurf *)   
   g.degree1 := rg.degree1;             
   g.degree2 := rg.degree2;             
   g.dim := rg.dim + 1;                 
   g.num_of_knots1 := rg.num_of_knots1; 
   g.num_of_knots2 := rg.num_of_knots2; 
   g.knot1 := rg.knot1;                 
   g.knot2 := rg.knot2;                 
                                        
   for i := 1 to g.num_of_knots1 - g.degree1- 1 do         
   for j := 1 to g.num_of_knots2 - g.degree2- 1 do         
   begin                                
      temp := product(rg.weight(.i,j.),rg.coeff(.i,j.));   
      g.coeff(.i,j.) := append(temp,rg.weight(.i,j.));     
   end;                                 
end;  (* procedure rbsurf_to_bsurf *)   
(**********************************************************************)        
                                   
(**********************************************************************)        
        
(*      name:  deboor_fix                                             *)        
        
(*  category:  conversion between unlike entities                     *)        
        
(*   purpose:  given a piecewise polynomial and a knot sequence,      *)        
(*             find the corresponding b-spline coefficients.          *)        
(*             this procedure applies a result of deboor and          *)        
(*             fix described on page 113 of deboor's book,            *)        
(*             a practical guide to splines , springer-verlag,        *)        
(*             berlin-heidelberg-new york 1978.                       *)        
        
(*     input:  f - a piecewise polynomial                             *)        
(*             number_of_knots - the number of knots and              *)        
(*                           t - the knot sequence                    *)        
        
(*    output:  alpha - the corresponding b-spline coefficients        *)        
(*                     for the given function with the given          *)        
(*                     knot sequence                                  *)        
        
(**********************************************************************)        
                                        
procedure deboor_fix(f : pp;            
        number_of_knots : integer;      
                      t : real44;       
             var alpha  : vec44);       
var                                     
    psi : realmd;                       
    i,j,k,r : integer;                  
    top : real;                         
    del : realmd;                       
    s : realmd;                         
    vec1 : vector;                      
    ppvalue : vecmd;                    
    tau : real;                         
    codegree : integer;                 
    order : integer;                    
                                        
begin (* procedure deboor_fix *)        
                                        
order := f.degree + 1;                  
                                        
for i := 1 to number_of_knots - order do
begin (* finding alpha(.i.)*)           
   tau := 0.5*( t(.i.) + t(.i + 1.) );  
   for j := 1 to f.degree do del(.j.) := t(.i+j.) - tau;   
                                        
   symmetric(f.degree,del,s);           
                                        
   for k := 0 to f.degree do            
   begin (* k loop *)                   
      top := power(-1,k)*factorial(k)*s(.f.degree - k.);   
      psi(.k.) := top/factorial(f.degree);                 
   end;  (* k loop *)                   
                                        
   alpha(.i.) := zerovec(f.dim);        
   eval_pp(f,tau,f.degree,ppvalue);     
   for r := 0 to f.degree do            
   begin (* r loop *)                   
      codegree := f.degree - r;         
      (*=========================================================               
      alpha(.i.) := alpha(.i.) +        
      power(-1,codegree)*psi(.codegree.)*ppvalue(.r.);     
      ==========================================================*)              
      vec1 := product(power(-1,codegree)*psi(.codegree.),ppvalue(.r.));         
      alpha(.i.) := sum(alpha(.i.),vec1);                  
                                        
   end;  (* r loop *)                   
end; (* finding alpha(.i.) *)           
/*                                      
   if idebug = 1 then                   
writeln(debug,' leaving procedure deboor_fix');            
*/                                      
end; (* procedure deboor_fix *)         
                                        
(**********************************************************************)        
                                   
(**********************************************************************)        
        
(*      name:  deboor_fix2                                            *)        
        
(*  category:  conversion between unlike entities                     *)        
        
(*   purpose:  given a piecewise polynomial surface,                  *)        
(*             and given two knot sequences,                          *)        
(*             find the corresponding b-spline coefficients.          *)        
(*             this is the tensor product version of a result of      *)        
(*             deboor and fix described on page 113 of deboor's book, *)        
(*             a practical guide to splines , springer-verlag,        *)        
(*             berlin-heidelberg-new york 1978.                       *)        
        
(*     input:  f - a piecewise polynomial surface                     *)        
(*            n1 - number of knots in 1st parameter                   *)        
(*            n2 - number of knots in 2nd parameter                   *)        
(*            t1 - 1st knot sequence                                  *)        
(*            t2 - 2nd knot sequence                                  *)        
        
(*    output:  alpha - the corresponding b-spline coefficients        *)        
(*                     for the given function with the given          *)        
(*                     knot sequences                                 *)        
        
(**********************************************************************)        
procedure deboor_fix2(f : ppsurf;       
                     n1 : integer;      
                     n2 : integer;      
                     t1 : real44;       
                     t2 : real44;       
              var alpha : vec2222);     
var                                     
   psi1,psi2 : realmd;                  
   i,j,k,l,m,n,r,s : integer;           
   tau1,tau2,top : real;                
   del1,del2 : realmd;                  
   s1,s2 : realmd;                      
   vec1 : vector;                       
   svalue : vecmdmd;                    
   order1,order2 : integer;             
   codegree1,codegree2 : integer;       
   c1,c2 : real;                        
begin (* procedure deboor_fix2 *)       
                                        
   order1 := f.degree1 + 1;             
   order2 := f.degree2 + 1;             
                                        
for i := 1 to n1 - order1 do            
for j := 1 to n2 - order2 do            
begin (* finding alpha(i,j) *)          
                                        
   tau1 := 0.5*(t1(.i.) + t1(.i + 1.)); 
   tau2 := 0.5*(t2(.j.) + t2(.j + 1.)); 
                                        
   for l := 1 to f.degree1 do del1(.l.) := t1(.i+l.) - tau1; 
   for l := 1 to f.degree2 do del2(.l.) := t2(.j+l.) - tau2; 
                                        
   symmetric(f.degree1,del1,s1);        
   symmetric(f.degree2,del2,s2);        
                                        
   for m := 0 to f.degree1 do           
   begin (* m loop *)                   
      top := power(-1,m)*factorial(m)*s1(.f.degree1 - m.); 
      psi1(.m.) := top/factorial(f.degree1);               
   end;  (* m loop *)                   
                                        
   for n := 0 to f.degree2 do           
   begin (* n loop *)                   
      top := power(-1,n)*factorial(n)*s2(.f.degree2 - n.); 
      psi2(.n.) := top/factorial(f.degree2);               
   end;  (* n loop *)                   
                                        
   alpha(.i,j.) := zerovec(f.dim);      
                                        
   eval_ppsurf(f,f.degree1,f.degree2,tau1,tau2,svalue);    
                                        
   for r := 0 to f.degree1 do           
   for s := 0 to f.degree2 do           
   begin (* (r,s)-loop *)               
                                        
      codegree1 := f.degree1 - r;       
      codegree2 := f.degree2 - s;       
                                        
      c1 := power(-1,codegree1)*psi1(.codegree1.);         
      c2 := power(-1,codegree2)*psi2(.codegree2.);         
                                        
      vec1 := product(c1*c2,svalue(.r,s.));                
                                        
      alpha(.i,j.) := sum(alpha(.i,j.),vec1);              
                                        
   end;  (* (r,s)-loop *)               
end;  (* finding alpha(i,j) *)          
end;  (* procedure deboor_fix2 *)       
(**********************************************************************)        
                                   
(**********************************************************************)        
        
(*      name:  pp_to_bspline                                          *)        
        
(*  category:  conversion between unlike entities                     *)        
        
(*   purpose:  converts a piecewise polynomial to a bspline by        *)        
(*             (1) -  checking derivative continuity to determine     *)        
(*                    the knot sequence                               *)        
(*             (2) -  applying the deboor-fix algorithm to obtain     *)        
(*                    the coefficients                                *)        
        
(*     input:  f - a piecewise polynomial                             *)        
        
(*    output:  g - the corresponding bspline                          *)        
        
(**********************************************************************)        
                                        
procedure pp_to_bspline(f : pp;         
                    var g : bspline);   
var                                     
   number_of_knots : integer;           
   t : real44;                          
   alpha : vec44;                       
   order : integer;                     
   tol : real;                          
begin (* procedure pp_to_bspline *)     
                                        
   order := f.degree + 1;               
   tol := exp(-30 + f.degree);          
   find_knots(f,order,tol,number_of_knots,t);              
                                        
   deboor_fix(f,number_of_knots,t,alpha);                  
                                        
      g.degree := f.degree;             
      g.dim := f.dim;                   
      g.number_of_knots := number_of_knots;                
      g.knot := t;                      
      g.coeff := alpha;                 
                                        
end;  (* procedure pp_to_bspline *)     
                                        
(**********************************************************************)        
                                   
(**********************************************************************)        
        
(*      name:  bspline_to_pp                                          *)        
        
(*  category:  conversion between unlike entities                     *)        
        
(*   purpose:  converts a bspline to a piecewise polynomial by        *)        
(*             (1)  finding the distinct knots of the bspline         *)        
(*                  to find the breakpoints of the piecewise poly.    *)        
(*             (2)  taking derivatives at the breakpoints to find     *)        
(*                  the coefficients of each polynomial.              *)        
        
(*     input:  g - a bspline curve                                    *)        
        
(*    output:  f - the corresponding piecewise polynomial             *)        
        
(**********************************************************************)        
                                        
procedure bspline_to_pp(g : bspline;    
                   var  f : pp);        
var                                     
   i,j : integer;                       
   n_distinct_knots : integer;          
   distinct_knot : real44;              
   e_lee : vecmd;                       
                                        
begin (* procedure bspline_to_pp *)     
                                        
   find_distinct_knots(g.number_of_knots,                  
                       g.knot,          
                       n_distinct_knots,
                       distinct_knot);  
                                        
   f.num_segs := n_distinct_knots - 1;  
   f.dim := g.dim;                      
   f.degree := g.degree;                
                                        
   for i := 1 to n_distinct_knots do    
   f.break_point(.i.) := distinct_knot(.i.);               
                                        
   with f do                            
   for i := 1 to num_segs do            
   begin (* defining i-th segment of f *)                  
      eval_bspline(g,distinct_knot(.i.),g.degree,e_lee);   
                                        
      segment(.i.).degree := g.degree;  
      segment(.i.).dim := g.dim;        
                                        
      for j:=0 to segment(.i.).degree do
      segment(.i.).coeff(.j.) := quotient(e_lee(.j.),factorial(j));             
                                        
   end;   (* defining i-th segment of f *)                 
                                        
end; (*procedure b_spline_to_pp *)      
                                        
(**********************************************************************)        
                                   
(**********************************************************************)        
        
(*      name:  rpp_to_rbspline                                        *)        
        
(*  category:  conversion between unlike entities                     *)        
        
(*   purpose:  given a rational piecewise polynomial curve            *)        
(*             express it in rational b-spline form                   *)        
        
(*     input:  rf - a rational piecewise polynomial curve             *)        
        
(*    output:  rg - the input curve expressed in rational b-spline    *)        
(*                  form                                              *)        
        
(**********************************************************************)        
procedure rpp_to_rbspline(rf : rational_pp;                
                      var rg : rational_bspline);          
var                                     
   f : pp;                              
   g : bspline;                         
begin (* procedure rpp_to_rbspline *)   
                                        
(* 1 *) rpp_to_pp(rf,f);                
                                        
(* 2 *) pp_to_bspline(f,g);             
                                        
(* 3 *) bspline_to_rbspline(g,rg);      
                                        
(* 4 *) rbs_fraction(rg);               
                                        
end;  (* procedure rpp_to_bspline *)    
(**********************************************************************)        
                                   
(**********************************************************************)        
        
(*      name:  rbspline_to_rpp                                        *)        
        
(*  category:  conversion between unlike entities                     *)        
        
(*   purpose:  given a rational b-spline curve, express it in         *)        
(*             rational piecewise polynomial form                     *)        
        
(*     input:  rg - a rational b-spline curve                         *)        
        
(*    output:  rf - the input curve, expressed in rational b-spline   *)        
(*                  form.                                             *)        
        
(**********************************************************************)        
procedure rbspline_to_rpp(rg : rational_bspline;           
                      var rf : rational_pp);               
var                                     
   i : integer;                         
   f : pp;                              
   g : bspline;                         
begin (* procedure rbspline_to_rpp *)   
                                        
(* 1 *) rbspline_to_bspline(rg,g);      
                                        
(* 2 *) bspline_to_pp(g,f);             
                                        
(* 3 *) pp_to_rpp(f,rf);                
                                        
(* 4 *) for i := 1 to rf.num_segs do    
        rp_fraction(rf.rpp_segment(.i.));                  
                                        
end;  (* procedure rbspline_to_rpp *)   
(**********************************************************************)        
                                   
(**********************************************************************)        
        
(*      name:  ppsurf_to_bsurf                                        *)        
        
(*  category:  conversion between unlike entities                     *)        
        
(*   purpose:  given a piecewise polynomial surface, find the         *)        
(*             equivalent b-spline surface.  the method used is:      *)        
        
(*             (1) - invoke procedure find_knots2 to find the         *)        
(*                   appropriate knot sequence in each parameter      *)        
(*                   that would be needed to represent the surface    *)        
(*                   as a tensor product patchwise bezier surface     *)        
(*                   in b-spline form.                                *)        
        
(*             (2) - invoke procedure deboor_fix2 to use the          *)        
(*                   tensor product version of the deboor-fix         *)        
(*                   algorithm to actually create the b-spline        *)        
(*                   surface corresponding to the input function      *)        
(*                   and to the knot sequence found in step (1).      *)        
        
(*             (3) - invoke remove_all_surface knots to the tensor    *)        
(*                   product patchwise bezier surface created in step *)        
(*                   (2) to obtain an equivalent b-spline surface     *)        
(*                   with a minimal number of knots.                  *)        
        
        
        
(*     input:  p - a piecewise polynomial surface                     *)        
        
        
(*    output:  b - the equivalent b-spline surface                    *)        
        
        
(**********************************************************************)        
                                        
procedure ppsurf_to_bsurf(p : ppsurf;   
                          var b : bsurf);                  
var                                     
   temp : bsurf;                        
   tol : real;                          
                                        
begin (* procedure ppsurf_to_rbsurf *)  
                                        
b.degree1 := p.degree1;                 
b.degree2 := p.degree2;                 
b.dim := p.dim;                         
                                        
find_knots2(p,b.num_of_knots1,b.num_of_knots2,b.knot1,b.knot2);                 
                                        
deboor_fix2(p,b.num_of_knots1,b.num_of_knots2,b.knot1,b.knot2,b.coeff);         
                                        
temp := b;                              
tol := exp(-20);                        
remove_all_surf_knots(temp,tol,b);      
                                        
end;  (* procedure ppsurf_to_rbsurf *)  
(**********************************************************************)        
                                   
(**********************************************************************)        
        
(*      name:  rppsurf_to_rbsurf                                      *)        
        
(*  category:  conversion between unlike entities                     *)        
        
(*   purpose:  convert a rational piecewise polynomial surface        *)        
(*             to a rational b-spline surface.                        *)        
        
(*     input:  rp - a rational piecewise polynomial surface           *)        
        
(*    output:  rb - the same surface expressed in a tensor product    *)        
(*                  rational b-spline form.                           *)        
        
(**********************************************************************)        
procedure rppsurf_to_rbsurf(rp : rppsurf;                  
                            var rb : rbsurf);              
var                                     
   p : ppsurf;                          
   b : bsurf;                           
                                        
begin (* procedure rppsurf_to_rbsurf *) 
                                        
(* 1 *) rppsurf_to_ppsurf(rp,p);        
                                        
(* 2 *) ppsurf_to_bsurf(p,b);           
                                        
(* 3 *) bsurf_to_rbsurf(b,rb);          
                                        
(* 4 *) rbsurf_fraction(rb);            
                                        
end;  (* procedure rppsurf_to_rbsurf *) 
(**********************************************************************)        
                                   
(**********************************************************************)        
        
(*      name:  bsurf_to_ppsurf                                        *)        
        
(*  category:  conversion between unlike entities                     *)        
        
(*   purpose:  converts a bspline surface to a piecewise polynomial   *)        
(*             surface by:*)            
(*             (1)  finding the distinct knots of the bspline         *)        
(*                  to find the breakpoints of the piecewise poly.    *)        
(*             (2)  taking derivatives at the breakpoints to find     *)        
(*                  the coefficients of each polynomial patch.        *)        
        
(*     input:  b - a bspline surface                                  *)        
        
(*    output:  p - the corresponding piecewise polynomial surface     *)        
        
(**********************************************************************)        
procedure bsurf_to_ppsurf(b : bsurf;    
                      var p : ppsurf);  
var                                     
   i,j,k,m,n : integer;                 
   n_distinct_knots : integer;          
   distinct_knot : real44;              
   svalue : vecmdmd;                    
begin (* procedure bsurf_to_ppsurf *)   
   find_distinct_knots(b.num_of_knots1, 
                       b.knot1,         
                       n_distinct_knots,
                       distinct_knot);  
   p.num_segs1 := n_distinct_knots - 1; 
   for i := 1 to n_distinct_knots do    
   p.break_point1(.i.) := distinct_knot(.i.);              
                                        
   find_distinct_knots(b.num_of_knots2, 
                       b.knot2,         
                       n_distinct_knots,
                       distinct_knot);  
   p.num_segs2 := n_distinct_knots - 1; 
   for i := 1 to n_distinct_knots do    
   p.break_point2(.i.) := distinct_knot(.i.);              
                                        
   p.degree1 := b.degree1;              
   p.degree2 := b.degree2;              
   p.dim := b.dim;                      
                                        
   with p do                            
   for i := 1 to num_segs1 do           
   for j := 1 to num_segs2 do           
   begin (* defining (i,j)-th patch *)  
      ppatch(.i,j.).degree1 := b.degree1;                  
      ppatch(.i,j.).degree2 := b.degree2;                  
      ppatch(.i,j.).dim := b.dim;       
                                        
      eval_bsurf(b,                     
                 p.degree1,             
                 p.degree2,             
                 p.break_point1(.i.),   
                 p.break_point2(.j.),   
                 svalue);               
                                        
      for m := 0 to p.degree1 do        
      for n := 0 to p.degree2 do        
      ppatch(.i,j.).coeff(.m,n.) :=     
      quotient(svalue(.m,n.),factorial(m)*factorial(n));   
                                        
   end;  (* defining (i,j)-th patch *)  
end;  (* procedure bsurf_to_ppsurf *)   
(**********************************************************************)        
                                   
(**********************************************************************)        
        
(*      name:  rbsurf_to_rppsurf                                      *)        
        
(*  category:  conversion between unlike entities                     *)        
        
(*   purpose:  convert from a rational b-spline surface to the        *)        
(*             corresponding rational piecewise polynomial            *)        
(*             surface                                                *)        
        
(*     input:  rb - a rational b-spline surface                       *)        
        
(*    output:  rp - the corresponding rational piecewise polynomial   *)        
(*                  surface                                           *)        
        
(**********************************************************************)        
procedure rbsurf_to_rppsurf(rb : rbsurf; var rp : rppsurf);
var                                     
   b : bsurf;                           
   p : ppsurf;                          
   i,j : integer;                       
                                        
begin (* procedure rbsurf_to_rppsurf *) 
                                        
(* 1 *) rbsurf_to_bsurf(rb,b);          
                                        
(* 2 *) bsurf_to_ppsurf(b,p);           
                                        
(* 3 *) ppsurf_to_rppsurf(p,rp);        
                                        
(* 4 *) for i := 1 to rp.num_segs1 do   
        for j := 1 to rp.num_segs2 do   
        rpsurf_fraction(rp.rppatch(.i,j.));                
                                        
end;  (* procedure rbsurf_to_rppsurf *) 
(**********************************************************************)        
                                   
(**********************************************************************)        
        
(*      name:  c_eval_and_compare                                     *)        
        
(*  category:  test utilities                                         *)        
        
(*   purpose:  compare a given rational piecewise polynomial curve    *)        
(*             and a rational b-spline curve by evaluating both       *)        
(*             of them at a designated number of points               *)        
        
(*     input:  f - a rational piecewise polynomial curve              *)        
(*             g - a rational b-spline curve                          *)        
(*             n - the curves are evaluated at each 1/n-th point      *)        
(*                 of the domain                                      *)        
        
(*    output:  no output parameters. for n + 1 evenly spaced values   *)        
(*             of x the vector valued quantities f(x),g(x) and        *)        
(*             f(x)-g(x) are written to the report file.  the maximum *)        
(*             and minimum norms, and the maximum value of            *)        
(*             !!f(x)-g(x)!! are written to the terminal and to the   *)        
(*             report file                                            *)        
        
(**********************************************************************)        
procedure c_eval_and_compare(f : rational_pp;              
                             g : rational_bspline;         
                             n : integer);                 
                                        
var                                     
   i,j,k : integer;                     
   u : real;                            
   max_error : real;                    
   rppnorm : real;                      
   rbsnorm : real;                      
   max_rppnorm : real;                  
   min_rppnorm : real;                  
   max_rbsnorm : real;                  
   min_rbsnorm : real;                  
   rpp_value,rbs_value : vector;        
   diff : vector;                       
   start_pt,end_pt : real;              
   del :real;                           
                                        
begin (* procedure c_eval_and_compare *)
   max_error := 0.0;                    
   max_rppnorm := 0.0;                  
   min_rppnorm := 999999999999.0;       
   max_rbsnorm := 0.0;                  
   min_rbsnorm := 999999999999.0;       
                                        
   del := 1.0/n;                        
   start_pt := f.break_point(.1.);      
   end_pt   := f.break_point(.f.num_segs + 1.);            
   for j:= 0 to n do                    
   begin (* function evaluation at u *) 
      u := (n-j)*(del)*start_pt + j*(del)*end_pt;          
      if j=0  then u:=start_pt;         
      if j=n then u:=end_pt;            
      rpp_value := eval_rpp(f,u);       
      rbs_value := eval_rbs(g,u);       
      writeln(report,' ');              
      diff := difference(rpp_value,rbs_value);             
      rppnorm := norm(rpp_value);       
      rbsnorm := norm(rbs_value);       
      max_error := max(norm(diff),max_error);              
      max_rppnorm := max(rppnorm,max_rppnorm);             
      min_rppnorm := min(rppnorm,min_rppnorm);             
      max_rbsnorm := max(rbsnorm,max_rbsnorm);             
      min_rbsnorm := min(rbsnorm,min_rbsnorm);             
                                        
      write(report,u:4:2,'  rpp = ');   
      writevec(report,rpp_value);       
      writeln(report);                  
      write(report,'      rbs = ');     
      writevec(report,rbs_value);       
      writeln(report);                  
      write(report,'     diff = ');     
      writevec(report,diff);            
      writeln(report);                  
                                        
   end;  (* function evaluation at u *) 
                                        
   cms('clrscrn',return_code);          
   writeln(report);                     
   writeln(report,'............................................');              
   writeln(report);                     
   writeln(output,' maximum error in difference ',         
   ' = ',max_error);                    
                                        
   writeln(report,' maximum error in difference ',         
   ' = ',max_error);                    
                                        
   writeln(output,' max_rppnorm = ',max_rppnorm:20:16);    
   writeln(output,' max_rbsnorm = ',max_rbsnorm:20:16);    
   writeln(output,' min_rppnorm = ',min_rppnorm:20:16);    
   writeln(output,' min_rbsnorm = ',min_rbsnorm:20:16);    
   writeln(report,' max_rppnorm = ',max_rppnorm:20:16);    
   writeln(report,' max_rbsnorm = ',max_rbsnorm:20:16);    
   writeln(report,' min_rppnorm = ',min_rppnorm:20:16);    
   writeln(report,' min_rbsnorm = ',min_rbsnorm:20:16);    
                                        
   writeln(report);                     
   writeln(report,'............................................');              
   writeln(report);                     
end;  (* procedure c_eval_and_compare *)
(**********************************************************************)        
                                   
(**********************************************************************)        
        
(*      name:  rpp_compare                                            *)        
        
(*  category:  test utilities                                         *)        
        
(*   purpose:  find the maximum difference in coefficients between    *)        
(*             two rational piecewise polynomial vector-valued        *)        
(*             parametric curves.                                     *)        
        
(*     input:  f1 - a rational piecewise polynomial curve             *)        
(*             f2 - a rational piecewise polynomial curve             *)        
        
(*    output:  no output parameters.  if f1 and f2 are "comparable"   *)        
(*             (see code below), the maximum distance between         *)        
(*             the coefficents of f1 and the coefficients of f2       *)        
(*             are written.                                           *)        
        
(**********************************************************************)        
                                        
procedure rpp_compare(f1 : rational_pp; f2 : rational_pp); 
var                                     
   i,j : integer;                       
   top_coe_diff,bot_coe_diff : real;    
begin (* procedure rpp_compare *)       
   top_coe_diff := 0.0;               
   bot_coe_diff := 0.0;               
                 
   if (f1.num_segs <> f2.num_segs) or 
      (f1.degree   <> f2.degree  ) or 
      (f1.dim      <> f2.dim     )    
   then          
   begin         
      writeln(output,' in rpp_compare, functions not comparable');              
      writeln(report,' in rpp_compare, functions not comparable');              
   end           
   else          
begin (*normal case *)                
   for i := 1 to f1.num_segs do       
   for j := 0 to f1.rpp_segment(.i.).degree do             
   begin (* finding maximum errors *) 

      top_coe_diff :=                 
      max(top_coe_diff,dist(f1.rpp_segment(.i.).top_coeff(.j.),                 
                            f2.rpp_segment(.i.).top_coeff(.j.)));               
      bot_coe_diff :=                 
      max(bot_coe_diff,abs(f1.rpp_segment(.i.).bot_coeff(.j.) -                 
                           f2.rpp_segment(.i.).bot_coeff(.j.)));                
   end;  (* finding maximum errors *) 
                 
writeln(output,'maximum difference between top coeffs = ',top_coe_diff);        
writeln(report,'maximum difference between top coeffs = ',top_coe_diff);        
                 
writeln(output,'maximum difference between bot coeffs = ',bot_coe_diff);        
writeln(report,'maximum difference between bot coeffs = ',bot_coe_diff);        
                 
   end;  (* normal case *)            
                 
end;  (* procedure rpp_compare *)     
(**********************************************************************)        
            
(**********************************************************************)        
        
(*      name:  rbs_compare                                            *)        
        
(*  category:  test utilities                                         *)        
        
(*   purpose:  find the maximum difference in coefficients between    *)        
(*             two rational b-spline vector-valued                    *)        
(*             parametric curves.                                     *)        
        
(*     input:  f1 - a rational b-spline curve                         *)        
(*             f2 - a rational b-spline curve                         *)        
        
(*    output:  no output parameters.  if f1 and f2 are "comparable"   *)        
(*             (see code below), the maximum distance between         *)        
(*             the coefficents of f1 and the coefficients of f2       *)        
(*             are written.  the same analysis is done with the       *)        
(*             weights.                                               *)        
        
(**********************************************************************)        
procedure rbs_compare(f1 : rational_bspline; f2 : rational_bspline);            
var              
   j,n1,n2: integer;                  
   coe_diff,weight_diff : real;       
begin (* procedure rbs_compare *)     
   coe_diff := 0.0;
   weight_diff := 0.0;                
   n1 := f1.number_of_knots - f1.degree - 1;               
   n2 := f2.number_of_knots - f2.degree - 1;               
   if (n1 <> n2) or
   (f1.degree <> f2.degree) or        
   (f1.dim <> f2.dim)                 
   then          
   begin         
      writeln(output,' in rbs_compare, splines not comparable');                
      writeln(report,' in rbs_compare, splines not comparable');                
   end           
   else          
begin (* normal case *)               
   for j := 1 to n1 do                
   begin (* finding maximum errors *) 
      coe_diff :=
      max(coe_diff,dist(f1.coeff(.j.),
                        f2.coeff(.j.)));

      weight_diff :=
      max(weight_diff,abs(f1.weight(.j.) -
                          f2.weight(.j.)));
   end;  (* finding maximum errors *) 

   writeln(output,' maximum difference between coeffs = ',coe_diff);
   writeln(report,' maximum difference between coeffs = ',coe_diff);

   writeln(output,' maximum difference between weights = ',weight_diff);
   writeln(report,' maximum difference between weights = ',weight_diff);

   end;  (* normal case *)
end;  (* procedure rbs_compare *)

(**********************************************************************)

(**********************************************************************)

(*      name:  rppsurf_compare                                        *)

(*  category:  test utilities                                         *)

(*   purpose:  find the maximum difference in coefficients between    *)
(*             two rational piecewise polynomial vector-valued        *)
(*             parametric surfaces                                    *)

(*     input:  f1 - a rational piecewise polynomial surface           *)
(*             f2 - a rational piecewise polynomial surface           *)

(*    output:  no output parameters.  if f1 and f2 are "comparable"   *)
(*             (see code below), the maximum distance between         *)
(*             the coefficents of f1 and the coefficients of f2       *)
(*             are written.                                           *)

(**********************************************************************)

procedure rppsurf_compare(f1 : rppsurf; f2 : rppsurf);
var
   i,j,k,l : integer;
   top_coe_diff,bot_coe_diff : real;  
begin (* procedure rppsurf_compare *) 
   top_coe_diff := 0.0;
   bot_coe_diff := 0.0;

   if (f1.num_segs1 <> f2.num_segs1) or 
      (f1.num_segs2 <> f2.num_segs2) or 
      (f1.degree1 <> f2.degree1) or   
      (f1.degree2 <> f2.degree2) or   
      (f1.dim <> f2.dim)
   then
   begin
      writeln(output,' in rppsurf_compare, functions not comparable');
      writeln(report,' in rppsurf_compare, functions not comparable');
   end
   else
begin (*normal case *)
   for i := 1 to f1.num_segs1 do
   for j := 1 to f1.num_segs2 do
   for k := 0 to f1.rppatch(.i,j.).degree1 do
   for l := 0 to f1.rppatch(.i,j.).degree2 do
   begin (* finding maximum errors *) 
      top_coe_diff :=
      max(top_coe_diff,dist(f1.rppatch(.i,j.).top_coeff(.k,l.),
                            f2.rppatch(.i,j.).top_coeff(.k,l.)));
      bot_coe_diff :=
      max(bot_coe_diff,abs(f1.rppatch(.i,j.).bot_coeff(.k,l.) -
                           f2.rppatch(.i,j.).bot_coeff(.k,l.)));
   end;  (* finding maximum errors *) 

writeln(output,'maximum difference between top coeffs = ',top_coe_diff);
writeln(report,'maximum difference between top coeffs = ',top_coe_diff);

writeln(output,'maximum difference between bot coeffs = ',bot_coe_diff);
writeln(report,'maximum difference between bot coeffs = ',bot_coe_diff);

   end;  (* normal case *)

end;  (* procedure rppsurf_compare *) 
(**********************************************************************)

(**********************************************************************)

(*      name:  rbsurf_compare                                         *)

(*  category:  test utilities                                         *)

(*   purpose:  find the maximum difference in coefficients between    *)
(*             two rational b-spline vector-valued                    *)
(*             parametric surfaces                                    *)

(*     input:  f1 - a rational b-spline surface                       *)
(*             f2 - a rational b-spline surface                       *)

(*    output:  no output parameters.  if f1 and f2 are "comparable"   *)
(*             (see code below), the maximum distance between         *)
(*             the coefficents of f1 and the coefficients of f2       *)
(*             are written.  the same analysis is done with the       *)
(*             weights.                                               *)

(**********************************************************************)
procedure rbsurf_compare(f1 : rbsurf; f2 : rbsurf);
var
   i,j,n11,n12,n21,n22: integer;
   coe_diff,weight_diff : real;
begin (* procedure rbsurf_compare *)  
   coe_diff := 0.0;
   weight_diff := 0.0;
   n11 := f1.num_of_knots1 - f1.degree1 - 1;
   n12 := f1.num_of_knots2 - f1.degree2 - 1;
   n21 := f2.num_of_knots1 - f2.degree1 - 1;
   n22 := f2.num_of_knots2 - f2.degree2 - 1;
   if (n11 <> n21) or
      (n12 <> n22) or
      (f1.degree1 <> f2.degree1) or   
      (f1.degree2 <> f2.degree2) or   
      (f1.dim <> f2.dim)
   then
   begin
      writeln(output,' in rbsurf_compare, splines not comparable');
      writeln(report,' in rbsurf_compare, splines not comparable');
   end
   else
begin (* normal case *)
   for i := 1 to n11 do
   for j := 1 to n22 do
   begin (* finding maximum errors *) 
      coe_diff :=
      max(coe_diff,dist(f1.coeff(.i,j.),
                        f2.coeff(.i,j.)));
      weight_diff :=
      max(weight_diff,abs(f1.weight(.i,j.) -
                          f2.weight(.i,j.)));
   end;  (* finding maximum errors *) 

   writeln(output,' maximum difference between coeffs = ',coe_diff);
   writeln(report,' maximum difference between coeffs = ',coe_diff);

   writeln(output,' maximum difference between weights = ',weight_diff);
   writeln(report,' maximum difference between weights = ',weight_diff);

   end;  (* normal case *)
end;  (* procedure rbsurf_compare *)  

(**********************************************************************)

(**********************************************************************)

(*      name:  s_eval_and_compare                                     *)

(*  category:  test utilities                                         *)

(*   purpose:  compare a given rational piecewise polynomial surface  *)
(*             and a rational b-spline surface by evaluating both     *)
(*             of them at a designated number of points               *)

(*     input:  p - a rational piecewise polynomial surface            *)
(*             b - a rational b-spline surface                        *)
(*             m - the curves are evaluated on a grid by taking       *)
(*             n - parameter values at increments of 1/m and 1/n      *)
(*                 of the domain in each parameter                    *)

(*    output:  no output parameters. for evenly spaced values of x    *)
(*             and y, the vector valued quantities p(x,y),b(x,y) and  *)
(*             p(x,y)-b(x,y) are written to the report file.  the     *)
(*             maximum and minimum norms, and the maximum value of    *)
(*             !!p(x,y)-b(x,y)!! are written to the terminal and      *)
(*             to the report file                                     *)

(**********************************************************************)
procedure s_eval_and_compare(p : rppsurf;
                             b : rbsurf;
                             m : integer;
                             n : integer);
var
   i,j ,k : integer;
   max_error : real;
   rppnorm : real; 
   rbsnorm : real; 
   max_rppnorm : real;
   min_rppnorm : real;
   max_rbsnorm : real;
   min_rbsnorm : real;
   u,v : real;   
   del1,del2 : real;
   start1,start2,end1,end2 : real;    
   rbs_value,rpp_value : vector;
   diff : vector;
begin (* procedure s_eval_and_compare *)

max_error := 0.0;
max_rppnorm := 0.0;
min_rppnorm := 999999999999.0;
max_rbsnorm := 0.0;
min_rbsnorm := 999999999999.0;

      del1 := 1.0/m;
      del2 := 1.0/n;
      start1 := b.knot1(.1.);
      start2 := b.knot2(.1.);
      end1   := b.knot1(.b.num_of_knots1.);
      end2   := b.knot2(.b.num_of_knots2.);

      for i := 0 to m do
      for j := 0 to n do
      begin (* evaluating at (i,j)-th point *)
         writeln(i,j);
         u := (m-i)*(del1)*start1 + i*(del1)*end1;
         v := (n-j)*(del2)*start2 + j*(del2)*end2;
         if i=0 then u := start1;
         if i=m then u := end1;
         if j=0 then v := start2;
         if j=n then v := end2;
         rbs_value := eval_rbsurf(b,u,v);
         rpp_value := eval_rppsurf(p,u,v);

writeln(report,' ');
diff := difference(rpp_value,rbs_value);
rppnorm := norm(rpp_value);
rbsnorm := norm(rbs_value);
max_error := max(norm(diff),max_error); 
max_rppnorm := max(rppnorm,max_rppnorm);
min_rppnorm := min(rppnorm,min_rppnorm);
max_rbsnorm := max(rbsnorm,max_rbsnorm);
min_rbsnorm := min(rbsnorm,min_rbsnorm);

write(report,u:4:2,' ',v:4:2,'  rpp = ');
writevec(report,rpp_value);
writeln(report); 
write(report,'           rbs = ');    
writevec(report,rbs_value);
writeln(report); 
write(report,'          diff = ');    
writevec(report,diff);
writeln(report); 

      end;  (* evaluating at (i,j)-th point *)

cms('clrscrn',return_code);
writeln(report); 
writeln(report,'............................................');
writeln(report); 
writeln(output,' maximum error in difference ',
' = ',max_error);

writeln(report,' maximum error in difference ',
' = ',max_error);

writeln(output,' max_rppnorm = ',max_rppnorm:20:16);
writeln(output,' max_rbsnorm = ',max_rbsnorm:20:16);
writeln(output,' min_rppnorm = ',min_rppnorm:20:16);
writeln(output,' min_rbsnorm = ',min_rbsnorm:20:16);
writeln(report,' max_rppnorm = ',max_rppnorm:20:16);
writeln(report,' max_rbsnorm = ',max_rbsnorm:20:16);
writeln(report,' min_rppnorm = ',min_rppnorm:20:16);
writeln(report,' min_rbsnorm = ',min_rbsnorm:20:16);

writeln(report); 
writeln(report,'............................................');
writeln(report); 
end;  (* procedure s_eval_and_compare *)
(**********************************************************************)

(**********************************************************************)

(*      name:  start_with_rpp_curve                                   *)

(*  category:  test drivers                                           *)

(*   purpose:  this is a "driver" procedure which is invoked          *)
(*             if the user indicated that he wished to start with     *)
(*             a rational piecewise polynomial curve.                 *)
(*             the appropriate input, conversion, output              *)
(*             and testing routines are called.                       *)

(*     input:  no input parameters                                    *)

(*    output:  no output parameters                                   *)

(**********************************************************************)
procedure start_with_rpp_curve;
var
   i,k,n : integer;
begin (* procedure start_with_rpp_curve *)
repeat (* until degree <= 0 *)

read_rational_pp(f);
if f.degree > 0 then
begin (* converting the input function *)
writeln(report,' ****************************************************');
writeln(report,' input data for this case ');
write_rpp(report,f);

rpp_to_rbspline(f,g);

writeln(report,' ----------------------------------------------------');
writeln(report,' ');
writeln(report,'output data for this case:  ');

cms('clrscrn',return_code);
writeln(output,' ----------------------------------------------------');
writeln(output,' ');
writeln(output,'output data for this case:  ');
   write_rbs(output,g);
   write_rbs(report,g);

writeln(' ');    
writeln(' enter n >= 1 to evaluate at each 1/n point of the interval');
writeln(' enter 0 otherwise.');
readln(n);
if n >= 1 then c_eval_and_compare(f,g,n);

writeln(' ');    
writeln(' enter k > 0 if you want to test k repeated ');   
writeln(' conversions between pp and b-spline.  ');
writeln(' enter 0 otherwise.  ');
readln(k);
if k > 0 then    
begin (* iterative testing *)
f0 := f;
g0 := g;
for i := 1 to k do 
begin (* i-th pair of conversions *)  
   writeln(' starting conversion pair number ',i:1);
   rpp_to_rbspline(f,g);
   rbspline_to_rpp(g,f);
end;  (* i-th pair of conversions *)  
cms('clrscrn',return_code);
writeln(' ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ');
writeln(report,' +++++++++++++++++++++++++++++++++++++++++++++++++ ');
writeln(' ');    
writeln(report,' ');
writeln(' after ',i:2,' conversions, the result is:');
writeln(report,' after ',i:2,' conversions, the result is:');
   write_rpp(report,f);
   rpp_compare(f0,f);

end;  (* iterative testing *)

end;  (* converting the input function *)
until f.degree <= 0;
end;  (* procedure start_with_rpp_curve *)
(**********************************************************************)

(**********************************************************************)

(*      name:  start_with_rbs_curve                                   *)

(*  category:  test drivers                                           *)

(*   purpose:  this is a "driver" procedure which is invoked          *)
(*             if the user indicated that he wished to start with     *)
(*             a rational b-spline curve.                             *)
(*             the appropriate input, conversion, output              *)
(*             and testing routines are called.                       *)

(*     input:  no input parameters                                    *)

(*    output:  no output parameters                                   *)

(**********************************************************************)
procedure start_with_rbs_curve;
var
   x : real;
   i,j,k,n : integer;
   gtemp,gnew:bspline;
   removed : boolean;
   tol : real;   
begin (* procedure start_with_rbs_curve *)
repeat (* until degree <= 0 *)

read_rational_bspline(g);
if g.degree > 0 then
begin (* converting the input function *)
writeln(' enter number of knots to add to this b-spline curve ');
readln(n_knots_to_add);
for j := 1 to n_knots_to_add do
begin
writeln(' enter knot to be added:  ');
readln(x);
rbspline_to_bspline(g,gtemp);
add_knot(gtemp,x,gnew);
bspline_to_rbspline(gnew,g)
end;
if g.number_of_knots > 2*(g.degree + 1) then
begin (* knot removal section *)
writeln
(' enter number of knots to try to remove from this b-spline curve ');
readln(n_knots_to_remove);
for j := 1 to n_knots_to_remove do    
begin
   writeln(' enter index of knot to be removed:  ');
   writeln(' index must be > ',g.degree+1:1,' and < ',
   g.number_of_knots - g.degree:1);   
   readln(i);    
   tol := exp(-20);
   rbspline_to_bspline(g,gtemp);
   remove_knot(gtemp,i,tol,removed,gnew);
   writeln(' removed =',removed);
   if removed then bspline_to_rbspline(gnew,g)
end;
end; (* knot removal section *)
writeln(report,' ****************************************************');
writeln(report,' input data for this case ');
write_rbs(report,g);

rbspline_to_rpp(g,f);

writeln(report,' ----------------------------------------------------');
writeln(report,' ');
writeln(report,'output data for this case:  ');

cms('clrscrn',return_code);
writeln(output,' ----------------------------------------------------');
writeln(output,' ');
writeln(output,'output data for this case:  ');
   write_rpp(output,f);
   write_rpp(report,f);

writeln(' ');    
writeln(' enter n >= 1 to evaluate at each 1/n point of the interval');
writeln(' enter 0 otherwise.');
readln(n);
if n >= 1 then c_eval_and_compare(f,g,n);

writeln(' ');    
writeln(' enter k > 0 if you want to test k repeated ');   
writeln(' conversions between pp and b-spline.  ');
writeln(' enter 0 otherwise.  ');
readln(k);
if k > 0 then    
begin (* iterative testing *)
f0 := f;
g0 := g;
for i := 1 to k do 
begin (* i-th pair of conversions *)  
   writeln(' starting conversion pair number ',i:1);
   rbspline_to_rpp(g,f);
   rpp_to_rbspline(f,g);
end;  (* i-th pair of conversions *)  
cms('clrscrn',return_code);
writeln(' ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ');
writeln(report,' +++++++++++++++++++++++++++++++++++++++++++++++++ ');
writeln(' ');    
writeln(report,' ');
writeln(' after ',i:2,' conversions, the result is:');
writeln(report,' after ',i:2,' conversions, the result is:');
   write_rbs(report,g);
   rbs_compare(g0,g);
end;  (* iterative testing *)

end;  (* converting the input function *)
until g.degree <= 0;
end;  (* procedure start_with_rbs_curve *)
(**********************************************************************)

(**********************************************************************)

(*      name:  start_with_rpp_surface                                 *)

(*  category:  test drivers                                           *)

(*   purpose:  this is a "driver" procedure which is invoked          *)
(*             if the user indicated that he wished to start with     *)
(*             a rational piecewise polynomial surface,               *)
(*             the appropriate input, conversion, output              *)
(*             and testing routines are called.                       *)

(*     input:  no input parameters                                    *)

(*    output:  no output parameters                                   *)

(**********************************************************************)
procedure start_with_rpp_surface;
var
   p0 : rppsurf; 
   i,j,k,m,n : integer;
begin (* procedure start_with_rpp_surface *)
   repeat (* until p.degree1 <= 0 or p.degree2 <= 0 *)
      read_rppsurf(p);
      if (p.degree1 > 0) and (p.degree2 > 0) then
      begin (* conversion *)
         p0 := p;
         writeln(report,'********************************************');
         writeln(report,' input data for this case:');
         write_rppsurf(report,p);
         rppsurf_to_rbsurf(p,b);
         writeln(report,'--------------------------------------------');
         writeln(report);
         writeln(report,' output data for this case ');
         write_rbsurf(report,b);
         cms('clrscrn',return_code);
         writeln(output,'--------------------------------------------');
         writeln(output);
         writeln(output,' output data for this case ');
         write_rbsurf(output,b);
         writeln(' ');
         writeln(' enter m >= 1, n >= 1 to evaluate on an m x n grid ');
         read(m); read(n);
         if (m>=1) and (n>=1) then
         s_eval_and_compare(p,b,m,n);
      end;  (* conversion *)

writeln(' ');
writeln(' enter k > 0 if you want to test k repeated ');
writeln(' conversions between pp and b-spline.  ');
writeln(' enter 0 otherwise.  ');
readln(k);
if k > 0 then
begin (* iterative testing *)
for i := 1 to k do
begin (* i-th pair of conversions *)
   writeln(' starting conversion pair number ',i:1);
   rppsurf_to_rbsurf(p,b);
   rbsurf_to_rppsurf(b,p);
end;  (* i-th pair of conversions *)
cms('clrscrn',return_code);
writeln(report,' +++++++++++++++++++++++++++++++++++++++++++++++++ ');
writeln(report,' ');
writeln(report,' after ',i:2,' conversions, the result is:');
   write_rppsurf(report,p);
writeln(output,' +++++++++++++++++++++++++++++++++++++++++++++++++ ');
writeln(output,' ');
writeln(output,' after ',i:2,' conversions, the result is:');
   write_rppsurf(report,p);
   rppsurf_compare(p0,p);
end;  (* iterative testing *)

   until (p.degree1 <= 0) or (p.degree2 <= 0)
end;  (* procedure start_with_rpp_surface *)
(**********************************************************************)

(**********************************************************************)

(*      name:  start_with_rbs_surface                                 *)

(*  category:  test drivers                                           *)

(*   purpose:  this is a "driver" procedure which is invoked          *)
(*             if the user indicated that he wished to start with     *)
(*             a rational b-spline surface,                           *)
(*             the appropriate input, conversion, output              *)
(*             and testing routines are called.                       *)

(*     input:  no input parameters                                    *)

(*    output:  no output parameters                                   *)

(**********************************************************************)
procedure start_with_rbs_surface;
var
   x,tol : real;
   i,j,k,m,n,param : integer;
   b0 : rbsurf;
   btemp,bnew :bsurf;
   removed : boolean;
begin (* procedure start_with_rbs_surface *)
   repeat (* until b.degree1 <= 0 or b.degree2 <= 0 *)
   read_rbsurf(b);
   if (b.degree1 > 0) and (b.degree2 > 0) then
   begin (* conversion *)
   b0 := b;
writeln(' enter number of knots to add to this b-spline surface');
readln(n_knots_to_add);
for j := 1 to n_knots_to_add do
begin (* adding knots *)
   writeln(' enter 1 or 2 for 1st or 2nd parameter to add knot');
   readln(param);
   writeln(' enter knot value to be added');
   readln(x);
   rbsurf_to_bsurf(b,btemp);
   add_surf_knot(btemp,param,x,bnew);
   bsurf_to_rbsurf(bnew,b);
end;  (* adding knots *)
writeln
(' enter number of knots to try to remove from this b-spline surface ');
readln(n_knots_to_remove);
for j := 1 to n_knots_to_remove do
begin (* removing knots *)
   writeln
   (' enter 1 or 2 for 1st or 2nd parameter to remove knot');
   readln(param);
   writeln(' enter index of knot to be removed');
   if param = 1 then
   writeln(' index must be > ',b.degree1+1:1,' and < ',
   b.num_of_knots1 - b.degree1:1);
   if param = 2 then
   writeln(' index must be > ',b.degree2+1:1,' and < ',
   b.num_of_knots2 - b.degree2:1);
   readln(i);
   tol := exp(-20);
   rbsurf_to_bsurf(b,btemp);
   remove_surf_knot(btemp,param,i,tol,removed,bnew);
   writeln(' removed =',removed);
   if removed then bsurf_to_rbsurf(bnew,b);
end;  (* removing knots *)
      writeln(report,'**************************************');
      writeln(report,' input data for this case ');
      write_rbsurf(report,b);
      rbsurf_to_rppsurf(b,p);
      writeln(report,'--------------------------------------');
      cms('clrscrn',return_code);
      writeln(output,'--------------------------------------');
      writeln(report,' ');
      writeln(output,' ');
      writeln(report,' output data for this case ');
      writeln(output,' output data for this case ');
      write_rppsurf(report,p);
      write_rppsurf(output,p);
      writeln(' ');
      writeln(' enter m >= 1, n >= 1 to evaluate on an m x n grid ');
      read(m); read(n);
      if (m>=1) and (n>=1) then
      s_eval_and_compare(p,b,m,n);

writeln(' ');
writeln(' enter k > 0 if you want to test k repeated ');
writeln(' conversions between pp and b-spline.  ');
writeln(' enter 0 otherwise.  ');
readln(k);
if k > 0 then
begin (* iterative testing *)
for i := 1 to k do
begin (* i-th pair of conversions *)
   writeln(' starting conversion pair number ',i:1);
   rbsurf_to_rppsurf(b,p);
   rppsurf_to_rbsurf(p,b);
end;  (* i-th pair of conversions *)
cms('clrscrn',return_code);
writeln(report,' +++++++++++++++++++++++++++++++++++++++++++++++++ ');
writeln(report,' ');
writeln(report,' after ',i:2,' conversions, the result is:');
   write_rbsurf(report,b);
writeln(output,' +++++++++++++++++++++++++++++++++++++++++++++++++ ');
writeln(output,' ');
writeln(output,' after ',i:2,' conversions, the result is:');
   rbsurf_compare(b0,b);
end;  (* iterative testing *)
   end;  (* conversion *)
   until (b.degree1 <= 0) or (b.degree2 <= 0);
end;  (* procedure start_with_rbs_surface *)
(**********************************************************************)
(**********************************************************************)
(**********************************************************************)

(**********************************************************************)

(*      name:  main program ppbspln                                   *)

(*  category:  test dirvers                                           *)

(*   purpose:  prompt the user for what type of curve or surface      *)
(*             he wants to start with, and call the appropriate       *)
(*             "start_with_...." procedure.                           *)

(*     input:  no input parameters.  the user is prompted for input.  *)

(*    output:  no output parameters.  the appropriate procedure is    *)
(*             invoked.                                               *)

(**********************************************************************)
begin (* main program ppbspln *)
termin(input); termout(output);
rewrite(debug,'name=ppbspln.debug.a');
rewrite(report,'name=ppbspln.report.a');
findeps(machine_epsilon);
cms('clrscrn',return_code);
writeln;writeln;
writeln(' enter 1 if you want debug on.');
writeln(' enter 0 otherwise.');
readln(idebug);
repeat (* until i < 1 or i > 4 *)
cms('clrscrn',return_code);
writeln;
writeln
(' enter 1 to start with a piecewise rational polynomial curve');
writeln
(' enter 2 to start with a rational b-spline curve');
writeln
(' enter 3 to start with a piecewise rational polynomial surface');
writeln
(' enter 4 to start with a rational b-spline surface');
writeln
(' enter 0 (or any other integer) to stop');
readln(i);
cms('clrscrn',return_code);
case i of
        1 : start_with_rpp_curve;
        2 : start_with_rbs_curve;
        3 : start_with_rpp_surface;
        4 : start_with_rbs_surface;
otherwise   writeln(' end program ppbspln ');
end;  (* case statement *)
until (i < 1) or (i > 4);

end. (* main program ppbspln *)