g2_splines.c

ViennaRNA-1.6.1

/*
 * 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		[e-mail]	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		[e-mail]	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);
}