binaries/src/ViennaRNA/RNAforester/g2-0.70/src/g2_splines.c

   1 /*****************************************************************************
   2 **  Copyright (C) 1998-2001  Ljubomir Milanovic & Horst Wagner
   3 **  This file is part of the g2 library
   4 **
   5 **  This library is free software; you can redistribute it and/or
   6 **  modify it under the terms of the GNU Lesser General Public
   7 **  License as published by the Free Software Foundation; either
   8 **  version 2.1 of the License, or (at your option) any later version.
   9 **
  10 **  This library is distributed in the hope that it will be useful,
  11 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 **  Lesser General Public License for more details.
  14 **
  15 **  You should have received a copy of the GNU Lesser General Public
  16 **  License along with this library; if not, write to the Free Software
  17 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  18 ******************************************************************************/
  19 /*
  20  * g2_splines.c
  21  * Tijs Michels
  22  * tijs@vimec.nl
  23  * 06/16/99
  24  */
  25
  26 #include <math.h>
  27 #include <stdio.h>
  28 #include "g2.h"
  29 #include "g2_util.h"
  30
  31 static void g2_split(int n, const double *points, double *x, double *y);
  32 static void g2_c_spline(int n, const double *points, int m, double *sxy);
  33 static void g2_c_b_spline(int n, const double *points, int m, double *sxy);
  34 static void g2_c_raspln(int n, const double *points, double tn, double *sxy);
  35 static void g2_c_newton(int n, const double *c1, const double *c2, int o, const double *xv, double *yv);
  36 static void g2_c_para_3(int n, const double *points, double *sxy);
  37 static void g2_c_para_5(int n, const double *points, double *sxy);
  38
  39 void g2_split(int n, const double *points, double *x, double *y)
  40 {
  41    int i;
  42    for (i = 0; i < n; i++) {
  43       x[i] = points[i+i];
  44       y[i] = points[i+i+1];
  45    }
  46 }
  47
  48 #define eps 1.e-12
  49
  50 void g2_c_spline(int n, const double *points, int m, double *sxy)
  51
  52 /*
  53  *      FUNCTIONAL DESCRIPTION:
  54  *
  55  *      Compute a curve of m points (sx[j],sy[j])
  56  *      -- j being a positive integer < m --
  57  *      passing through the n data points (x[i],y[i])
  58  *      -- i being a positive integer < n --
  59  *      supplied by the user.
  60  *      The procedure to determine sy[j] involves
  61  *      Young's method of successive over-relaxation.
  62  *
  63  *      FORMAL ARGUMENTS:
  64  *
  65  *      n                       number of data points
  66  *      points                  data points (x[i],y[i])
  67  *      m                       number of interpolated points; m = (n-1)*o+1
  68  *                              for o curve points for every data point
  69  *      sxy                     interpolated points (sx[j],sy[j])
  70  *
  71  *      IMPLICIT INPUTS:        NONE
  72  *      IMPLICIT OUTPUTS:       NONE
  73  *      SIDE EFFECTS:           NONE
  74  *
  75  *      REFERENCES:
  76  *
  77  *      1. Ralston and Wilf, Mathematical Methods for Digital Computers,
  78  *         Vol. II, John Wiley and Sons, New York 1967, pp. 156-158.
  79  *      2. Greville, T.N.E., Ed., Proceedings of An Advanced Seminar
  80  *         Conducted by the Mathematics Research Center, U.S. Army,
  81  *         University of Wisconsin, Madison. October 7-9, 1968. Theory
  82  *         and Applications of Spline Functions, Academic Press,
  83  *         New York / London 1969, pp. 156-167.
  84  *
  85  *      AUTHORS:
  86  *
  87  *      Josef Heinen    04/06/88        <J.Heinen@KFA-Juelich.de>
  88  *      Tijs Michels    06/16/99        <t.michels@vimec.nl>
  89  */
  90
  91 {
  92    int i, j;
  93    double *x, *y, *g, *h;
  94    double k, u, delta_g;
  95
  96    if (n < 3) {
  97       fputs("\nERROR calling function \"g2_c_spline\":\n"
  98             "number of data points input should be at least three\n", stderr);
  99       return;
 100    }
 101    if ((m-1)%(n-1)) {
 102       fputs("\nWARNING from function \"g2_c_spline\":\n"
 103             "number of curve points output for every data point input "
 104             "is not an integer\n", stderr);
 105    }
 106
 107    x = (double *) g2_malloc(n*4*sizeof(double));
 108    y = x + n;
 109    g = y + n;
 110    h = g + n; /* for the constant copy of g */
 111    g2_split(n, points, x, y);
 112
 113    n--; /* last value index */
 114    k =  x[0]; /* look up once */
 115    u = (x[n] - k) / (m - 1); /* calculate step outside loop */
 116    for (j = 0; j < m; j++)      sxy[j+j] = j * u + k; /* x-coordinates */
 117
 118    for (i = 1; i < n; i++) {
 119       g[i] = 2. * ((y[i+1] - y[i]) / (x[i+1] - x[i]) -
 120                    (y[i] - y[i-1]) / (x[i] - x[i-1]))
 121         / (x[i+1] - x[i-1]); /* whereas g[i] will later be changed repeatedly */
 122       h[i] = 1.5 * g[i];     /* copy h[i] of g[i] will remain constant */
 123    }
 124
 125    k = 0.;
 126
 127    do {
 128       for (u = 0., i = 1; i < n; i++) {
 129          delta_g = .5 * (x[i] - x[i-1]) / (x[i+1] - x[i-1]);
 130          delta_g = (h[i] -
 131                     g[i] -
 132                     g[i-1] * delta_g -      /* 8. - 4 * sqrt(3.) */
 133                     g[i+1] * (.5 - delta_g)) * 1.0717967697244907832;
 134          g[i] += delta_g;
 135
 136          if (fabs(delta_g) > u) u = fabs(delta_g);
 137       } /* On loop termination u holds the largest delta_g. */
 138
 139       if (k == 0.)      k = u * eps;
 140         /* Only executed once, at the end of pass one. So k preserves
 141          * the largest delta_g of pass one, multiplied by eps.
 142          */
 143    } while (u > k);
 144
 145    m += m, i = 1, j = 0;
 146    do {
 147       u = sxy[j++]; /* x-coordinate */
 148
 149       while (x[i] < u)  i++;
 150
 151       if (--i > n)      i = n;
 152
 153       k = (u - x[i]) / (x[i+1] - x[i]); /* calculate outside loop */
 154       sxy[j++] = y[i] +
 155         (y[i+1] - y[i]) * k +
 156         (u - x[i]) * (u - x[i+1]) *
 157         ((2. - k) * g[i] +
 158          (1. + k) * g[i+1]) / 6.; /* y-coordinate */
 159    } while (j < m);
 160    g2_free(x);
 161 }
 162
 163 void g2_spline(int id, int n, double *points, int o)
 164
 165 /*
 166  *      FORMAL ARGUMENTS:
 167  *
 168  *      id                      device id
 169  *      n                       number of data points
 170  *      points                  data points (x[i],y[i])
 171  *      o                       number of interpolated points per data point
 172  *
 173  *      Given an array of n data points {x[1], y[1], ... x[n], y[n]} plot a
 174  *      spline curve on device id with o interpolated points per data point.
 175  *      So the larger o, the more fluent the curve.
 176  */
 177
 178 {
 179    int m;
 180    double *sxy;
 181
 182    m = (n-1)*o+1;
 183    sxy = (double*)g2_malloc(m*2*sizeof(double));
 184
 185    g2_c_spline(n, points, m, sxy);
 186    g2_poly_line(id, m, sxy);
 187
 188    g2_free(sxy);
 189 }
 190
 191 void g2_filled_spline(int id, int n, double *points, int o)
 192
 193 /*
 194  *      FORMAL ARGUMENTS:
 195  *
 196  *      id                      device id
 197  *      n                       number of data points
 198  *      points                  data points (x[i],y[i])
 199  *      o                       number of interpolated points per data point
 200  */
 201
 202 {
 203    int m;
 204    double *sxy;
 205
 206    m = (n-1)*o+1;
 207    sxy = (double*)g2_malloc((m+1)*2*sizeof(double));
 208
 209    g2_c_spline(n, points, m, sxy);
 210    sxy[m+m] = points[n+n-2];
 211    sxy[m+m+1] = points[1];
 212    g2_filled_polygon(id, m+1, sxy);
 213    g2_free(sxy);
 214 }
 215
 216 void g2_c_b_spline(int n, const double *points, int m, double *sxy)
 217
 218 /*
 219  * g2_c_b_spline takes n input points. It uses parameter t
 220  * to compute sx(t) and sy(t) respectively
 221  */
 222
 223 {
 224    int i, j;
 225    double *x, *y;
 226    double t, bl1, bl2, bl3, bl4;
 227    double interval, xi_3, yi_3, xi, yi;
 228
 229    if (n < 3) {
 230       fputs("\nERROR calling function \"g2_c_b_spline\":\n"
 231             "number of data points input should be at least three\n", stderr);
 232       return;
 233    }
 234    x = (double *) g2_malloc(n*2*sizeof(double));
 235    y = x + n;
 236    g2_split(n, points, x, y);
 237
 238    m--; /* last value index */
 239    n--; /* last value index */
 240    interval = (double)n / m;
 241
 242    for (m += m, i = 2, j = 0; i <= n+1; i++) {
 243       if (i == 2) {
 244          xi_3 = 2 * x[0] - x[1];
 245          yi_3 = 2 * y[0] - y[1];
 246       } else {
 247          xi_3 = x[i-3];
 248          yi_3 = y[i-3];
 249       }
 250       if (i == n+1) {
 251          xi = 2 * x[n] - x[n-1];
 252          yi = 2 * y[n] - y[n-1];
 253       } else {
 254          xi = x[i];
 255          yi = y[i];
 256       }
 257
 258       t = fmod(j * interval, 1.);
 259
 260       while (t < 1. && j < m) {
 261          bl1 = (1. - t);
 262          bl2 = t * t;   /* t^2 */
 263          bl4 = t * bl2; /* t^3 */
 264          bl3 = bl4 - bl2;
 265
 266          bl1 = bl1 * bl1 * bl1;
 267          bl2 = 3. * (bl3 - bl2) + 4.;
 268          bl3 = 3. * (  t - bl3) + 1.;
 269
 270          sxy[j++] = (bl1 * xi_3 + bl2 * x[i-2] + bl3 * x[i-1] + bl4 * xi) / 6.; /* x-coordinate */
 271          sxy[j++] = (bl1 * yi_3 + bl2 * y[i-2] + bl3 * y[i-1] + bl4 * yi) / 6.; /* y-coordinate */
 272
 273          t += interval;
 274       }
 275    }
 276    sxy[m]   = x[n];
 277    sxy[m+1] = y[n];
 278    g2_free(x);
 279 }
 280
 281 void g2_b_spline(int id, int n, double *points, int o)
 282
 283 /*
 284  *      FORMAL ARGUMENTS:
 285  *
 286  *      id                      device id
 287  *      n                       number of data points
 288  *      points                  data points (x[i],y[i])
 289  *      o                       number of interpolated points per data point
 290  */
 291
 292 {
 293    int m;
 294    double *sxy;
 295
 296    m = (n-1)*o+1;
 297    sxy = (double*)g2_malloc(m*2*sizeof(double));
 298
 299    g2_c_b_spline(n, points, m, sxy);
 300    g2_poly_line(id, m, sxy);
 301
 302    g2_free(sxy);
 303 }
 304
 305 void g2_filled_b_spline(int id, int n, double *points, int o)
 306
 307 /*
 308  *      FORMAL ARGUMENTS:
 309  *
 310  *      id                      device id
 311  *      n                       number of data points
 312  *      points                  data points (x[i],y[i])
 313  *      o                       number of interpolated points per data point
 314  */
 315
 316 {
 317    int m;
 318    double *sxy;
 319
 320    m = (n-1)*o+1;
 321    sxy = (double*)g2_malloc((m+1)*2*sizeof(double));
 322
 323    g2_c_b_spline(n, points, m, sxy);
 324    sxy[m+m] = points[n+n-2];
 325    sxy[m+m+1] = points[1];
 326    g2_filled_polygon(id, m+1, sxy);
 327
 328    g2_free(sxy);
 329 }
 330
 331 /*
 332  *      FUNCTION g2_c_raspln
 333  *
 334  *      FUNCTIONAL DESCRIPTION:
 335  *
 336  *      This function draws a piecewise cubic polynomial through
 337  *      the specified data points. The (n-1) cubic polynomials are
 338  *      basically parametric cubic Hermite polynomials through the
 339  *      n specified data points with tangent values at the data
 340  *      points determined by a weighted average of the slopes of
 341  *      the secant lines. A tension parameter "tn" is provided to
 342  *      adjust the length of the tangent vector at the data points.
 343  *      This allows the "roundness" of the curve to be adjusted.
 344  *      For further information and references on this technique see:
 345  *
 346  *      D. Kochanek and R. Bartels, Interpolating Splines With Local
 347  *      Tension, Continuity and Bias Control, Computer Graphics,
 348  *      18(1984)3.
 349  *
 350  *      AUTHORS:
 351  *
 352  *      Dennis Mikkelson        distributed in GPLOT    Jan 7, 1988     F77
 353  *      Tijs Michels            t.michels@vimec.nl      Jun 7, 1999     C
 354  *
 355  *      FORMAL ARGUMENTS:
 356  *
 357  *      n       number of data points, n > 2
 358  *      points  double array holding the x and y-coords of the data points
 359  *      tn      double parameter in [0.0, 2.0]. When tn = 0.0,
 360  *              the curve through the data points is very rounded.
 361  *              As tn increases the curve is gradually pulled tighter.
 362  *              When tn = 2.0, the curve is essentially a polyline
 363  *              through the given data points.
 364  *      sxy     double array holding the coords of the spline curve
 365  *
 366  *      IMPLICIT INPUTS:        NONE
 367  *      IMPLICIT OUTPUTS:       NONE
 368  *      SIDE EFFECTS:           NONE
 369  */
 370
 371 #define nb 40
 372 /*
 373  * Number of straight connecting lines of which each polynomial consists.
 374  * So between one data point and the next, (nb-1) points are placed.
 375  */
 376
 377 void g2_c_raspln(int n, const double *points, double tn, double *sxy)
 378 {
 379    int i, j;
 380    double *x, *y;
 381    double bias, tnFactor, tangentL1, tangentL2;
 382    double D1x, D1y, D2x, D2y, t1x, t1y, t2x, t2y;
 383    double h1[nb+1];     /* Values of the Hermite basis functions */
 384    double h2[nb+1];     /* at nb+1 evenly spaced points in [0,1] */
 385    double h3[nb+1];
 386    double h4[nb+1];
 387
 388    x = (double *) g2_malloc(n*2*sizeof(double));
 389    y = x + n;
 390    g2_split(n, points, x, y);
 391
 392 /*
 393  * First, store the values of the Hermite basis functions in a table h[ ]
 394  * so no time is wasted recalculating them
 395  */
 396    for (i = 0; i < nb+1; i++) {
 397       double t, tt, ttt;
 398       t = (double) i / nb;
 399       tt  = t * t;
 400       ttt = t * tt;
 401       h1[i] =  2. * ttt - 3. * tt + 1.;
 402       h2[i] = -2. * ttt + 3. * tt;
 403       h3[i] =       ttt - 2. * tt + t;
 404       h4[i] =       ttt -      tt;
 405    }
 406
 407 /*
 408  * Set local tnFactor based on input parameter tn
 409  */
 410    if (tn <= 0.) {
 411       tnFactor = 2.;
 412       fputs("g2_c_raspln: Using Tension Factor 0.0: very rounded", stderr);
 413    }
 414    else if (tn >= 2.) {
 415       tnFactor = 0.;
 416       fputs("g2_c_raspln: Using Tension Factor 2.0: not rounded at all", stderr);
 417    }
 418    else                 tnFactor = 2. - tn;
 419
 420    D1x = D1y = 0.; /* first point has no preceding point */
 421    for (j = 0; j < n - 2; j++) {
 422       t1x = x[j+1] - x[j];
 423       t1y = y[j+1] - y[j];
 424       t2x = x[j+2] - x[j+1];
 425       t2y = y[j+2] - y[j+1];
 426       tangentL1 = t1x * t1x + t1y * t1y;
 427       tangentL2 = t2x * t2x + t2y * t2y;
 428       if (tangentL1 + tangentL2 == 0) bias = .5;
 429       else bias = tangentL2 / (tangentL1 + tangentL2);
 430       D2x = tnFactor * (bias  * t1x + (1 - bias) * t2x);
 431       D2y = tnFactor * (bias  * t1y + (1 - bias) * t2y);
 432       for (i = 0; i < nb; i++) {
 433         sxy[2 * nb * j + i + i] =
 434            h1[i] * x[j] + h2[i] * x[j+1] + h3[i] * D1x + h4[i] * D2x;
 435         sxy[2 * nb * j + i + i + 1] =
 436            h1[i] * y[j] + h2[i] * y[j+1] + h3[i] * D1y + h4[i] * D2y;
 437       }
 438       D1x = D2x; /* store as preceding point in */
 439       D1y = D2y; /* the next pass */
 440    }
 441
 442 /*
 443  * Do the last subinterval as a special case since no point follows the
 444  * last point
 445  */
 446    for (i = 0; i < nb+1; i++) {
 447       sxy[2 * nb * (n-2) + i + i] =
 448         h1[i] * x[n-2] + h2[i] * x[n-1] + h3[i] * D1x;
 449       sxy[2 * nb * (n-2) + i + i + 1] =
 450         h1[i] * y[n-2] + h2[i] * y[n-1] + h3[i] * D1y;
 451    }
 452    g2_free(x);
 453 }
 454
 455 void g2_raspln(int id, int n, double *points, double tn)
 456
 457 /*
 458  *      FORMAL ARGUMENTS:
 459  *
 460  *      id                      device id
 461  *      n                       number of data points
 462  *      points                  data points (x[i],y[i])
 463  *      tn                      tension factor [0.0, 2.0]
 464  *                              0.0  very rounded
 465  *                              2.0  not rounded at all
 466  */
 467
 468 {
 469    int m;
 470    double *sxy;         /*      coords of the entire spline curve */
 471    m = (n-1)*nb+1;
 472    sxy = (double *) g2_malloc(m*2*sizeof(double));
 473
 474    g2_c_raspln(n, points, tn, sxy);
 475    g2_poly_line(id, m, sxy);
 476
 477    g2_free(sxy);
 478 }
 479
 480 void g2_filled_raspln(int id, int n, double *points, double tn)
 481
 482 /*
 483  *      FORMAL ARGUMENTS:
 484  *
 485  *      id                      device id
 486  *      n                       number of data points
 487  *      points                  data points (x[i],y[i])
 488  *      tn                      tension factor [0.0, 2.0]
 489  *                              0.0  very rounded
 490  *                              2.0  not rounded at all
 491  */
 492
 493 {
 494    int m;
 495    double *sxy;         /*      coords of the entire spline curve */
 496    m = (n-1)*nb+2;
 497    sxy = (double *) g2_malloc(m*2*sizeof(double));
 498
 499    g2_c_raspln(n, points, tn, sxy);
 500    sxy[(n+n-2) * nb + 2] = points[n+n-2];
 501    sxy[(n+n-2) * nb + 3] = points[1];
 502    g2_filled_polygon(id, m, sxy);
 503
 504    g2_free(sxy);
 505 }
 506
 507 /* ---- And now for a rather different approach ---- */
 508
 509 /*
 510  *      FUNCTION g2_c_newton
 511  *
 512  *      FUNCTIONAL DESCRIPTION:
 513  *
 514  *      Use Newton's Divided Differences to calculate an interpolation
 515  *      polynomial through the specified data points.
 516  *      This function is called by
 517  *              g2_c_para_3 and
 518  *              g2_c_para_5.
 519  *
 520  *      Dennis Mikkelson        distributed in GPLOT    Jan  5, 1988    F77
 521  *      Tijs Michels            t.michels@vimec.nl      Jun 16, 1999    C
 522  *
 523  *      FORMAL ARGUMENTS:
 524  *
 525  *      n       number of entries in c1 and c2, 4 <= n <= MaxPts
 526  *              for para_3      (degree 3)      n = 4
 527  *              for para_5      (degree 5)      n = 6
 528  *              for para_i      (degree i)      n = (i + 1)
 529  *      c1      double array holding at most MaxPts values giving the
 530  *              first  coords of the points to be interpolated
 531  *      c2      double array holding at most MaxPts values giving the
 532  *              second coords of the points to be interpolated
 533  *      o       number of points at which the interpolation
 534  *              polynomial is to be evaluated
 535  *      xv      double array holding o points at which to
 536  *              evaluate the interpolation polynomial
 537  *      yv      double array holding upon return the values of the
 538  *              interpolation polynomial at the corresponding points in xv
 539  *
 540  *              yv is the OUTPUT
 541  *
 542  *      IMPLICIT INPUTS:        NONE
 543  *      IMPLICIT OUTPUTS:       NONE
 544  *      SIDE EFFECTS:           NONE
 545  */
 546
 547 #define MaxPts 21
 548 #define xstr(s) __str(s)
 549 #define __str(s) #s
 550
 551 /*
 552  * Maximum number of data points allowed
 553  * 21 would correspond to a polynomial of degree 20
 554  */
 555
 556 void g2_c_newton(int n, const double *c1, const double *c2,
 557                  int o, const double *xv, double *yv)
 558 {
 559    int i, j;
 560    double p, s;
 561    double ddt[MaxPts][MaxPts];          /* Divided Difference Table */
 562
 563    if (n < 4) {
 564       fputs("g2_c_newton: Error! Less than 4 points passed "
 565             "to function g2_c_newton\n", stderr);
 566       return;
 567    }
 568
 569    if (n > MaxPts) {
 570       fputs("g2_c_newton: Error! More than " xstr(MaxPts) " points passed "
 571             "to function g2_c_newton\n", stderr);
 572       return;
 573    }
 574
 575 /* First, build the divided difference table */
 576
 577    for (i = 0; i < n; i++)      ddt[i][0] = c2[i];
 578    for (j = 1; j < n; j++) {
 579       for (i = 0; i < n - j; i++)
 580         ddt[i][j] = (ddt[i+1][j-1] - ddt[i][j-1]) / (c1[i+j] - c1[i]);
 581    }
 582
 583 /* Next, evaluate the polynomial at the specified points */
 584
 585    for (i = 0; i < o; i++) {
 586       for (p = 1., s = ddt[0][0], j = 1; j < n; j++) {
 587          p *= xv[i] - c1[j-1];
 588          s += p * ddt[0][j];
 589       }
 590       yv[i] = s;
 591    }
 592 }
 593
 594 /*
 595  *      FUNCTION: g2_c_para_3
 596  *
 597  *      FUNCTIONAL DESCRIPTION:
 598  *
 599  *      This function draws a piecewise parametric interpolation
 600  *      polynomial of degree 3 through the specified data points.
 601  *      The effect is similar to that obtained using DISSPLA to
 602  *      draw a curve after a call to the DISSPLA routine PARA3.
 603  *      The curve is parameterized using an approximation to the
 604  *      curve's arc length. The basic interpolation is done
 605  *      using function g2_c_newton.
 606  *
 607  *      Dennis Mikkelson        distributed in GPLOT    Jan  7, 1988    F77
 608  *      Tijs Michels            t.michels@vimec.nl      Jun 17, 1999    C
 609  *
 610  *      FORMAL ARGUMENTS:
 611  *
 612  *      n       number of data points through which to draw the curve
 613  *      points  double array containing the x and y-coords of the data points
 614  *
 615  *      IMPLICIT INPUTS:        NONE
 616  *      IMPLICIT OUTPUTS:       NONE
 617  *      SIDE EFFECTS:           NONE
 618  */
 619
 620 /*
 621  * #undef  nb
 622  * #define nb 40
 623  * Number of straight connecting lines of which each polynomial consists.
 624  * So between one data point and the next, (nb-1) points are placed.
 625  */
 626
 627 void g2_c_para_3(int n, const double *points, double *sxy)
 628 {
 629 #define dgr     (3+1)
 630 #define nb2     (nb*2)
 631    int i, j;
 632    double x1t, y1t;
 633    double o, step;
 634    double X[nb2];               /* x-coords of the current curve piece */
 635    double Y[nb2];               /* y-coords of the current curve piece */
 636    double t[dgr];               /* data point parameter values */
 637    double Xpts[dgr];            /* x-coords data point subsection */
 638    double Ypts[dgr];            /* y-coords data point subsection */
 639    double s[nb2];               /* parameter values at which to interpolate */
 640
 641    /* Do first TWO subintervals first */
 642
 643    g2_split(dgr, points, Xpts, Ypts);
 644
 645    t[0] = 0.;
 646    for (i = 1; i < dgr; i++) {
 647       x1t = Xpts[i] - Xpts[i-1];
 648       y1t = Ypts[i] - Ypts[i-1];
 649       t[i] = t[i-1] + sqrt(x1t * x1t + y1t * y1t);
 650    }
 651
 652    step = t[2] / nb2;
 653    for (i = 0; i < nb2; i++)    s[i] = i * step;
 654
 655    g2_c_newton(dgr, t, Xpts, nb2, s, X);
 656    g2_c_newton(dgr, t, Ypts, nb2, s, Y);
 657    for (i = 0; i < nb2; i++) {
 658       sxy[i+i]   = X[i];
 659       sxy[i+i+1] = Y[i];
 660    }
 661
 662    /* Next, do later central subintervals */
 663
 664    for (j = 1; j < n - dgr + 1; j++) {
 665       g2_split(dgr, points + j + j, Xpts, Ypts);
 666
 667       for (i = 1; i < dgr; i++) {
 668          x1t = Xpts[i] - Xpts[i-1];
 669          y1t = Ypts[i] - Ypts[i-1];
 670          t[i] = t[i-1] + sqrt(x1t * x1t + y1t * y1t);
 671       }
 672
 673       o = t[1]; /* look up once */
 674       step = (t[2] - o) / nb;
 675       for (i = 0; i < nb; i++)  s[i] = i * step + o;
 676
 677       g2_c_newton(dgr, t, Xpts, nb, s, X);
 678       g2_c_newton(dgr, t, Ypts, nb, s, Y);
 679
 680       for (i = 0; i < nb; i++) {
 681          sxy[(j + 1) * nb2 + i + i]     = X[i];
 682          sxy[(j + 1) * nb2 + i + i + 1] = Y[i];
 683       }
 684    }
 685
 686    /* Now do last subinterval */
 687
 688    o = t[2];
 689    step = (t[3] - o) / nb;
 690    for (i = 0; i < nb; i++)     s[i] = i * step + o;
 691
 692    g2_c_newton(dgr, t, Xpts, nb, s, X);
 693    g2_c_newton(dgr, t, Ypts, nb, s, Y);
 694
 695    for (i = 0; i < nb; i++) {
 696       sxy[(n - dgr + 2) * nb2 + i + i]     = X[i];
 697       sxy[(n - dgr + 2) * nb2 + i + i + 1] = Y[i];
 698    }
 699    sxy[(n - 1) * nb2]     = points[n+n-2];
 700    sxy[(n - 1) * nb2 + 1] = points[n+n-1];
 701 }
 702
 703 /*
 704  *      FORMAL ARGUMENTS:
 705  *
 706  *      id                      device id
 707  *      n                       number of data points
 708  *      points                  data points (x[i],y[i])
 709  */
 710
 711 void g2_para_3(int id, int n, double *points)
 712 {
 713    int m;
 714    double *sxy;         /*      coords of the entire spline curve */
 715    m = (n-1)*nb+1;
 716    sxy = (double *) g2_malloc(m*2*sizeof(double));
 717
 718    g2_c_para_3(n, points, sxy);
 719    g2_poly_line(id, m, sxy);
 720
 721    g2_free(sxy);
 722 }
 723
 724 /*
 725  *      FORMAL ARGUMENTS:
 726  *
 727  *      id                      device id
 728  *      n                       number of data points
 729  *      points                  data points (x[i],y[i])
 730  */
 731
 732 void g2_filled_para_3(int id, int n, double *points)
 733 {
 734    int m;
 735    double *sxy;         /*      coords of the entire spline curve */
 736    m = (n-1)*nb+2;
 737    sxy = (double *) g2_malloc(m*2*sizeof(double));
 738
 739    g2_c_para_3(n, points, sxy);
 740    sxy[m+m-2] = points[n+n-2];
 741    sxy[m+m-1] = points[1];
 742    g2_filled_polygon(id, m, sxy);
 743
 744    g2_free(sxy);
 745 }
 746
 747 /*
 748  *      FUNCTION: g2_c_para_5
 749  *
 750  *      As g2_c_para_3, but now plot a polynomial of degree 5
 751  */
 752
 753 /*
 754  * #undef  nb
 755  * #define nb 40
 756  * Number of straight connecting lines of which each polynomial consists.
 757  * So between one data point and the next, (nb-1) points are placed.
 758  */
 759
 760 void g2_c_para_5(int n, const double *points, double *sxy)
 761 {
 762 #undef  dgr
 763 #define dgr     (5+1)
 764 #define nb3     (nb*3)
 765    int i, j;
 766    double x1t, y1t;
 767    double o, step;
 768    double X[nb3];               /* x-coords of the current curve piece */
 769    double Y[nb3];               /* y-coords of the current curve piece */
 770    double t[dgr];               /* data point parameter values */
 771    double Xpts[dgr];            /* x-coords data point subsection */
 772    double Ypts[dgr];            /* y-coords data point subsection */
 773    double s[nb3];               /* parameter values at which to interpolate */
 774
 775    /* Do first THREE subintervals first */
 776
 777    g2_split(dgr, points, Xpts, Ypts);
 778
 779    t[0] = 0.;
 780    for (i = 1; i < dgr; i++) {
 781       x1t = Xpts[i] - Xpts[i-1];
 782       y1t = Ypts[i] - Ypts[i-1];
 783       t[i] = t[i-1] + sqrt(x1t * x1t + y1t * y1t);
 784    }
 785
 786    step = t[3] / nb3;
 787    for (i = 0; i < nb3; i++)    s[i] = i * step;
 788
 789    g2_c_newton(dgr, t, Xpts, nb3, s, X);
 790    g2_c_newton(dgr, t, Ypts, nb3, s, Y);
 791    for (i = 0; i < nb3; i++) {
 792       sxy[i+i]   = X[i];
 793       sxy[i+i+1] = Y[i];
 794    }
 795
 796    /* Next, do later central subintervals */
 797
 798    for (j = 1; j < n - dgr + 1; j++) {
 799       g2_split(dgr, points + j + j, Xpts, Ypts);
 800
 801       for (i = 1; i < dgr; i++) {
 802          x1t = Xpts[i] - Xpts[i-1];
 803          y1t = Ypts[i] - Ypts[i-1];
 804          t[i] = t[i-1] + sqrt(x1t * x1t + y1t * y1t);
 805       }
 806
 807       o = t[2]; /* look up once */
 808       step = (t[3] - o) / nb;
 809       for (i = 0; i < nb; i++)  s[i] = i * step + o;
 810
 811       g2_c_newton(dgr, t, Xpts, nb, s, X);
 812       g2_c_newton(dgr, t, Ypts, nb, s, Y);
 813
 814       for (i = 0; i < nb; i++) {
 815          sxy[(j + 2) * nb2 + i + i]     = X[i];
 816          sxy[(j + 2) * nb2 + i + i + 1] = Y[i];
 817       }
 818    }
 819
 820    /* Now do last TWO subinterval */
 821
 822    o = t[3];
 823    step = (t[5] - o) / nb2;
 824    for (i = 0; i < nb2; i++)    s[i] = i * step + o;
 825
 826    g2_c_newton(dgr, t, Xpts, nb2, s, X);
 827    g2_c_newton(dgr, t, Ypts, nb2, s, Y);
 828
 829    for (i = 0; i < nb2; i++) {
 830       sxy[(n - dgr + 3) * nb2 + i + i]     = X[i];
 831       sxy[(n - dgr + 3) * nb2 + i + i + 1] = Y[i];
 832    }
 833    sxy[(n - 1) * nb2]     = points[n+n-2];
 834    sxy[(n - 1) * nb2 + 1] = points[n+n-1];
 835 }
 836
 837 /*
 838  *      FORMAL ARGUMENTS:
 839  *
 840  *      id                      device id
 841  *      n                       number of data points
 842  *      points                  data points (x[i],y[i])
 843  */
 844
 845 void g2_para_5(int id, int n, double *points)
 846 {
 847    int m;
 848    double *sxy;         /*      coords of the entire spline curve */
 849    m = (n-1)*nb+1;
 850    sxy = (double *) g2_malloc(m*2*sizeof(double));
 851
 852    g2_c_para_5(n, points, sxy);
 853    g2_poly_line(id, m, sxy);
 854
 855    g2_free(sxy);
 856 }
 857
 858 /*
 859  *      FORMAL ARGUMENTS:
 860  *
 861  *      id                      device id
 862  *      n                       number of data points
 863  *      points                  data points (x[i],y[i])
 864  */
 865
 866 void g2_filled_para_5(int id, int n, double *points)
 867 {
 868    int m;
 869    double *sxy;         /*      coords of the entire spline curve */
 870    m = (n-1)*nb+2;
 871    sxy = (double *) g2_malloc(m*2*sizeof(double));
 872
 873    g2_c_para_5(n, points, sxy);
 874    sxy[m+m-2] = points[n+n-2];
 875    sxy[m+m-1] = points[1];
 876    g2_filled_polygon(id, m, sxy);
 877
 878    g2_free(sxy);
 879 }
 880