#include <stdio.h>
#include <math.h>

//::::::::::::::
//aitken.c
//::::::::::::::
float aitken(int degree,float *coeff,float t)
/*
     uses  Aitken's algorithm to compute one coordinate
     value of a  Lagrange curve. Has to be called 
     for each coordinate  (x,y, and/or z) of data set.
     Input:   degree: degree of curve.
              coeff:  array with coordinates to be interpolated.
              t:      parameter value where to interpolate.
      Output: coordinate value.

      Note: we  assume a uniform knot sequence!
*/
 {
	int r,i;
	float t1;
	float coeffa[30];
   
	for (i=0; i<=degree; i++) coeffa[i]=coeff[i]; /* save input */
     

	for (r=1; r<= degree; r++)
	for (i=0; i<= degree - r; i++)
	{
	t1=(degree*t-i)/((float) r);
		coeffa[i]= (1.0-t1)* coeffa[i]  +   t1 * coeffa[i+1] ;
	}     

	return (coeffa[0]);
}

//::::::::::::::
//area.c
//::::::::::::::
float area(float *p1,float *p2,float *p3)
/* find area of 2D triangle p1,p2,p3 */

	//float p1[2], p2[2], p3[2];
{
	return( (p2[0]-p1[0])*(p3[1]-p1[1])-(p2[1]-p1[1])*(p3[0]-p1[0]))/2.0;
}
//::::::::::::::
//bessel_ends.c
//::::::::::::::
void bessel_ends(float *data, float *knot,int l)
/*	Computes B-spline points data[1] and data[l+]
	according to Bessel end condition.

input:	data:	sequence of data coordinates data[0] to data[l+2].
		Note that data[1] and data[l+1] are expected to
		be empty, as they will be filled by this routine.
	knot:	knot sequence
	l:	number of intervals
 
output:	data:	completed, as above.
*/
	//float data[], knot[];
	//int l;
{
	float alpha, beta; int i;
	
	if (l==1)
	{/*  This is not really Bessel, but then what do you do
	     when you have only one interval? -- make it linear!
	 */
		data[1]= (2.0*data[0] + data[3])/3.0;
		data[2]= (2.0*data[3] + data[0])/3.0;
	}

	else if (l==2)
	{
		/* beginning:    */
		alpha= (knot[2]-knot[1])/(knot[2]-knot[0]);
		beta = 1.0 - alpha;

		data[1]=(data[2]-alpha*alpha*data[0]-beta*beta*data[4])
		/(2.0*alpha*beta);
		data[1]=2.0*(alpha*data[0]+beta*data[1])/3.0 + data[0]/3.0;

		/* end:  */
		alpha= (knot[2]-knot[1])/(knot[2]-knot[0]);
		beta = 1.0 - alpha;

		data[3]=(data[2]-alpha*alpha*data[0]-beta*beta*data[4])
		/(2.0*alpha*beta);

		data[3]=2.0*(alpha*data[3]+beta*data[4])/3.0+data[4]/3.0;	
	}
	
	else

	{


		/* beginning:    */
		alpha= (knot[2]-knot[1])/(knot[2]-knot[0]);
		beta = 1.0 - alpha;

		data[1]=(data[2]-alpha*alpha*data[0]-beta*beta*data[3])
		/(2.0*alpha*beta);
		data[1]=2.0*(alpha*data[0]+beta*data[1])/3.0 + data[0]/3.0;


		/* end:  */
		alpha= (knot[l]-knot[l-1])/(knot[l]-knot[l-2]);
		beta = 1.0 - alpha;

		data[l+1]=(data[l]-alpha*alpha*data[l-1]-beta*beta*data[l+2])
		/(2.0*alpha*beta);

		data[l+1]=2.0*(alpha*data[l+1]+beta*data[l+2])/3.0+data[l+2]/3.0;

	}

}
//::::::::::::::
//bez_to_points.
//::::::::::::::
void bez_to_points(int degree,int npoints,float *coeff, float *points)
/*	Converts Bezier curve into point sequence. Works on
	one coordinate only.
 
	Input:   degree:  degree of curve.
	         npoints: # of coordinates to be generated. (counting
		          from 0!)
	         coeff:   coordinates of control polygon.
	Output:  points:  coordinates of points on curve.

	Remark: For a 2D curve, this routine needs to be called twice,
		once for the x-coordinates and once for y.
*/

	//int degree, npoints;
	//float coeff[],points[];

{	
	float t,delt;
	int i;
	float hornbez();
	int k;

	delt=1.0/(float)npoints;
	t=0.0;
	for(i=0; i<=npoints; i++) 
	{
		points[i]=hornbez(degree,coeff,t);
		t = t+delt;
	}
}
//::::::::::::::
//bez_to_power.c
//::::::::::::::
void bezier_to_power(int degree, float *bez, float *coeff)
/*	Converts Bezier form to power (monomial) form. Works on
	one coordinate only.
 
	Input:   degree:  degree of curve.
	         bez:	  coefficients of Bezier form
	Output:  coeff:   coefficients of power form.

	Remark: For a 2D curve, this routine needs to be called twice,
		once for the x-coordinates and once for y.
*/

	//int degree;
	//float coeff[], bez[];

{	
	float  i_factorial, n_r;
	int i;

	differences(degree,bez,coeff); /* compute forward differences    */
				       /* and store them in coeff.       */

					/* Note that i_factorial is int. */
					/* For high degrees: danger!     */
	coeff[0]=bez[0];
	i_factorial=1;	
	n_r=1;
	for (i=1; i<=degree; i++)
	{
		i_factorial=i_factorial*i;
		n_r= n_r*(degree-i+1);

		coeff[i]=n_r*coeff[i]/i_factorial;
	}


}


//::::::::::::::
//bspl_kappas.c
//::::::::::::::
extern FILE *outfile;

