/* * Do the last subinterval as a special case since no point follows the * last point */   for (i = 0; i < nb+1; i++) {      sxy[2 * nb * (n-2) + i + i] =	h1[i] * x[n-2] + h2[i] * x[n-1] + h3[i] * D1x;      sxy[2 * nb * (n-2) + i + i + 1] =	h1[i] * y[n-2] + h2[i] * y[n-1] + h3[i] * D1y;   }   g2_free(x);}void g2_raspln(int id, int n, double *points, double tn)/* *	FORMAL ARGUMENTS: * *	id			device id *	n			number of data points *	points			data points (x[i],y[i]) *	tn			tension factor [0.0, 2.0] *				0.0  very rounded *				2.0  not rounded at all */{   int m;   double *sxy;		/*	coords of the entire spline curve */   m = (n-1)*nb+1;   sxy = (double *) g2_malloc(m*2*sizeof(double));   g2_c_raspln(n, points, tn, sxy);   g2_poly_line(id, m, sxy);   g2_free(sxy);}void g2_filled_raspln(int id, int n, double *points, double tn)/* *	FORMAL ARGUMENTS: * *	id			device id *	n			number of data points *	points			data points (x[i],y[i]) *	tn			tension factor [0.0, 2.0] *				0.0  very rounded *				2.0  not rounded at all */{   int m;   double *sxy;		/*	coords of the entire spline curve */   m = (n-1)*nb+2;   sxy = (double *) g2_malloc(m*2*sizeof(double));   g2_c_raspln(n, points, tn, sxy);   sxy[(n+n-2) * nb + 2] = points[n+n-2];   sxy[(n+n-2) * nb + 3] = points[1];   g2_filled_polygon(id, m, sxy);   g2_free(sxy);}/* ---- And now for a rather different approach ---- *//* *	FUNCTION g2_c_newton * *	FUNCTIONAL DESCRIPTION: * *	Use Newton's Divided Differences to calculate an interpolation *	polynomial through the specified data points. *	This function is called by *		g2_c_para_3 and *		g2_c_para_5. * *	Dennis Mikkelson	distributed in GPLOT	Jan  5, 1988	F77 *	Tijs Michels		t.michels@vimec.nl	Jun 16, 1999	C * *	FORMAL ARGUMENTS: * *	n	number of entries in c1 and c2, 4 <= n <= MaxPts *		for para_3	(degree 3)	n = 4 *		for para_5	(degree 5)	n = 6 *		for para_i	(degree i)	n = (i + 1) *	c1	double array holding at most MaxPts values giving the *		first  coords of the points to be interpolated *	c2	double array holding at most MaxPts values giving the *		second coords of the points to be interpolated *	o	number of points at which the interpolation *		polynomial is to be evaluated *	xv	double array holding o points at which to *		evaluate the interpolation polynomial *	yv	double array holding upon return the values of the *		interpolation polynomial at the corresponding points in xv * *		yv is the OUTPUT * *	IMPLICIT INPUTS:	NONE *	IMPLICIT OUTPUTS:	NONE *	SIDE EFFECTS:		NONE */#define MaxPts 21#define xstr(s) __str(s)#define __str(s) #s/* * Maximum number of data points allowed * 21 would correspond to a polynomial of degree 20 */void g2_c_newton(int n, const double *c1, const double *c2,		 int o, const double *xv, double *yv){   int i, j;   double p, s;   double ddt[MaxPts][MaxPts];		/* Divided Difference Table */   if (n < 4) {      fputs("g2_c_newton: Error! Less than 4 points passed "	    "to function g2_c_newton\n", stderr);      return;   }   if (n > MaxPts) {      fputs("g2_c_newton: Error! More than " xstr(MaxPts) " points passed "	    "to function g2_c_newton\n", stderr);      return;   }/* First, build the divided difference table */   for (i = 0; i < n; i++)	ddt[i][0] = c2[i];   for (j = 1; j < n; j++) {      for (i = 0; i < n - j; i++)	ddt[i][j] = (ddt[i+1][j-1] - ddt[i][j-1]) / (c1[i+j] - c1[i]);   }/* Next, evaluate the polynomial at the specified points */   for (i = 0; i < o; i++) {      for (p = 1., s = ddt[0][0], j = 1; j < n; j++) {	 p *= xv[i] - c1[j-1];	 s += p * ddt[0][j];      }      yv[i] = s;   }}/* *	FUNCTION: g2_c_para_3 * *	FUNCTIONAL DESCRIPTION: * *	This function draws a piecewise parametric interpolation *	polynomial of degree 3 through the specified data points. *	The effect is similar to that obtained using DISSPLA to *	draw a curve after a call to the DISSPLA routine PARA3. *	The curve is parameterized using an approximation to the *	curve's arc length. The basic interpolation is done *	using function g2_c_newton. * *	Dennis Mikkelson	distributed in GPLOT	Jan  7, 1988	F77 *	Tijs Michels		t.michels@vimec.nl	Jun 17, 1999	C * *	FORMAL ARGUMENTS: * *	n	number of data points through which to draw the curve *	points	double array containing the x and y-coords of the data points * *	IMPLICIT INPUTS:	NONE *	IMPLICIT OUTPUTS:	NONE *	SIDE EFFECTS:		NONE *//* * #undef  nb * #define nb 40 * Number of straight connecting lines of which each polynomial consists. * So between one data point and the next, (nb-1) points are placed. */void g2_c_para_3(int n, const double *points, double *sxy){#define dgr	(3+1)#define nb2	(nb*2)   int i, j;   double x1t, y1t;   double o, step;   double X[nb2];		/* x-coords of the current curve piece */   double Y[nb2];		/* y-coords of the current curve piece */   double t[dgr];		/* data point parameter values */   double Xpts[dgr];		/* x-coords data point subsection */   double Ypts[dgr];		/* y-coords data point subsection */   double s[nb2];		/* parameter values at which to interpolate */   /* Do first TWO subintervals first */   g2_split(dgr, points, Xpts, Ypts);   t[0] = 0.;   for (i = 1; i < dgr; i++) {      x1t = Xpts[i] - Xpts[i-1];      y1t = Ypts[i] - Ypts[i-1];      t[i] = t[i-1] + sqrt(x1t * x1t + y1t * y1t);   }   step = t[2] / nb2;   for (i = 0; i < nb2; i++)	s[i] = i * step;   g2_c_newton(dgr, t, Xpts, nb2, s, X);   g2_c_newton(dgr, t, Ypts, nb2, s, Y);   for (i = 0; i < nb2; i++) {      sxy[i+i]   = X[i];      sxy[i+i+1] = Y[i];   }   /* Next, do later central subintervals */   for (j = 1; j < n - dgr + 1; j++) {      g2_split(dgr, points + j + j, Xpts, Ypts);      for (i = 1; i < dgr; i++) {	 x1t = Xpts[i] - Xpts[i-1];	 y1t = Ypts[i] - Ypts[i-1];	 t[i] = t[i-1] + sqrt(x1t * x1t + y1t * y1t);      }      o = t[1]; /* look up once */      step = (t[2] - o) / nb;      for (i = 0; i < nb; i++)	s[i] = i * step + o;      g2_c_newton(dgr, t, Xpts, nb, s, X);      g2_c_newton(dgr, t, Ypts, nb, s, Y);      for (i = 0; i < nb; i++) {	 sxy[(j + 1) * nb2 + i + i]     = X[i];	 sxy[(j + 1) * nb2 + i + i + 1] = Y[i];      }   }   /* Now do last subinterval */   o = t[2];   step = (t[3] - o) / nb;   for (i = 0; i < nb; i++)	s[i] = i * step + o;   g2_c_newton(dgr, t, Xpts, nb, s, X);   g2_c_newton(dgr, t, Ypts, nb, s, Y);   for (i = 0; i < nb; i++) {      sxy[(n - dgr + 2) * nb2 + i + i]     = X[i];      sxy[(n - dgr + 2) * nb2 + i + i + 1] = Y[i];   }   sxy[(n - 1) * nb2]     = points[n+n-2];   sxy[(n - 1) * nb2 + 1] = points[n+n-1];}/* *	FORMAL ARGUMENTS: * *	id			device id *	n			number of data points *	points			data points (x[i],y[i]) */void g2_para_3(int id, int n, double *points){   int m;   double *sxy;		/*	coords of the entire spline curve */   m = (n-1)*nb+1;   sxy = (double *) g2_malloc(m*2*sizeof(double));   g2_c_para_3(n, points, sxy);   g2_poly_line(id, m, sxy);   g2_free(sxy);}/* *	FORMAL ARGUMENTS: * *	id			device id *	n			number of data points *	points			data points (x[i],y[i]) */void g2_filled_para_3(int id, int n, double *points){   int m;   double *sxy;		/*	coords of the entire spline curve */   m = (n-1)*nb+2;   sxy = (double *) g2_malloc(m*2*sizeof(double));   g2_c_para_3(n, points, sxy);   sxy[m+m-2] = points[n+n-2];   sxy[m+m-1] = points[1];   g2_filled_polygon(id, m, sxy);   g2_free(sxy);}/* *	FUNCTION: g2_c_para_5 * *	As g2_c_para_3, but now plot a polynomial of degree 5 *//* * #undef  nb * #define nb 40 * Number of straight connecting lines of which each polynomial consists. * So between one data point and the next, (nb-1) points are placed. */void g2_c_para_5(int n, const double *points, double *sxy){#undef	dgr#define dgr	(5+1)#define nb3	(nb*3)   int i, j;   double x1t, y1t;   double o, step;   double X[nb3];		/* x-coords of the current curve piece */   double Y[nb3];		/* y-coords of the current curve piece */   double t[dgr];		/* data point parameter values */   double Xpts[dgr];		/* x-coords data point subsection */   double Ypts[dgr];		/* y-coords data point subsection */   double s[nb3];		/* parameter values at which to interpolate */   /* Do first THREE subintervals first */   g2_split(dgr, points, Xpts, Ypts);   t[0] = 0.;   for (i = 1; i < dgr; i++) {      x1t = Xpts[i] - Xpts[i-1];      y1t = Ypts[i] - Ypts[i-1];      t[i] = t[i-1] + sqrt(x1t * x1t + y1t * y1t);   }   step = t[3] / nb3;   for (i = 0; i < nb3; i++)	s[i] = i * step;   g2_c_newton(dgr, t, Xpts, nb3, s, X);   g2_c_newton(dgr, t, Ypts, nb3, s, Y);   for (i = 0; i < nb3; i++) {      sxy[i+i]   = X[i];      sxy[i+i+1] = Y[i];   }   /* Next, do later central subintervals */   for (j = 1; j < n - dgr + 1; j++) {      g2_split(dgr, points + j + j, Xpts, Ypts);      for (i = 1; i < dgr; i++) {	 x1t = Xpts[i] - Xpts[i-1];	 y1t = Ypts[i] - Ypts[i-1];	 t[i] = t[i-1] + sqrt(x1t * x1t + y1t * y1t);      }      o = t[2]; /* look up once */      step = (t[3] - o) / nb;      for (i = 0; i < nb; i++)	s[i] = i * step + o;      g2_c_newton(dgr, t, Xpts, nb, s, X);      g2_c_newton(dgr, t, Ypts, nb, s, Y);      for (i = 0; i < nb; i++) {	 sxy[(j + 2) * nb2 + i + i]     = X[i];	 sxy[(j + 2) * nb2 + i + i + 1] = Y[i];      }   }   /* Now do last TWO subinterval */   o = t[3];   step = (t[5] - o) / nb2;   for (i = 0; i < nb2; i++)	s[i] = i * step + o;   g2_c_newton(dgr, t, Xpts, nb2, s, X);   g2_c_newton(dgr, t, Ypts, nb2, s, Y);   for (i = 0; i < nb2; i++) {      sxy[(n - dgr + 3) * nb2 + i + i]     = X[i];      sxy[(n - dgr + 3) * nb2 + i + i + 1] = Y[i];   }   sxy[(n - 1) * nb2]     = points[n+n-2];   sxy[(n - 1) * nb2 + 1] = points[n+n-1];}/* *	FORMAL ARGUMENTS: * *	id			device id *	n			number of data points *	points			data points (x[i],y[i]) */void g2_para_5(int id, int n, double *points){   int m;   double *sxy;		/*	coords of the entire spline curve */   m = (n-1)*nb+1;   sxy = (double *) g2_malloc(m*2*sizeof(double));   g2_c_para_5(n, points, sxy);   g2_poly_line(id, m, sxy);   g2_free(sxy);}/* *	FORMAL ARGUMENTS: * *	id			device id *	n			number of data points *	points			data points (x[i],y[i]) */void g2_filled_para_5(int id, int n, double *points){   int m;   double *sxy;		/*	coords of the entire spline curve */   m = (n-1)*nb+2;   sxy = (double *) g2_malloc(m*2*sizeof(double));   g2_c_para_5(n, points, sxy);   sxy[m+m-2] = points[n+n-2];   sxy[m+m-1] = points[1];   g2_filled_polygon(id, m, sxy);   g2_free(sxy);}