--- /dev/null
+/*****************************************************************************
+** Copyright (C) 1998-2001 Ljubomir Milanovic & Horst Wagner
+** This file is part of the g2 library
+**
+** This library is free software; you can redistribute it and/or
+** modify it under the terms of the GNU Lesser General Public
+** License as published by the Free Software Foundation; either
+** version 2.1 of the License, or (at your option) any later version.
+**
+** This library is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+** Lesser General Public License for more details.
+**
+** You should have received a copy of the GNU Lesser General Public
+** License along with this library; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+******************************************************************************/
+/*
+ * g2_splines.c
+ * Tijs Michels
+ * tijs@vimec.nl
+ * 06/16/99
+ */
+
+#include <math.h>
+#include <stdio.h>
+#include "g2.h"
+#include "g2_util.h"
+
+static void g2_split(int n, const double *points, double *x, double *y);
+static void g2_c_spline(int n, const double *points, int m, double *sxy);
+static void g2_c_b_spline(int n, const double *points, int m, double *sxy);
+static void g2_c_raspln(int n, const double *points, double tn, double *sxy);
+static void g2_c_newton(int n, const double *c1, const double *c2, int o, const double *xv, double *yv);
+static void g2_c_para_3(int n, const double *points, double *sxy);
+static void g2_c_para_5(int n, const double *points, double *sxy);
+
+void g2_split(int n, const double *points, double *x, double *y)
+{
+ int i;
+ for (i = 0; i < n; i++) {
+ x[i] = points[i+i];
+ y[i] = points[i+i+1];
+ }
+}
+
+#define eps 1.e-12
+
+void g2_c_spline(int n, const double *points, int m, double *sxy)
+
+/*
+ * FUNCTIONAL DESCRIPTION:
+ *
+ * Compute a curve of m points (sx[j],sy[j])
+ * -- j being a positive integer < m --
+ * passing through the n data points (x[i],y[i])
+ * -- i being a positive integer < n --
+ * supplied by the user.
+ * The procedure to determine sy[j] involves
+ * Young's method of successive over-relaxation.
+ *
+ * FORMAL ARGUMENTS:
+ *
+ * n number of data points
+ * points data points (x[i],y[i])
+ * m number of interpolated points; m = (n-1)*o+1
+ * for o curve points for every data point
+ * sxy interpolated points (sx[j],sy[j])
+ *
+ * IMPLICIT INPUTS: NONE
+ * IMPLICIT OUTPUTS: NONE
+ * SIDE EFFECTS: NONE
+ *
+ * REFERENCES:
+ *
+ * 1. Ralston and Wilf, Mathematical Methods for Digital Computers,
+ * Vol. II, John Wiley and Sons, New York 1967, pp. 156-158.
+ * 2. Greville, T.N.E., Ed., Proceedings of An Advanced Seminar
+ * Conducted by the Mathematics Research Center, U.S. Army,
+ * University of Wisconsin, Madison. October 7-9, 1968. Theory
+ * and Applications of Spline Functions, Academic Press,
+ * New York / London 1969, pp. 156-167.
+ *
+ * AUTHORS:
+ *
+ * Josef Heinen 04/06/88 <J.Heinen@KFA-Juelich.de>
+ * Tijs Michels 06/16/99 <t.michels@vimec.nl>
+ */
+
+{
+ int i, j;
+ double *x, *y, *g, *h;
+ double k, u, delta_g;
+
+ if (n < 3) {
+ fputs("\nERROR calling function \"g2_c_spline\":\n"
+ "number of data points input should be at least three\n", stderr);
+ return;
+ }
+ if ((m-1)%(n-1)) {
+ fputs("\nWARNING from function \"g2_c_spline\":\n"
+ "number of curve points output for every data point input "
+ "is not an integer\n", stderr);
+ }
+
+ x = (double *) g2_malloc(n*4*sizeof(double));
+ y = x + n;
+ g = y + n;
+ h = g + n; /* for the constant copy of g */
+ g2_split(n, points, x, y);
+
+ n--; /* last value index */
+ k = x[0]; /* look up once */
+ u = (x[n] - k) / (m - 1); /* calculate step outside loop */
+ for (j = 0; j < m; j++) sxy[j+j] = j * u + k; /* x-coordinates */
+
+ for (i = 1; i < n; i++) {
+ g[i] = 2. * ((y[i+1] - y[i]) / (x[i+1] - x[i]) -
+ (y[i] - y[i-1]) / (x[i] - x[i-1]))
+ / (x[i+1] - x[i-1]); /* whereas g[i] will later be changed repeatedly */
+ h[i] = 1.5 * g[i]; /* copy h[i] of g[i] will remain constant */
+ }
+
+ k = 0.;
+
+ do {
+ for (u = 0., i = 1; i < n; i++) {
+ delta_g = .5 * (x[i] - x[i-1]) / (x[i+1] - x[i-1]);
+ delta_g = (h[i] -
+ g[i] -
+ g[i-1] * delta_g - /* 8. - 4 * sqrt(3.) */
+ g[i+1] * (.5 - delta_g)) * 1.0717967697244907832;
+ g[i] += delta_g;
+
+ if (fabs(delta_g) > u) u = fabs(delta_g);
+ } /* On loop termination u holds the largest delta_g. */
+
+ if (k == 0.) k = u * eps;
+ /* Only executed once, at the end of pass one. So k preserves
+ * the largest delta_g of pass one, multiplied by eps.
+ */
+ } while (u > k);
+
+ m += m, i = 1, j = 0;
+ do {
+ u = sxy[j++]; /* x-coordinate */
+
+ while (x[i] < u) i++;
+
+ if (--i > n) i = n;
+
+ k = (u - x[i]) / (x[i+1] - x[i]); /* calculate outside loop */
+ sxy[j++] = y[i] +
+ (y[i+1] - y[i]) * k +
+ (u - x[i]) * (u - x[i+1]) *
+ ((2. - k) * g[i] +
+ (1. + k) * g[i+1]) / 6.; /* y-coordinate */
+ } while (j < m);
+ g2_free(x);
+}
+
+void g2_spline(int id, int n, double *points, int o)
+
+/*
+ * FORMAL ARGUMENTS:
+ *
+ * id device id
+ * n number of data points
+ * points data points (x[i],y[i])
+ * o number of interpolated points per data point
+ *
+ * Given an array of n data points {x[1], y[1], ... x[n], y[n]} plot a
+ * spline curve on device id with o interpolated points per data point.
+ * So the larger o, the more fluent the curve.
+ */
+
+{
+ int m;
+ double *sxy;
+
+ m = (n-1)*o+1;
+ sxy = (double*)g2_malloc(m*2*sizeof(double));
+
+ g2_c_spline(n, points, m, sxy);
+ g2_poly_line(id, m, sxy);
+
+ g2_free(sxy);
+}
+
+void g2_filled_spline(int id, int n, double *points, int o)
+
+/*
+ * FORMAL ARGUMENTS:
+ *
+ * id device id
+ * n number of data points
+ * points data points (x[i],y[i])
+ * o number of interpolated points per data point
+ */
+
+{
+ int m;
+ double *sxy;
+
+ m = (n-1)*o+1;
+ sxy = (double*)g2_malloc((m+1)*2*sizeof(double));
+
+ g2_c_spline(n, points, m, sxy);
+ sxy[m+m] = points[n+n-2];
+ sxy[m+m+1] = points[1];
+ g2_filled_polygon(id, m+1, sxy);
+ g2_free(sxy);
+}
+
+void g2_c_b_spline(int n, const double *points, int m, double *sxy)
+
+/*
+ * g2_c_b_spline takes n input points. It uses parameter t
+ * to compute sx(t) and sy(t) respectively
+ */
+
+{
+ int i, j;
+ double *x, *y;
+ double t, bl1, bl2, bl3, bl4;
+ double interval, xi_3, yi_3, xi, yi;
+
+ if (n < 3) {
+ fputs("\nERROR calling function \"g2_c_b_spline\":\n"
+ "number of data points input should be at least three\n", stderr);
+ return;
+ }
+ x = (double *) g2_malloc(n*2*sizeof(double));
+ y = x + n;
+ g2_split(n, points, x, y);
+
+ m--; /* last value index */
+ n--; /* last value index */
+ interval = (double)n / m;
+
+ for (m += m, i = 2, j = 0; i <= n+1; i++) {
+ if (i == 2) {
+ xi_3 = 2 * x[0] - x[1];
+ yi_3 = 2 * y[0] - y[1];
+ } else {
+ xi_3 = x[i-3];
+ yi_3 = y[i-3];
+ }
+ if (i == n+1) {
+ xi = 2 * x[n] - x[n-1];
+ yi = 2 * y[n] - y[n-1];
+ } else {
+ xi = x[i];
+ yi = y[i];
+ }
+
+ t = fmod(j * interval, 1.);
+
+ while (t < 1. && j < m) {
+ bl1 = (1. - t);
+ bl2 = t * t; /* t^2 */
+ bl4 = t * bl2; /* t^3 */
+ bl3 = bl4 - bl2;
+
+ bl1 = bl1 * bl1 * bl1;
+ bl2 = 3. * (bl3 - bl2) + 4.;
+ bl3 = 3. * ( t - bl3) + 1.;
+
+ sxy[j++] = (bl1 * xi_3 + bl2 * x[i-2] + bl3 * x[i-1] + bl4 * xi) / 6.; /* x-coordinate */
+ sxy[j++] = (bl1 * yi_3 + bl2 * y[i-2] + bl3 * y[i-1] + bl4 * yi) / 6.; /* y-coordinate */
+
+ t += interval;
+ }
+ }
+ sxy[m] = x[n];
+ sxy[m+1] = y[n];
+ g2_free(x);
+}
+
+void g2_b_spline(int id, int n, double *points, int o)
+
+/*
+ * FORMAL ARGUMENTS:
+ *
+ * id device id
+ * n number of data points
+ * points data points (x[i],y[i])
+ * o number of interpolated points per data point
+ */
+
+{
+ int m;
+ double *sxy;
+
+ m = (n-1)*o+1;
+ sxy = (double*)g2_malloc(m*2*sizeof(double));
+
+ g2_c_b_spline(n, points, m, sxy);
+ g2_poly_line(id, m, sxy);
+
+ g2_free(sxy);
+}
+
+void g2_filled_b_spline(int id, int n, double *points, int o)
+
+/*
+ * FORMAL ARGUMENTS:
+ *
+ * id device id
+ * n number of data points
+ * points data points (x[i],y[i])
+ * o number of interpolated points per data point
+ */
+
+{
+ int m;
+ double *sxy;
+
+ m = (n-1)*o+1;
+ sxy = (double*)g2_malloc((m+1)*2*sizeof(double));
+
+ g2_c_b_spline(n, points, m, sxy);
+ sxy[m+m] = points[n+n-2];
+ sxy[m+m+1] = points[1];
+ g2_filled_polygon(id, m+1, sxy);
+
+ g2_free(sxy);
+}
+
+/*
+ * FUNCTION g2_c_raspln
+ *
+ * FUNCTIONAL DESCRIPTION:
+ *
+ * This function draws a piecewise cubic polynomial through
+ * the specified data points. The (n-1) cubic polynomials are
+ * basically parametric cubic Hermite polynomials through the
+ * n specified data points with tangent values at the data
+ * points determined by a weighted average of the slopes of
+ * the secant lines. A tension parameter "tn" is provided to
+ * adjust the length of the tangent vector at the data points.
+ * This allows the "roundness" of the curve to be adjusted.
+ * For further information and references on this technique see:
+ *
+ * D. Kochanek and R. Bartels, Interpolating Splines With Local
+ * Tension, Continuity and Bias Control, Computer Graphics,
+ * 18(1984)3.
+ *
+ * AUTHORS:
+ *
+ * Dennis Mikkelson distributed in GPLOT Jan 7, 1988 F77
+ * Tijs Michels t.michels@vimec.nl Jun 7, 1999 C
+ *
+ * FORMAL ARGUMENTS:
+ *
+ * n number of data points, n > 2
+ * points double array holding the x and y-coords of the data points
+ * tn double parameter in [0.0, 2.0]. When tn = 0.0,
+ * the curve through the data points is very rounded.
+ * As tn increases the curve is gradually pulled tighter.
+ * When tn = 2.0, the curve is essentially a polyline
+ * through the given data points.
+ * sxy double array holding the coords of the spline curve
+ *
+ * IMPLICIT INPUTS: NONE
+ * IMPLICIT OUTPUTS: NONE
+ * SIDE EFFECTS: NONE
+ */
+
+#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_raspln(int n, const double *points, double tn, double *sxy)
+{
+ int i, j;
+ double *x, *y;
+ double bias, tnFactor, tangentL1, tangentL2;
+ double D1x, D1y, D2x, D2y, t1x, t1y, t2x, t2y;
+ double h1[nb+1]; /* Values of the Hermite basis functions */
+ double h2[nb+1]; /* at nb+1 evenly spaced points in [0,1] */
+ double h3[nb+1];
+ double h4[nb+1];
+
+ x = (double *) g2_malloc(n*2*sizeof(double));
+ y = x + n;
+ g2_split(n, points, x, y);
+
+/*
+ * First, store the values of the Hermite basis functions in a table h[ ]
+ * so no time is wasted recalculating them
+ */
+ for (i = 0; i < nb+1; i++) {
+ double t, tt, ttt;
+ t = (double) i / nb;
+ tt = t * t;
+ ttt = t * tt;
+ h1[i] = 2. * ttt - 3. * tt + 1.;
+ h2[i] = -2. * ttt + 3. * tt;
+ h3[i] = ttt - 2. * tt + t;
+ h4[i] = ttt - tt;
+ }
+
+/*
+ * Set local tnFactor based on input parameter tn
+ */
+ if (tn <= 0.) {
+ tnFactor = 2.;
+ fputs("g2_c_raspln: Using Tension Factor 0.0: very rounded", stderr);
+ }
+ else if (tn >= 2.) {
+ tnFactor = 0.;
+ fputs("g2_c_raspln: Using Tension Factor 2.0: not rounded at all", stderr);
+ }
+ else tnFactor = 2. - tn;
+
+ D1x = D1y = 0.; /* first point has no preceding point */
+ for (j = 0; j < n - 2; j++) {
+ t1x = x[j+1] - x[j];
+ t1y = y[j+1] - y[j];
+ t2x = x[j+2] - x[j+1];
+ t2y = y[j+2] - y[j+1];
+ tangentL1 = t1x * t1x + t1y * t1y;
+ tangentL2 = t2x * t2x + t2y * t2y;
+ if (tangentL1 + tangentL2 == 0) bias = .5;
+ else bias = tangentL2 / (tangentL1 + tangentL2);
+ D2x = tnFactor * (bias * t1x + (1 - bias) * t2x);
+ D2y = tnFactor * (bias * t1y + (1 - bias) * t2y);
+ for (i = 0; i < nb; i++) {
+ sxy[2 * nb * j + i + i] =
+ h1[i] * x[j] + h2[i] * x[j+1] + h3[i] * D1x + h4[i] * D2x;
+ sxy[2 * nb * j + i + i + 1] =
+ h1[i] * y[j] + h2[i] * y[j+1] + h3[i] * D1y + h4[i] * D2y;
+ }
+ D1x = D2x; /* store as preceding point in */
+ D1y = D2y; /* the next pass */
+ }
+
+/*
+ * 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);
+}
+
git://source.jalview.org / jabaws.git / blobdiff