void bspl_kappas(float *bspl_x, float *bspl_y,float *bspl_w,float *knot,int l,int dense)
/*	writes curvatures of cubic rational  B-spline curve into 
	a file.
input:
	bspl_x,bspl_y:  2D rat. B-spline polygon
	bspl_w:	the  B-spline weights
	knot:   the knot sequence
	dense:	how many curvature values to compute per interval
	l:      no. of intervals
output:
	written into file outfile
*/

	//float bspl_x[],bspl_y[],bspl_w[],knot[];
	//int  dense,l;
{
	float bez_x[300],bez_y[300], bez_w[300];
	float bleftx[4],blefty[4],brightx[4],brighty[4];
	float wleft[4], wright[4],coeffx[4],coeffy[4],weight[4];
	float dist,t,delt,h,u,diff;
	int i,j,i3;
	float curvature_0(),curvature_1(),abs();

	/* first, convert B-spline to Bezier: */

	ratbspline_to_bezier(bspl_x,bspl_y,bspl_w,knot,l,
	                     bez_x,bez_y,bez_w);

	/* Now plot kappas for each interval. */

	fprintf(outfile,"%d  %d\n",dense*l,dense);
	for(i=0; i<l; i++)
	{	i3=i*3;
		u= knot[i];  /*starting value for i-th segment*/
		diff=knot[i+1]-knot[i];
		for(j=0;j<=3;j++)
		{       /* create i-th cubic Bezier curve: */ 
			coeffx[j]=bez_x[i3+j]; 
			coeffy[j]=bez_y[i3+j];
			weight[j]=bez_w[i3+j];
		}


		delt=1.0/(float)dense;

		for(t=0.0; t<0.5; t=t+delt)
		{
		subdiv(3,coeffx,weight,t,bleftx,brightx,wleft,wright);
		subdiv(3,coeffy,weight,t,blefty,brighty,wleft,wright);
		h=curvature_0(brightx,brighty,wright,3);

		/* to print file, we'll have to use the GLOBAL
		   parameter u!   */

		u=knot[i]+t*diff;
		fprintf(outfile,"%2.3f  %f\n",u,h);
		}

		for(t=0.5; t<0.999; t=t+delt)
		{
		subdiv(3,coeffx,weight,t,bleftx,brightx,wleft,wright);
		subdiv(3,coeffy,weight,t,blefty,brighty,wleft,wright);

		h=curvature_0(bleftx,blefty,wleft,3);
		h=0.0-h;
		/* minus sign since order of polygon 
		   traversal is reversed!
		*/
		u=knot[i]+t*diff;
		fprintf(outfile,"%2.3f    %f\n",u,h);
		}
	}
	/* print very last kappa: */
	t=1.0;
	h=curvature_1(coeffx,coeffy,weight,3);
	fprintf(outfile,"%2.3f    %f\n",knot[l],h);


}
//::::::::::::::
//c2_spline.c
//::::::::::::::
void c2_spline(float *knot,int l,float *data_x,float *data_y,float *bspl_x,float *bspl_y)
/* Finds the C2 cubic spline interpolant to
the data points in data_x, data_y.
Input:	knot:	the knot sequence  knot[0], ..., knot[l]
	l:	the number of intervals
	data_x, data_y: the data points data_x[0], ...,
		data[l+2]. Attention:
		data_x[1] and data_x[l+1] are filled
		by  Bessel end conditions and are
		thus ignored on input. Same for data_y.
Output:	bspl_x, bspl_y: the B-spline control polygon of
		the interpolant. Dimensions:
		bspl_x[0], ..., bspl_x[l+2]. Same for
		bspl_y.
*/
	//float knot[],data_x[],data_y[],bspl_x[],bspl_y[];
	//int l;
{
	float alpha[100], beta[100], gamma[100], up[100], low[100];
	set_up_system(knot,l,alpha,beta,gamma);

	l_u_system(alpha,beta,gamma,l,up,low);

	bessel_ends(data_x,knot,l);
	bessel_ends(data_y,knot,l);

	solve_system(up,low,gamma,l,data_x,bspl_x);
	solve_system(up,low,gamma,l,data_y,bspl_y);

}
//::::::::::::::
//check_flat.c
//::::::::::::::


int check_flat(float *bx,float *by,int degree,float tol)
/* Checks if a polygon is flat. If all points
   are closer  than tol to the connection of the 
   two endpoints, then it is flat. Crashes if the endpoints  
   are identical.

Input:	bx,by, degree: the Bezier curve
	tol:            tolerance

*/
	//float bx[],by[];
	//int degree; float tol;
{
	float dist;
	int i;
	float height();

	for (i=1; i<degree; i++)
	{
		dist = height(bx[i],by[i],bx[0],by[0],
			bx[degree],by[degree]);
		if (fabs(dist) > tol)  return(0);
	}
	
	return(1);
}
//::::::::::::::
//conic_weight.c
//::::::::::::::

float conic_weight(float *b0, float *b1, float *b2, float *p)
/* 
Input:	b0,b1,b2:	conic control polygon
	p:		point on conic
Output:	weight of b1 (assuming standard form).

For method, check Farin  14.5, 14.6.
*/

	//float b0[],b1[],b2[],p[];
{
	float area(),length_2();
	float tau0,  tau2 , t,t1, d, dist, weight;
	float b11[2], b10[2];

	tau0=area(b1,p,b2);
	tau2=area(b1,b0,p);
	return(area(b0,b2,p)/(2.0*sqrt(tau0*tau2))); /* weight of b1         */

}
//::::::::::::::
//curvature_0.c
//::::::::::::::

float curvature_0(float *bez_x,float *bez_y, float *weight,int degree)
/* computes curvature of rational  Bezier curve at t=0 */
	//float bez_x[], bez_y[], weight[];
	//int degree;
{
	float b0[2],b1[2],b2[2];
	float dist;
	float area();

	b0[0]=bez_x[0]; b1[0]=bez_x[1]; b2[0]=bez_x[2];
	b0[1]=bez_y[0]; b1[1]=bez_y[1]; b2[1]=bez_y[2];

	dist = sqrt( (b1[0]-b0[0])*(b1[0]-b0[0])+
	             (b1[1]-b0[1])*(b1[1]-b0[1])  );


	return (2.0*(degree-1)*weight[0]*weight[2]*area(b0,b1,b2)
	       /(degree*weight[1]*weight[1]*dist*dist*dist));
}
//::::::::::::::
//deboor.c
//::::::::::::::
float deboor(int degree,float *coeff,float *knot,float u,int i)
/*
	uses de Boor algorithm to compute one
	coordinate on B-spline curve for param. value u in interval i.
input:  degree:	polynomial degree of each piece of curve
	coeff:	B-spline control points
	knot:	knot sequence
	u:	evaluation abscissa
	i:	u's interval: u[i]<= u < u[i+1]
output:	coordinate value. 
*/

	//float coeff[],knot[];
	//float u;
	//int degree,i;
{
	int k,j;
	float t1,t2;
	float coeffa[30];  /* might need adjustment! */
	
	for (j=i-degree+1; j<=i+1; j++)coeffa[j]=coeff[j];

	for (k=1; k<= degree; k++)
	for ( j=i+1 ;j>=i-degree+k+1; j--)
	{
		t1= (knot[j+degree-k] - u )/(knot[j+degree-k]-knot[j-1]);
		t2= 1.0-t1;

		coeffa[j]=t1* coeffa[j-1]+t2* coeffa[j];
	}	
	return (coeffa[i+1]);
}




//::::::::::::::
//decas.c
//::::::::::::::
float decas(int degree,float *coeff,float t)
/*
     uses  de Casteljau to compute one coordinate
     value of a  Bezier curve. Has to be called 
     for each coordinate  (x,y, and/or z) of a control polygon.
     Input:   degree: degree of curve.
              coeff:  array with coefficients of curve.
              t:      parameter value.
      Output: coordinate value.
*/

	//float coeff[];
	//float t;
	//int degree;
{
	int r,i;
	float t1;
	float coeffa[30];  /* an auxiliary array. Change dim. if too small*/

	t1 = 1.0 - t;
	for(i=0; i<=degree; i++) coeffa[i]=coeff[i]; /* save input */

	for (r=1; r<= degree; r++)
		for (i=0; i<= degree - r; i++)
		{
			coeffa[i]= t1* coeffa[i]  +   t * coeffa[i+1] ;
		}

	return (coeffa[0]);
}


//::::::::::::::
//differences.c
//::::::::::::::
void differences(int degree,float *coeff,float *diffs)
/*
     Computes all forward differences Delta^i(b_0).
     Has to be called 
     for each coordinate  (x,y, and/or z) of a control polygon.
     Input:   degree: length (from 0) of coeff.
              coeff:  array of coefficients.
      Output: diffs:  diffs[i]= Delta^i(coeff[0]).
*/
//float coeff[],diffs[];
//int degree; 
{
	int r,i;
	float diffs1[20];

	for(i=0; i<=degree; i++) diffs1[degree-i]=coeff[i]; 

	for (r=1; r<= degree; r++)
	for (i=0; i<= degree - r; i++)
	diffs1[i]=  diffs1[i] - diffs1[i+1] ;  

	for(i=0;i<=degree;i++) diffs[i]=diffs1[degree-i];   

}

