/* ARC LENGTH
*
*/

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

/* #define DEBUG */

/* ------------------------------------------------------------------------
STRUCTURES
*/

// the structure to hold entries of the table consisting of
// parameter value (u) and estimated length (length) 
typedef struct table_entry_structure {
	double u,length;
} table_entry_td;

// the structure to hold an interval of the curve, defined by 
// starting and ending parameter values and the estimated 
// length (to be filled in by the adaptive integration
// procedure
typedef struct interval_structure {
	double		u1,u2;
	double		length;
}  interval_td;

// coefficients for a 2D cubic curve
typedef struct cubic_curve_structure {
	double		ax,bx,cx,dx;
	double		ay,by,cy,dy;
} cubic_curve_td;

// polynomial function structure; a quadric function is generated
// from a cubic curve during the arclength computation
typedef struct polynomial_structure {
	double		*coeff;
	int		degree;
} polynomial_td;

float adaptive_integration(cubic_curve_td *curve, double u1, double u2, double tolerance);
double subdivide(interval_td *full_interval, polynomial_td *func, 
		double total_length, double tolerance);
void add_table_entry(double u, double length);
double integrate_func(polynomial_td *func,interval_td *interval);
double evaluate_polynomial(polynomial_td *poly, double u);

/* ------------------------------------------------------------------------
ADAPTIVE INTEGRATION
this is the high level call used whenver a curve's length is to be computed
 */
float adaptive_integration(cubic_curve_td *curve, double u1, double u2, double tolerance)
{
  double subdivide();
  polynomial_td	func;
  interval_td full_interval;
  double  total_length;
  double  integrate_func();
  double  temp;

  func.degree = 4;
  func.coeff = (double *)malloc(sizeof(double)*5);
  func.coeff[4] = 9*(curve->ax*curve->ax + curve->ay*curve->ay);
  func.coeff[3] = 12*(curve->ax*curve->bx + curve->ay*curve->by);
  func.coeff[2] = (6*(curve->ax*curve->cx + curve->ay*curve->cy) +
		  4*(curve->bx*curve->bx + curve->by*curve->by)   );
  func.coeff[1] = 4*(curve->bx*curve->cx + curve->by*curve->cy);
  func.coeff[0] = curve->cx*curve->cx + curve->cy*curve->cy;

  full_interval.u1 = u1; full_interval.u2 = u2;
  temp = integrate_func(&func,&full_interval);
  printf("\nInitial guess = %lf; %lf:%lf",temp,u1,u2);
  full_interval.length = temp;
  total_length = subdivide(&full_interval,&func,0.0,tolerance);
  printf("\n total length = %lf\n",total_length);
  
  return total_length;
}

/* ------------------------------------------------------------------------
SUBDIVIDE
'total_length' is the length of the curve up to, but not including, the 'full_interval' 
if the difference between the interval and the sum of its halfs is less than 'tolerance', 
     stop the recursive subdivision
'func' is a polynomial function
*/
double subdivide(interval_td *full_interval, polynomial_td *func, 
		double total_length, double tolerance)
{

  interval_td	left_interval, right_interval;
  double  left_length,right_length;
  double  midu;
  double  subdivide();
  double integrate_func();
  double  temp;
  void add_table_entry();

  midu = (full_interval->u1+full_interval->u2)/2; 
  left_interval.u1 = full_interval->u1; left_interval.u2 = midu;
  right_interval.u1 = midu; right_interval.u2 = full_interval->u2;

  left_length = integrate_func(func, &left_interval);
  right_length = integrate_func(func,&right_interval);

  temp = fabs(full_interval->length - (left_length+right_length));
  if (temp > tolerance) {
    left_interval.length = left_length;
    right_interval.length = right_length;
    total_length = subdivide(&left_interval, func, total_length, tolerance);
    total_length = subdivide(&right_interval, func, total_length, tolerance);
    return(total_length);
  }
  else {
    total_length = total_length + left_length;
    add_table_entry(midu,total_length);
    total_length = total_length + right_length;
    add_table_entry(full_interval->u2,total_length);
    return(total_length);
  }
}

/* ------------------------------------------------------------------------
ADD TABLE ENTRY
adds an entry of the form (parametric value, length) to the table being constructed
*/
void add_table_entry(double u, double length)
{
  /* add entry of (u, length) */
  printf("\ntable entry:  u: %lf, length: %lf",u,length);
}

/* ------------------------------------------------------------------------
INTEGRATE FUNCTION
use gaussian quadrature to integrate square root of given function in thegiven interval
*/
double integrate_func(polynomial_td *func,interval_td *interval)
{
  double x[5]={.1488743389,.4333953941,.6794095682,.8650633666, .9739065285};
  double w[5]={.2966242247,.2692667193,.2190863625,.1494513491,.0666713443};
  double		length, midu, dx, diff;
  int		i;
  double  evaluate_polynomial();
  double  u1,u2;

  u1 = interval->u1;
  u2 = interval->u2;

  midu = (u1+u2)/2.0;
  diff = (u2-u1)/2.0;
  length = 0.0;
  for (i=0; i<5; i++) {
    dx = diff*x[i];
    length += w[i]*(sqrt(evaluate_polynomial(func,midu+dx)) + 
		    sqrt(evaluate_polynomial(func,midu-dx)));
  }
  length *= diff;

  return (length);
}

/* ------------------------------------------------------------------------
EVALUATE POLYNOMIAL
evaluate a polynomial
 */
double evaluate_polynomial(polynomial_td *poly, double u)
{
  double  w;
  int    i;
  double  value;

  value = 0.0;
  w = 1.0;
  for (i=0; i<=poly->degree; i++) {
    value += poly->coeff[i]*w;
    w *= u;
  }
  return value;
}