//::::::::::::::
//dist.c
//::::::::::::::


float dist(float *a,float *b)
/* finds distance between 2D points a   and b.

Input: a,b:	points;
*/
	//float a[], b[];
{
	return(sqrt((a[0]-b[0])*(a[0]-b[0])+
	            (a[1]-b[1])*(a[1]-b[1]) ) );
}


//::::::::::::::
//fair_bspline.c
//::::::::::::::
void fair_kr_bspline(float *bspl,float *knot,int l,int from,int to)
/* Fairs a cubic rational B-spline curve by knot removal/reinsertion.
 
Input:	bspl:      cubic B-spline control polygon (one coord.)
	knot:      knot sequence
	l:         no. of intervals
	from, to:  from where to where to fair
Output:	same as input, but hopefully fairer.
*/


	//float bspl[], knot[];
	//int l,from,to;
{
	int i, i3,  l2;
	float left,right, middle;
	float Knot[100]; /*to simulate multiple knots at ends */

	Knot[0]=Knot[1]=knot[0];

	for(i=0; i<=l; i++)Knot[i+2]=knot[i];
	Knot[l+3]=Knot[l+4]=knot[l];

	/* With new knot vector, bspl[i] is associated with Knot[i+1] */
	
	l2=l+2; 

	if(from <= 2) from = 2;
	if(to >=l) to=l;

	for (i=from; i<=to; i++) 
	{
		left=((Knot[i+2]-Knot[i-2])*bspl[i-1]
                       -(Knot[i+2]-Knot[i+1]  )*bspl[i-2]  )
                       /(Knot[i+1]-Knot[i-2] );

		
		right=((Knot[i+4]-Knot[i])*bspl[i+1]
                       -(Knot[i+1]-Knot[i]  )*bspl[i+2] )
                       /(Knot[i+4]-Knot[i+1] );

		
		bspl[i]=((Knot[i+3]-Knot[i+1])*left
		          +(Knot[i+1]-Knot[i-1])*right )
		          /(Knot[i+3]-Knot[i-1]) ;

	
	
	}
   
}
 
//::::::::::::::
//fair_surf.c
//::::::::::::::
void fair_surf(float *bspl[20],int lu,int lv,float *knot_u, float *knot_v)
/* 	Fairs B-spline control net.
Input:	bspl:   B-spline control net (one coordinate only)
	lu,lv:  no. of intervals in u- and v-direction
	knots_u, knots_v:   knot vectors in u- and v-direction
Output:	as input

*/
	//float bspl[20][20],knot_u[],knot_v[];
	//int  lu,lv;
{
	int  i,j,lu2,lv2;
	float b[20];

	lu2=lu+2; lv2=lv+2; 
 

	for(i=0;  i<= lu2; i++)
	{

		for(j=0; j<=lv2; j++)  b[j]=bspl[i][j]; /* grab i-th row*/
		fair_bspline(b,knot_v,lv,0,lv);         /* fair row     */
		for(j=0; j<=lv2; j++)  bspl[i][j]=b[j]; /* put row back */
	}

	for(j=0;  j<= lv2; j++)
	{

		for(i=0; i<=lu2; i++)  b[i]=bspl[i][j]; /* grab j-th col*/
		fair_bspline(b,knot_u,lu,0,lu);         /* fair col     */
		for(i=0; i<=lu2; i++)  bspl[i][j]=b[i]; /* put col back */
	}


}
//::::::::::::::
//gamma_spline.c
//::::::::::::::
void  gamma_spline(float *bspl,int l,float *gamma,float *bez)
/*From B-spline polygon and gammas, find interior
  Bezier points b3i+-1.
Input:	bspl:	control polygon (one coordinate only)
	l:	number of cubic pieces
	gamma:	gamma_i's
Output:	bez:	piecewise Bezier polygon,  but without the 
		junction points b3i
*/
	//int l;
	//float bspl[],gamma[],bez[];
{
	int i,i3,l3;

	bez[0]=bspl[0];
	bez[1]=bspl[1];
	if(l>1)
	{
	bez[2]=(gamma[2]*bspl[1]
		+bspl[2]    )/(gamma[2]+1.0);
	}   
	if(l>2)
	{
	bez[4]=( (1.0+gamma[3])*bspl[2]
	        +gamma[2]*bspl[3]  )/(1.0+gamma[2]+gamma[3]);
	}
/* Now the main part:   */

	for(i=2; i<l-1; i++)
	{		
	i3=3*i;
	bez[i3-1]=( gamma[i+1]*bspl[i]
	         +  (gamma[i]+1.0)*bspl[i+1])/(1.0+gamma[i]+gamma[i+1]);

	bez[i3+1]=( (gamma[i+2]+1.0)*bspl[i+1]
	         +  gamma[i+1]*bspl[i+2]  )/(1.0+gamma[i+1]+gamma[i+2]);

	}
	
/* The last points need special attention:  */

	l3=l*3;
	if(l>2)
	{
	bez[l3-4]=( gamma[l]*bspl[l-1]
	         +(1.0+gamma[l-1])*bspl[l]  )/(1.0+gamma[l-1]+gamma[l]);
	}
	if(l>1)
	{
	bez[l3-2]=( bspl[l]
		 +  gamma[l]*bspl[l+1]  )/(gamma[l]+1.0);
	}
	bez[l3-1]= bspl[l+1];
	bez[l3]  = bspl[l+2]; 

}
//::::::::::::::
//height.c
//::::::::::::::


float height(float px,float py,float ax,float ay,float bx,float by)
/* find the height of point (px,py) over
   straight line thru (ax,ay) and (bx,by). 
*/
	//float px,py,ax,ay,bx,by;
{
	float det,length;
	det = px*ay + ax*by + bx*py -ay*bx - by*px - py*ax;

	length = sqrt( (bx-ax)*(bx-ax) + (by-ay)*(by-ay) );

	if(length < 0.000001)   /* height doesn't make sense, hence: */
		return( sqrt( (px-ax)*(px-ax) + (py-ay)*(py-ay) ) );
	else
	return (det/length);
}

//::::::::::::::
//hornbez.c
//::::::::::::::
float hornbez(int degree,float *coeff,float t)
/*
     uses  Horner's scheme to compute one coordinate
     value of a  Bezier curve. Has to be called 
     for each coordinate  (x,y, and/or z) of a control polygon.
     Input:   degree: degree of curve.
              coeff:  array with coefficients of curve.
              t:      parameter value.
      Output: coordinate value.
*/
	//int degree;
	//float coeff[];
	//float t;
{
	int i;
	int n_choose_i;          /* shouldn't be too large! */
	float fact,t1,aux;

	t1=1.0 - t;  fact=1.0;
	n_choose_i=1;

	aux=coeff[0]*t1;          /* starting the evaluation
                                  loop */
	for(i=1; i<degree; i++)
	{
		fact=fact*t;
		n_choose_i=n_choose_i*(degree-i+1)/i;  /* always int! */
		aux=(aux + fact*n_choose_i*coeff[i])*t1;
	}
	aux = aux + fact*t*coeff[degree];

	return aux;
}

//::::::::::::::
//horner.c
//::::::::::::::
float horner(int degree,float *coeff,float t)
/*
     uses  Horner's scheme to compute one coordinate
     value of a curve in power form. Has to be called 
     for each coordinate  (x,y, and/or z) of a control polygon.
     Input:   degree: degree of curve.
              coeff:  array with coefficients of curve.
              t:      parameter value.
      Output: coordinate value.
*/
	//int degree;
	//float coeff[];
	//float t;
{
	int i;
	float aux;

	aux=coeff[degree];      

	for(i=degree; i>0; i--)
		aux= t*aux + coeff[i-1];


	return(aux);
}


//:::::::::::::::::
//intersect_lines.c
//:::::::::::::::::

void intersect_lines(float *a,float *b,float *v,float *w,float *p)
/* finds the intersection p of the two straight lines a + tv and
	b+tw. No check for infinite p is made.
*/
//	float a[2],b[2],v[2],w[2],p[2];

{
	float D,t;
	float diff[2];

	diff[0]=a[0]-b[0];  diff[1]=a[1]-b[1];

	D= w[1]*v[0] - w[0]*v[1]; 


	t= (diff[1]*w[0] - diff[0]*w[1]  )/D;

	p[0]=a[0]+t*v[0];  p[1]=a[1]+t*v[1];

}

//::::::::::::::::::
//intersect_points.c
//::::::::::::::::::

void intersect_points(float *a,float *b,float *c,float *d,float *p)
/* finds the intersection p of the two straight lines [a,b] and
	[c,d]. No check for infinite p is made.
*/
	//float a[2],b[2],c[2],d[2],p[2];

{
	float diff1[2],diff2[2];

	diff1[0]=b[0]-a[0]; diff1[1]=b[1]-a[1];
	diff2[0]=d[0]-c[0]; diff2[1]=d[1]-c[1];

	intersect_lines(a,d,diff1,diff2,p);

}


//::::::::::::::
//intersect.c
//::::::::::::::

extern FILE *outfile, *psfile;

void intersect(float *bx,float *by,float *w,int degree,float tol)
/* Intersects bezier curve with x-axis  by  adaptive subdivision. 
   Subdivision is controlled by tolerance tol. There is
   no check for stack depth! Intersection points are not found  in
   `natural' order.  Results are written into file outfile.
	Input: bx,by,w:	rational Bezier curve
		degree:	its degree
		tol: 	accuracy for results
 */

	//int degree; float tol;
	//float bx[],by[],w[];
{
	float a[2],b[2],c[2],d[2],p[2];
	float bxleft[20],byleft[20],bxright[20],
		byright[20],wleft[20],wright[20];
	float box[4],xx,yy;
	int j;
	int check_flat();
	
	/********* does minmaxbox contain x-axis? **********/

	minmax(bx,by,degree,box);



	if (box[2]*box[3] > 0.0) return;

	/********* is polygon flat enough? ********/

	if( check_flat(bx,by,degree,tol)== 1 )
		{ 
		a[0]=bx[0];      a[1]=by[0];        /* intersect curve       */
		b[0]=bx[degree]; b[1]=by[degree];   /* endpoints with        */
		c[0]=0.0;        c[1]=0.0;          /* x- axis. Intersection */
		d[0]=1.0;	 d[1]=0.0;          /* is in p[0], p[1].     */
		intersect_points(a,b,c,d,p);

		fprintf(outfile,"%f  %f\n",p[0],p[1]);
		}
	/******** if not flat enough, subdivide: ********/
	else
		{

		subdiv(degree,bx,w,0.5,bxleft,bxright,wleft,wright);
		subdiv(degree,by,w,0.5,byleft,byright,wleft,wright);

	/******** ... and repeat for each half: **********/

		intersect(bxleft, byleft, wleft, degree,tol);
		intersect(bxright,byright,wright,degree,tol);
		}
}
//::::::::::::::
//l_u_system.c
//::::::::::::::
void l_u_system(float *alpha,float *beta,float *gamma,int l,float *up,float *low)
/*	perform LU decomposition of tridiagonal system with
	lower diagonal alpha, diagonal beta, upper diagonal gamma.

input:	alpha,beta,gamma: the coefficient matrix entries
	l:	matrix size [0,l]x[0,l]
	low:	L-matrix entries
	up:	U-matrix entries
*/
	//float alpha[],beta[],gamma[],low[],up[];
	//int l;
{
	int i;

	up[0]=beta[0];
	for(i=1; i<=l; i++)
	{
		low[i]=alpha[i]/up[i-1];
		up[i] =beta[i]-low[i]*gamma[i-1];
	}
}
//::::::::::::::
//length_2.c
//::::::::::::::

float length_2(float *a)
/* computes length of vector a
input:	a
output:	length
*/
//float a[];
{float x;
	x=sqrt(a[0]*a[0]+a[1]*a[1]);
	return(x);
}

//::::::::::::::
//minmax.c
//::::::::::::::
 float minmax(float *bx,float *by,int degree,float *value)
	/* find minmax box of 2D polygon bx,by */

	//int degree;
	//float bx[],by[],value[];
{
	int i;
	float xmin,xmax,ymin,ymax;

	xmin= bx[0]; ymin = by[0];
	xmax= bx[0]; ymax = by[0] ;

	for (i=0; i<=degree; i++)
	{
		if (bx[i]<xmin)xmin=bx[i];
		if (bx[i]>xmax)xmax=bx[i];

		if (by[i]<ymin)ymin=by[i];
		if (by[i]>ymax)ymax=by[i];
	
	}
		value[0]=xmin;
		value[1]=xmax;
		value[2]=ymin;
		value[3]=ymax;

}

	

//::::::::::::::
//minmax_surf.c
//::::::::::::::
 minmax_surf(float *bx,float *by,int degree_u,int degree_v,float *value)
	/* find minmax box of 2D polygon bx,by */

	//int degree_u,degree_v;
	//float bx[20][20],by[20][20],value[];
{
	int i,j;
	float xmin,xmax,ymin,ymax;

	xmin= bx[0][0]; ymin = by[0][0];
	xmax= bx[0][0]; ymax = by[0][0] ;

	for (i=0; i<=degree_u; i++)
	for (j=0; j<=degree_v; j++)
	{
		if (bx[i][j]<xmin)xmin=bx[i][j];
		if (bx[i][j]>xmax)xmax=bx[i][j];

		if (by[i][j]<ymin)ymin=by[i][j];
		if (by[i][j]>ymax)ymax=by[i][j];
	
	}
		value[0]=xmin;
		value[1]=xmax;
		value[2]=ymin;
		value[3]=ymax;

}

	
//::::::::::::::
//netcoons.c
//::::::::::::::
void netcoons(float *net[20],int rows,int columns)
/* Uses bilinear Coons blending to complete a control
net of which only the four boundary polygons are given.
Works for one coordinate only.
Input:	net:	control net - only boundaries are considered as input.
	rows, columns: dimensions of net.
Output:	net:	the completed net, with the old boundaries.
*/
	//float net[20][20]; 
	//int rows, columns; /*If either of these is >20, adjust dimensions*/
 //                          /*in previous statement!                      */
{
	int i,j;
	float u,v,u1,v1;

	for(i=0; i<=rows;    i++)
	for(j=0; j<=columns; j++)
	{

		u=((float)i)/rows;    u1=1.0-u;
		v=((float)j)/columns; v1=1.0-v;

		net[i][j]=u1*net[0][j] + u*net[rows][j]
			 +v1*net[i][0] + v*net[i][columns]
			 -u1*v1*net[0][0]    - u1*v*net[0][columns]
			 -u *v1*net[rows][0] - u *v*net[rows][columns];

	}
}

//::::::::::::::
//parameters.c
//::::::::::::::


void parameters(float *data_x,float *data_y,int l,float *knot)
/*	Finds a centripetal parametrization for a given set
	of 2D data points. 
Input:	data_x, data_y:	input points, numbered from 0 to l+2. 
	l:		number of intervals.
Output:	knot:		knot sequence. Note: not (knot[l]=1.0)!
Note:	data_x[1], data_x[l+1] are not used! Same for data_y.
*/
	//float data_x[], data_y[],knot[];
	//int l;

{
	int i,  l1,l2;     /* In the following, special care must be    */
	float delta;       /* at the ends because of the data structure */
	                   /* used. See note above.                     */
	l1= l-1; l2 = l+2;
	knot[0] = 0.0;     /* initialize -- arbitrary */

	delta=sqrt((data_x[2]-data_x[0])*(data_x[2]-data_x[0])
		         +(data_y[2]-data_y[0])*(data_y[2]-data_y[0]));
	knot[1]= sqrt(delta);                /* leave out this sqrt
		                                if you want chord length.*/

	for(i=2; i<l; i++)
	{
	        delta=sqrt((data_x[i+1]-data_x[i])*(data_x[i+1]-data_x[i])
		         +(data_y[i+1]-data_y[i])*(data_y[i+1]-data_y[i]));

		knot[i]=knot[i-1]+sqrt(delta);/*leave out this sqrt
		                                if you want chord length.*/
	}

	delta=sqrt((data_x[l2]-data_x[l])*(data_x[l2]-data_x[l])
		         +(data_y[l2]-data_y[l])*(data_y[l2]-data_y[l]));
	knot[l]= knot[l-1]+sqrt(delta);      /* leave out this sqrt
		                                if you want chord length.*/

}

//::::::::::::::
//plot_ratsurf.c
//::::::::::::::

extern FILE *psfile;    /* for postscript output */

void plot_ratsurf(float *bx[20],float *by[20],float *bw[20],
	int degree_u,int degree_v,int u_points,int v_points,float scale_x,float scale_y)
/*	plots v_points isoparametric 
        curves of the  rat Bez surface, each with u_points
        points on it.
	
Input:  bx, by:		 arrays with x- and y- coordinates of
	       		 control net.
        degree_u,degree_v: degrees in u- and v- direction
	u_points,v_points: plot resolution
	scale_x,scale_y:   scale factor for postscript.

Output: postscript file 
*/

	//int degree_u,degree_v,u_points,v_points;
	//float bx[20][20], by[20][20],bw[20][20],scale_x,scale_y;

{
	int i,j,u_points1;
	float u,v,del_u,del_v;
	float aux_x[5],aux_y[5],aux_w[5], arrax[5],array[5],arraw[5];
	float points_x[300], points_y[300];
	int iv;
	float ratbez(),decas(),hornbez();


	del_v= 1.0/((float)v_points);
	u_points1=u_points+1;
	v=0.0;

	for(iv=0; iv<= v_points; iv++)/* loop over all curves to be drawn*/
	{

		for(i=0;i<=degree_u; i++)/* create aux. arrays*/
		{
		for(j=0; j<= degree_v; j++)
			{
			arraw[j]=bw[i][j];
			arrax[j]=bx[i][j]*arraw[j];
			array[j]=by[i][j]*arraw[j];
			}

		aux_w[i]=hornbez(degree_v,arraw,v);
		aux_x[i]=hornbez(degree_v,arrax,v)/aux_w[i];
		aux_y[i]=hornbez(degree_v,array,v)/aux_w[i];

		}
		
		ratbez_to_points(degree_u,u_points,aux_x,aux_w,points_x);
		ratbez_to_points(degree_u,u_points,aux_y,aux_w,points_y);

	/******************* now ps: *******************************/	
			
		fprintf(psfile,"newpath\n");

		fprintf(psfile,"%f %f moveto\n",
			scale_x*points_x[0], scale_y*points_y[0]);


		for(j=1; j<= u_points; j++)
			fprintf(psfile,"%f %f lineto\n",scale_x*points_x[j], scale_y*points_y[j]);

		fprintf(psfile,"stroke\n");


		v = v+del_v;

	}
}


//::::::::::::::
//psplot_curve.c
//::::::::::::::


extern FILE *psfile;


void psplot_curve(float *bez_x,float *bez_y,float *bez_w,float *bspl_x,float *bspl_y,
	int l,int degree,float *value,int pol,int dense,float scale_x,float scale_y)
/* plots piecewise rational cubic  Bezier curve, dense points each segment,
  into the postscript-file   psfile. B-spline polygon can also be plotted.
Input:	bez_x, bez_y, bez_w: piecewise control polygons and weights
	bspl_x,bspl_y:B-spline controll polygon 
	l:            number of curve pieces
	degree:       degree of each piece
	value:	      minmax-box of polygon. xmin,xmax,...
	pol:	      which polygon to plot:
		      0: nothing, 1: Bezier, 2: B-spline,
	dense:	      how many points per segment
*/
	//float bez_x[],bez_y[],bez_w[],value[],bspl_x[],bspl_y[],scale_x, scale_y;
	//int l,degree,pol,dense;
{
	int i,k,j, ldeg,l2,l31,i3;
	float weight[5];

	float points_x[70],points_y[70], q_x[60], q_y[60];
	float aux_x[5], aux_y[5],aux_w[5];
	

	l2 = l+2;

	/* macro to mark B-spline vertex:*/
		fprintf(psfile,"/blackbox\n");
		fprintf(psfile,"{0 setgray\n");
		fprintf(psfile,"-3.5 -3.5 rlineto\n");
		fprintf(psfile,"7 0 rlineto\n");
		fprintf(psfile,"0 7 rlineto\n");
		fprintf(psfile,"-7 0 rlineto\n");
		fprintf(psfile,"0 -7 rlineto\n");
		fprintf(psfile,"closepath fill} def\n");
	

		fprintf(psfile,"/whitebox\n");
		fprintf(psfile,"{1 setgray\n");
		fprintf(psfile,"-3 -3 rlineto\n");
		fprintf(psfile,"6 0 rlineto\n");
		fprintf(psfile,"0 6 rlineto\n");
		fprintf(psfile,"-6 0 rlineto\n");
		fprintf(psfile,"0 -6 rlineto\n");
		fprintf(psfile,"closepath fill} def\n");
	
	/* macro to mark weight point:*/
		fprintf(psfile,"/diamond\n");
		fprintf(psfile,"{0  -5 rlineto\n");
		fprintf(psfile,"5 5 rlineto\n");
		fprintf(psfile,"-5 5 rlineto\n");
		fprintf(psfile,"-5 -5 rlineto\n");
		fprintf(psfile,"5 -5 rlineto\n");
		fprintf(psfile,"closepath} def\n");


/* plot curve */

	fprintf(psfile,"3 setlinewidth \n");
	

	/* for each segment, convert to points  and plot: */
	ldeg=degree*l;
	for(i=0; i< ldeg; i+=degree)
	{

	for(k=0; k<=degree; k++)
			{
			aux_x[k]=bez_x[i+k];
			aux_y[k]=bez_y[i+k];
			aux_w[k]=bez_w[i+k];
			}


		ratbez_to_points(degree,dense,aux_x,aux_w,points_x);
		ratbez_to_points(degree,dense,aux_y,aux_w,points_y);

		fprintf(psfile,"newpath \n");

		fprintf(psfile,"%f %f moveto\n",
			scale_x*(points_x[0]-value[0]),
			scale_y*(points_y[0]-value[2]));
		
		for( j=1; j<=dense ; j++)
		{
			fprintf(psfile,"%f %f lineto\n",
			scale_x*(points_x[j]-value[0]),
			scale_y*(points_y[j]-value[2]));				}
 		fprintf(psfile,"stroke \n");

	}

	



	/**** fill in junction points  black: **********/
	fprintf(psfile,"1 setlinewidth \n");


	for(j=0; j<= ldeg; j=j+degree)
			{
			fprintf(psfile,"%f %f %f %f %f arc fill\n",
			scale_x*(bez_x[j]-value[0]), 
			scale_y*(bez_y[j]-value[2]),
		 	3.0, 0.0, 360.0);
			}


if (pol == 0)
	{
		fprintf(psfile,"showpage \n");
		 return;
	}


if (pol>0) /* plot B-spline polygon: */
	{
		fprintf(psfile,"0.0 setgray\n ");
		fprintf(psfile,"2 setlinewidth\n ");

		fprintf(psfile,"newpath \n");
		fprintf(psfile,"%f %f moveto\n",
			scale_x*(bez_x[0]-value[0]), 
			scale_y*(bez_y[0]-value[2]));

		for(j=1; j<=  l2; j++)
			{
			fprintf(psfile,"%f %f lineto\n",
			scale_x*(bspl_x[j]-value[0]), 
			scale_y*(bspl_y[j]-value[2]));
			}
		fprintf(psfile,"stroke \n");



	/***** plot  control vertices: ************/


	for(j=0; j<=l2; j++)
	{
		fprintf(psfile,"%f %f moveto\n",
		scale_x*bspl_x[j], scale_y*bspl_y[j]);
		fprintf(psfile,"blackbox \n");

		if (j!=0 && j!=l2)
		{
			fprintf(psfile,"%f %f moveto\n",
			scale_x*(bspl_x[j]-value[0]),
			 scale_y*(bspl_y[j]-value[2]));
			fprintf(psfile,"whitebox \n");
		}
	}

	}

		
	
	/***** plot  bezier control polygon if pol=1: ************/
if (pol == 1)
	{
		fprintf(psfile,"0.0 setgray\n ");
		fprintf(psfile,"0.5 setlinewidth\n ");

		fprintf(psfile,"newpath \n");
		fprintf(psfile,"%f %f moveto\n",
			scale_x*(bez_x[0]-value[0]),
			 scale_y*(bez_y[0]-value[2]));

		for(j=0; j<=  ldeg; j++)
			{
			fprintf(psfile,"%f %f lineto\n",
			scale_x*(bez_x[j]-value[0]), 
			scale_y*(bez_y[j]-value[2]));
			}
		fprintf(psfile,"stroke \n");



	/***** plot  control vertices: ************/

	fprintf(psfile,"0.0 setgray\n ");

	for(j=1; j< ldeg ; j++)
			{
		if(j !=(j/degree)*degree)
			fprintf(psfile,"%f %f %f %f %f arc stroke\n",
			scale_x*(bez_x[j]-value[0]), 
			scale_y*(bez_y[j]-value[2]),
		 	3.0, 0.0, 360.0);
			}

	fprintf(psfile,"1.0 setgray\n ");

	for(j=1; j< ldeg; j++)
			{
		if( j!=(j/degree)*degree)	
		fprintf(psfile,"%f %f %f %f %f arc fill\n",
			scale_x*(bez_x[j]-value[0]), 
			scale_y*(bez_y[j]-value[2]),
		 	2.80, 0.0, 360.0);
			}


/*		fprintf(psfile,"showpage\n");  
		return;
*/
	


	/***** plot  weight points: ************

	fprintf(psfile,"0.3 setgray\n ");
	l31 = 3*l -1;

	for (j=0; j<= l31; j++)
	{
		q_x[j]=(bez_w[j]*bez_x[j] + bez_w[j+1]*bez_x[j+1])
			/(bez_w[j] + bez_w[j+1]);

		q_y[j]=(bez_w[j]*bez_y[j] + bez_w[j+1]*bez_y[j+1])
			/(bez_w[j] + bez_w[j+1]);
	}
	for(j=0; j<=l31; j++)
	{
		fprintf(psfile,"%f %f moveto\n",
		scale_x*(q_x[j]-value[0]), scale_y*(q_y[j]-value[2]));
		fprintf(psfile,"diamond fill\n");

	}
	
	************************** weight points plotted  */
	}


	fprintf(psfile,"showpage\n");
	return;

}		
//::::::::::::::
//psplot_net.c
//::::::::::::::



extern FILE *psfile;

void psplot_net(int lu,int lv,float *bx[20],float *by[20],int step_u,int step_v,float scale_x,float scale_y,float *value)
/*	plots  control net     	into postscript-file.
Input:	lu,lv:		dimensions of net
	bx,by:		net vertices
	step_u,step_v:	subnet sizes (e.g. both=3 for pw bicubic net)
	scale_x,scale_y:scale factors for ps
	value:		window size in world coords
*/
	//int lu,lv,step_u,step_v;
	//float bx[20][20], by[20][20],value[4],scale_x,scale_y;

{
	int lv1;
	int i,j;

	
		fprintf(psfile,"0.0 setgray \n");
		fprintf(psfile,"1 setlinewidth\n");


	lv1= lv+1;
	for(i=0; i<= lu; i++)
	{
		fprintf(psfile,"newpath\n");
		fprintf(psfile," %f %f moveto\n",
			scale_x*(bx[i][0]-value[0]),
			scale_y*(by[i][0]-value[2]));
		for(j=1; j<= lv; j++)
			fprintf(psfile,"%f %f lineto\n",
			scale_x*(bx[i][j]-value[0]), 
			scale_y*(by[i][j]-value[2]));
			
		fprintf(psfile,"stroke\n");
	}

	for(j=0;j<= lv;j++)
	{
		fprintf(psfile,"newpath\n");
		fprintf(psfile," %f %f moveto\n",
			scale_x*(bx[0][j]-value[0]),
			scale_y*(by[0][j]-value[2]));

		for(i=1; i<= lu; i++) 
			fprintf(psfile," %f %f lineto\n",
			scale_x*(bx[i][j]-value[0]),
			scale_y*(by[i][j]-value[2]));
		
		fprintf(psfile,"stroke\n");
	}

	/***** now mark control points *****/


	for(i=0; i<= lu; i++)
	{
		for(j=0; j<= lv; j++)
			fprintf(psfile,"%f %f %f %f %f arc stroke\n",
			scale_x*(bx[i][j]-value[0]),
			 scale_y*(by[i][j]-value[2]),
			3.0, 0.0, 360.0);		
	}


	fprintf(psfile,"1 setgray\n"); /* to white out circles */

	for(i=0; i<= lu; i++)
	{
		for(j=0; j<= lv; j++)
			fprintf(psfile,"%f %f %f %f %f arc fill\n",
			scale_x*(bx[i][j]-value[0]),
			 scale_y*(by[i][j]-value[2]),
			2.9, 0.0, 360.0);		
	}


	/* subnet corners in black: */
	fprintf(psfile,"0.0 setgray\n");


	for(i=0; i<= lu; i+= step_u)
	{
		for(j=0; j<= lv; j+= step_v)
			fprintf(psfile,"%f %f %f %f %f arc fill\n",
			scale_x*(bx[i][j]-value[0]), 
			scale_y*(by[i][j]-value[2]),
			2.9, 0.0, 360.0);		
	}
	


}

//::::::::::::::
//ratbez.c
//::::::::::::::


float ratbez(int degree,float *coeff,float *weight,float t)
/*
	uses rational de casteljau to compute
	point on ratbez curve for param. value t.
*/
	//float coeff[],weight[];
	//float t;
	//int degree; 
{
	int r,i;
	float t1,t2,ww1,ww2;
	float coeffa[20], weighta[20];
	t1 = 1.0 - t;

	
	for (i=0;i<=degree; i++)
	{
		coeffa[i]=coeff[i];
		weighta[i]=weight[i];
	}

	for (r=1; r<= degree; r++)
	for (i=0; i<= degree - r; i++)
	{
	/*	t1= (i+r-degree*t)/(float)r;
		t2= 1.0-t1;  ... makes it interpolatory ...
	*/
		t1=1.0-t; t2=t;
		ww1 = weighta[i]; ww2= weighta[i+1];
		weighta[i] = t1*ww1 + t2*ww2;
		coeffa[i]  = ( t1*ww1 * coeffa[i]
			    + t2*ww2 * coeffa[i+1] ) / weighta[i];

	}	
	return (coeffa[0]);
}
//::::::::::::::
//ratbez_kappas.c
//::::::::::::::


extern FILE *outfile;

void ratbez_kappas(float *bez_x,float *bez_y,float *bez_w,int l,int dense)
/*	writes curvatures of  piecewise cubic rational curve into 
	a file.
input:
	bez_x,bez_y:  2D rat. cubic spline polygon
	bez_w:	the  weights
	dense:	how many curvature values to compute per interval
	l:      no. of intervals
output:		 written into file outfile
Note:		this plots curvatures vs. a uniform u-spacing
*/

	//float bez_x[],bez_y[],bez_w[];
	//int  dense,l;
{
	float bleftx[4],blefty[4],brightx[4],brighty[4];
	float wleft[4], wright[4],coeffx[4],coeffy[4],weight[4];
	float dist,t,delt,h,u,diff;
	int i,j,i3;
	float curvature_0(),curvature_1(),abs();

	/* Plot kappas for each interval. */
	fprintf(outfile,"%d  %d\n",dense*l,dense);
	for(i=0; i<l; i++)
	{	i3=i*3;
		u= (float)i;  /*starting value for i-th segment*/
		diff=1.0;
		for(j=0;j<=3;j++)
		{       /* create i-th cubic Bezier curve: */ 
			coeffx[j]=bez_x[i3+j]; 
			coeffy[j]=bez_y[i3+j];
			weight[j]=bez_w[i3+j];
		}


		delt=1.0/(float)dense;

		for(t=0.0; t<0.5; t=t+delt)
		{
		subdiv(3,coeffx,weight,t,bleftx,brightx,wleft,wright);
		subdiv(3,coeffy,weight,t,blefty,brighty,wleft,wright);
		h=curvature_0(brightx,brighty,wright,3);
		/* to print file, we'll have to use the GLOBAL
		   parameter u!   */
		u=(float)i+t*diff;
		fprintf(outfile,"%2.3f  %f\n",u,h);
		}

		for(t=0.5; t<0.999; t=t+delt)
		{
		subdiv(3,coeffx,weight,t,bleftx,brightx,wleft,wright);
		subdiv(3,coeffy,weight,t,blefty,brighty,wleft,wright);

		h=curvature_0(bleftx,blefty,wleft,3);
		h=0.0-h;
		/* minus sign since order of polygon 
		   traversal is reversed!
		*/
		u=(float)i+t*diff;
		fprintf(outfile,"%2.3f    %f\n",u,h);
		}
	}
	/* print very last kappa: */
	t=1.0;
	h=curvature_1(coeffx,coeffy,weight,3);
	fprintf(outfile,"%2.3f    %f\n",(float)l,h);


}
//:::::::::::::::::
//ratbez_to_points.c
//:::::::::::::::::
void  ratbez_to_points(int degree,int npoints,float *coeff,float *coeff_w,float *points)
/*	Converts rational Bezier curve into point sequence. Works on
	one coordinate only.
 
	Input:   degree:  degree of curve.
	         npoints: # of coordinates to be generated. (counting
		          from 0!)
	         coeff:   coordinates of control polygon.
	Output:  points:  coordinates of points on curve.

	Remark: For a 2D curve, this routine needs to be called twice,
		once for the x-coordinates and once for y.
*/

	//int degree, npoints;
	//float coeff[],coeff_w[],points[];

{	
	float t,delt;
	int i;
	float ratbez();

	delt=1.0/(float)npoints;
	t=0.0;
	for(i=0; i<=npoints; i++) 
	{
	points[i]=ratbez(degree,coeff,coeff_w,t);  
	t = t+delt;
	}
}



//::::::::::::::
//reparam.c
//::::::::::::::
void reparam(float *wold,int degree,float s,float *wnew)
/*	reparametrizes ratbez curve: only the weights,
	stored in wold, are changed. New weights are in
	wnew. Parametrization is determined by shoulder
	point s. For s=0.5, nothing changes. Also,
	s should be in (0,1).
*/

	//float wold[],wnew[];
	//int degree;	float s;
{
	int i;
	float c;

	c = s / (1.0 - s);

	wnew[0] = wold[0];

	for (i=1; i<= degree; i++)
		wnew[i] = c*wnew[i-1]*wold[i]/wold[i-1];
}
//::::::::::::::
//set_up_system.c
//::::::::::::::
void set_up_system(float *knot,int l,float *alpha,float *beta,float *gamma)
/*	given the knot sequence, the linear system 
	for clamped end condition B-spline interpolation
	is set up.
input:	knot:	knot sequence
	points:	points to be interpolated
	l:	number of intervals
output:	alpha,beta,gamma: 1-D arrays that constitute
		the elements of the interpolation matrix.

	Note: no data points needed so far! 
*/
	//int l;
	//float knot[],alpha[],beta[],gamma[];
{
	int i,l1;
	float delta_im2,delta_im1, delta_i, delta_ip1,sum;
	l1=l-1;

	/* some special cases: */
	if(l==1)
	{alpha[0]=0.0; alpha[1]=0.0; beta[0]=1.0; beta[1]=1.0;
	 gamma[0]=0.0; gamma[1]=0.0; return;}

	if(l==2)
	{
		beta[0]=1.0;
		delta_im1=(knot[1]-knot[0]);
		delta_i  =(knot[2]-knot[1]);
		delta_ip1=(knot[3]-knot[2]);
		sum = delta_im1+delta_i;

		alpha[1]=delta_i*delta_i/sum;
		beta[1] =(delta_i*delta_im1)/sum
                + delta_im1*delta_i/ sum;
		gamma[1]=delta_im1*delta_im1/sum;

		alpha[1]=alpha[1]/sum;
		beta[1] =beta[1]/sum;
		gamma[1]=gamma[1]/sum;

		beta[2]=1.0;
		alpha[2]=0.0;
		gamma[2]=0.0;
		return;
	}

	/* the rest does the cases l>2.  */

	delta_im1=(knot[1]-knot[0]);
	delta_i  =(knot[2]-knot[1]);
	delta_ip1=(knot[3]-knot[2]);
	sum = delta_im1+delta_i;

	beta[0]=1.0; gamma[0]=0.0;

	alpha[1]=delta_i*delta_i/sum;
	beta[1] =(delta_i*delta_im1)/sum
                + delta_im1*(delta_i+delta_ip1)
		/ (sum+delta_ip1);
	gamma[1]=delta_im1*delta_im1/(sum+delta_ip1);

	alpha[1]=alpha[1]/sum;
	beta[1] =beta[1]/sum;
	gamma[1]=gamma[1]/sum;

	/*Now for the main loop:   */
	for(i=2; i<l1; i++)
	{
		/* compute delta_i_minus_2,...  */
		delta_im2=(knot[i-1]-knot[i-2]);
		delta_im1=(knot[i]  -knot[i-1]);
		delta_i  =(knot[i+1]-knot[i]);
		delta_ip1=(knot[i+2]-knot[i+1]);
	
		sum = delta_im1+delta_i;

		alpha[i]=delta_i*delta_i/(delta_im2 +sum);
		beta[i] =delta_i*(delta_im2+delta_im1) 
			/(delta_im2 + sum)
			+
			 delta_im1*(delta_i+delta_ip1)
			/(sum + delta_ip1);
		gamma[i]=delta_im1*delta_im1
			/(sum + delta_ip1);
		
		alpha[i]=alpha[i]/sum;
		beta[i] =beta[i]/sum;
		gamma[i]=gamma[i]/sum;
	}

	/*  special care at the end:  */
	delta_im2=knot[l-2]-knot[l-3];
	delta_im1=knot[l1]-knot[l-2];
	delta_i  =knot[l]-knot[l1];
	sum=delta_im1+delta_i;

	alpha[l1]=delta_i*delta_i/(delta_im2+sum);
	beta[l1] =delta_i*(delta_im2+delta_im1)/(delta_im2  + sum)
                 +
		  delta_im1*delta_i / sum;
	gamma[l1]=delta_im1*delta_im1/sum;

	alpha[l1]=alpha[l1]/sum;
	beta[l1] =beta[l1]/sum;
	gamma[l1]=gamma[l1]/sum;


	alpha[l]=0.0; beta[l]=1.0; gamma[l]=0.0;


}

//::::::::::::::
//solve_system.c
//::::::::::::::
 void solve_system(float *up,float *low,float *gamma,int l,float *rhs,float *d)
/*	solve  tridiagonal linear system
	of size (l+1)(l+1) whose LU decompostion has entries up and low,
	and whose right hand side is rhs, and whose original matrix
	had upper diagonal gamma. Solution is d[0],...,d[l+2];
Input: up,low,gamma:  as above.
       l:             size of system: l+1 eqs in l+1 unknowns.
       rhs:           right hand side, i.e, data points with end
                      `tangent Bezier points' in rhs[1] and rhs[l+1].
Output:d:             solution vector.

Note:	Both rhs and d are from 0 to l+2.
*/
	//float up[],low[],gamma[],rhs[],d[];
	//int l;
{
	int i;
	float aux[100];

	d[0] = rhs[0];
	d[1] = rhs[1];

	/* forward substitution:  */
	aux[0]=rhs[1];
	for(i=1; i<=l; i++) aux[i]=rhs[i+1]-low[i]*aux[i-1];
	
	/* backward substitution:  */
	d[l+1]=aux[l]/up[l];
	for(i=l-1; i>0; i--) d[i+1]=(aux[i]-gamma[i]*d[i+2])/up[i];
	d[l+2]=rhs[l+2];
}
		

//::::::::::::::
//subdiv.c
//::::::::::::::


void subdiv(int degree,float *coeff,float *weight,float t,float *bleft,float *bright,float *wleft,float *wright)
/*
	subdivides ratbez curve at parameter value t.
	Output: left and right polygon with respective weights.
	Ordering of right polygon is reversed.
*/
	//float coeff[],weight[],bleft[],
	//	bright[],wleft[],wright[];
	//float t;
	//int degree; 
{
	int r,i;
	float t1,ww1,ww2;
	t1 = 1.0 - t;
	
/* 
	first, obtain right subpolygon from rat de Casteljau
*/

	for (i=0;i <= degree; i++) wright[i] = weight[i];
 
	for (i=0;i <= degree; i++) bright[i] = coeff[i];

	for (r=1; r<= degree; r++)
	for (i=0; i<= degree - r; i++)
	{
		ww1 = wright[i]; ww2= wright[i+1];
		wright[i] = t1*wright[i] + t*wright[i+1];
		bright[i]= ( t1*ww1 * bright[i]
			  +  t*ww2 * bright[i+1] ) / wright[i];
	}	

/* 
	use same as above in order to get left half. Idea:
	reverse ordering; then the above yields left half.
 */

	t = 1.0 - t; t1 = 1.0 - t;
	for (i=0;i <= degree; i++) wleft[degree -i] = weight[i];
 
	for (i=0;i <= degree; i++) bleft[degree-i] = coeff[i];

	for (r=1; r<= degree; r++)
	for (i=0; i<= degree - r; i++)
	{
		ww1 = wleft[i]; ww2= wleft[i+1];
		wleft[i] = t1*wleft[i] + t*wleft[i+1];
		bleft[i]= ( t1*ww1 * bleft[i]
			  +  t*ww2 * bleft[i+1] ) / wleft[i];
	}	

}
//::::::::::::::
//tri_decas.c
//::::::::::::::
/*
 * Routine: tri_decast()
 *
 * Function: Triangular de Casteljau algorithm for an n^th
 *           degree triangular Bezier patch. 
 *           Algorithm is applied once for a given (u,v,w) and works on one
 *           coordinate only. 
 *
 * Input:    bpts[i] ................. Bezier points (of one coordinate)
 *                                     as a linear array (see below). 
 *                                     i=0...tri_num  
 *           tri_num ................. Based on the degree of the patch.
 *                                     (n+1)(n+2)/2
 *           ndeg .................... Degree (n) of the patch. 
 *           u[i] .................... Barycentric coordinates (u,v,w) of 
 *                                     evaluation point.  i=0,2
 *           b[i]..................... A working array with dimension >= 
 *                                     to bpts[].
 *
 * Output:   patch_pt ................ One coordinate of the point on 
 *                                     the patch evaluated at (u,v,w).
 *           b[] ..................... Contents have been changed. 
 *
 * Linear array structure: It is assumed that the usual (i,k,j) structure
 *                         of the Bezier net has been put into a linear
 *                         array in the following manner.
 *                         (E.g., for n=3)
 *                         b_(300) --> bpts[0] (u=1)
 *                         b_(030) --> bpts[6] (v=1)
 *                         b_(003) --> bpts[9] (w=1)
 */



tri_decast(float *bpts, int tri_num, int ndeg, float *u, float *b, float *patch_pt ) 
    //float	bpts[];
    //int         tri_num;
    //int         ndeg;
    //float       u[3];
    //float       b[];
    //float       *patch_pt;
{

    int         i, j, k, l, m;
    int         r;

    /* ------------------------------ */


    /* 
     * To avoid writing over the original control points, 
     * copy into a working array.
     */

    for (i=0; i<tri_num; i++) 
        b[i] = bpts[i]; 


    /*
     * de Casteljau algorithm applied to the linear array. 
     * The element at the "top" of each mini-triangle in the
     * algorithm is over-written at each step.
     */

    for (r=1; r<=ndeg; r++) {

	m = -1;
	for (i=0; i<=ndeg-r; i++) {

	    for (l=0; l<=i; l++) {

		m = m + 1;
		b[m] = u[0]*b[m] 
                     + u[1]*b[m+1+i] 
                     + u[2]*b[m+2+i]; 
	    }
	}
    } 

    *patch_pt = b[0];

}