git.maemo.org Git - gnuplot/blob - src/interpol.c

   1 #ifndef lint
   2 static char *RCSid() { return RCSid("$Id: interpol.c,v 1.31 2004/07/25 12:25:01 broeker Exp $"); }
   3 #endif
   4
   5 /* GNUPLOT - interpol.c */
   6
   7 /*[
   8  * Copyright 1986 - 1993, 1998, 2004   Thomas Williams, Colin Kelley
   9  *
  10  * Permission to use, copy, and distribute this software and its
  11  * documentation for any purpose with or without fee is hereby granted,
  12  * provided that the above copyright notice appear in all copies and
  13  * that both that copyright notice and this permission notice appear
  14  * in supporting documentation.
  15  *
  16  * Permission to modify the software is granted, but not the right to
  17  * distribute the complete modified source code.  Modifications are to
  18  * be distributed as patches to the released version.  Permission to
  19  * distribute binaries produced by compiling modified sources is granted,
  20  * provided you
  21  *   1. distribute the corresponding source modifications from the
  22  *    released version in the form of a patch file along with the binaries,
  23  *   2. add special version identification to distinguish your version
  24  *    in addition to the base release version number,
  25  *   3. provide your name and address as the primary contact for the
  26  *    support of your modified version, and
  27  *   4. retain our contact information in regard to use of the base
  28  *    software.
  29  * Permission to distribute the released version of the source code along
  30  * with corresponding source modifications in the form of a patch file is
  31  * granted with same provisions 2 through 4 for binary distributions.
  32  *
  33  * This software is provided "as is" without express or implied warranty
  34  * to the extent permitted by applicable law.
  35 ]*/
  36
  37
  38 /*
  39  * C-Source file identification Header
  40  *
  41  * This file belongs to a project which is:
  42  *
  43  * done 1993 by MGR-Software, Asgard  (Lars Hanke)
  44  * written by Lars Hanke
  45  *
  46  * Contact me via:
  47  *
  48  *  InterNet: mgr@asgard.bo.open.de
  49  *      FIDO: Lars Hanke @ 2:243/4802.22   (as long as they keep addresses)
  50  *
  51  **************************************************************************
  52  *
  53  *   Project: gnuplot
  54  *    Module:
  55  *      File: interpol.c
  56  *
  57  *   Revisor: Lars Hanke
  58  *   Revised: 26/09/93
  59  *  Revision: 1.0
  60  *
  61  **************************************************************************
  62  *
  63  * LEGAL
  64  *  This module is part of gnuplot and distributed under whatever terms
  65  *  gnuplot is or will be published, unless exclusive rights are claimed.
  66  *
  67  * DESCRIPTION
  68  *  Supplies 2-D data interpolation and approximation routines
  69  *
  70  * IMPORTS
  71  *  plot.h
  72  *    - cp_extend()
  73  *    - structs: curve_points, coordval, coordinate
  74  *
  75  *  setshow.h
  76  *    - samples, axis array[] variables
  77  *    - plottypes
  78  *
  79  *  proto.h
  80  *    - solve_tri_diag()
  81  *    - typedef tri_diag
  82  *
  83  * EXPORTS
  84  *  gen_interp()
  85  *  sort_points()
  86  *  cp_implode()
  87  *
  88  * BUGS and TODO
  89  *  I would really have liked to use Gershon Elbers contouring code for
  90  *  all the stuff done here, but I failed. So I used my own code.
  91  *  If somebody is able to consolidate Gershon's code for this purpose
  92  *  a lot of gnuplot users would be very happy - due to memory problems.
  93  *
  94  **************************************************************************
  95  *
  96  * HISTORY
  97  * Changes:
  98  *  Nov 24, 1995  Markus Schuh (M.Schuh@meteo.uni-koeln.de):
  99  *      changed the algorithm for csplines
 100  *      added algorithm for approximation csplines
 101  *      copied point storage and range fix from plot2d.c
 102  *
 103  *  Dec 12, 1995 David Denholm
 104  *      oops - at the time this is called, stored co-ords are
 105  *      internal (ie maybe log of data) but min/max are in
 106  *      user co-ordinates.
 107  *      Work with min and max of internal co-ords, and
 108  *      check at the end whether external min and max need to
 109  *      be increased. (since samples_1 is typically 100 ; we
 110  *      dont want to take more logs than necessary)
 111  *      Also, need to take into account which axes are active
 112  *
 113  *  Jun 30, 1996 Jens Emmerich
 114  *      implemented handling of UNDEFINED points
 115  */
 116
 117 #include "interpol.h"
 118
 119 #include "alloc.h"
 120 #include "axis.h"
 121 #include "contour.h"
 122 #include "graphics.h"
 123 #include "misc.h"
 124 #include "plot2d.h"
 125 /*  #include "setshow.h" */
 126 #include "util.h"
 127
 128
 129 /* in order to support multiple axes, and to simplify ranging in
 130  * parametric plots, we use arrays to store some things. For 2d plots,
 131  * elements are z=0,y1=1,x1=2,z2=4,y2=5,x2=6 these are given symbolic
 132  * names in plot.h
 133  */
 134
 135
 136 /*
 137  * IMHO, code is getting too cluttered with repeated chunks of
 138  * code. Some macros to simplify, I hope.
 139  *
 140  * do { } while(0) is comp.lang.c recommendation for complex macros
 141  * also means that break can be specified as an action, and it will
 142  *
 143  */
 144
 145
 146 /* store VALUE or log(VALUE) in STORE, set TYPE as appropriate Do
 147  * OUT_ACTION or UNDEF_ACTION as appropriate. Adjust range provided
 148  * type is INRANGE (ie dont adjust y if x is outrange). VALUE must not
 149  * be same as STORE */
 150 /* FIXME 20010610: UNDEF_ACTION is completely unused ??? Furthermore,
 151  * this is so similar to STORE_WITH_LOG_AND_UPDATE_RANGE() from axis.h
 152  * that the two should probably be merged.  */
 153 #define STORE_AND_FIXUP_RANGE(store, value, type, min, max, auto,       \
 154                               out_action, undef_action)                 \
 155 do {                                                                    \
 156     store=value;                                                        \
 157     if (type != INRANGE)                                                \
 158         break;  /* don't set y range if x is outrange, for example */   \
 159     if ((value) < (min)) {                                              \
 160        if ((auto) & AUTOSCALE_MIN)                                      \
 161            (min) = (value);                                             \
 162        else {                                                           \
 163            (type) = OUTRANGE;                                           \
 164            out_action;                                                  \
 165            break;                                                       \
 166        }                                                                \
 167     }                                                                   \
 168     if ((value) > (max)) {                                              \
 169        if ((auto) & AUTOSCALE_MAX)                                      \
 170            (max) = (value);                                             \
 171        else {                                                           \
 172            (type) = OUTRANGE;                                           \
 173            out_action;                                                  \
 174        }                                                                \
 175     }                                                                   \
 176 } while(0)
 177
 178 #define UPDATE_RANGE(TEST,OLD,NEW,AXIS)         \
 179 do {                                            \
 180     if (TEST)                                   \
 181         (OLD) = AXIS_DE_LOG_VALUE(AXIS,NEW);    \
 182 } while(0)
 183
 184 #define spline_coeff_size 4
 185 typedef double spline_coeff[spline_coeff_size];
 186 typedef double five_diag[5];
 187
 188 static int next_curve __PROTO((struct curve_points * plot, int *curve_start));
 189 static int num_curves __PROTO((struct curve_points * plot));
 190 static double *cp_binomial __PROTO((int points));
 191 static void eval_bezier __PROTO((struct curve_points * cp, int first_point,
 192                                  int num_points, double sr, coordval * px,
 193                                  coordval *py, double *c));
 194 static void do_bezier __PROTO((struct curve_points * cp, double *bc, int first_point, int num_points, struct coordinate * dest));
 195 static int solve_five_diag __PROTO((five_diag m[], double r[], double x[], int n));
 196 static spline_coeff *cp_approx_spline __PROTO((struct curve_points * plot, int first_point, int num_points));
 197 static spline_coeff *cp_tridiag __PROTO((struct curve_points * plot, int first_point, int num_points));
 198 static void do_cubic __PROTO((struct curve_points * plot, spline_coeff * sc, int first_point, int num_points, struct coordinate * dest));
 199 static void do_freq __PROTO((struct curve_points *plot, int first_point, int num_points));
 200 int compare_points __PROTO((SORTFUNC_ARGS p1, SORTFUNC_ARGS p2));
 201
 202
 203 /*
 204  * position curve_start to index the next non-UNDEFINDED point,
 205  * start search at initial curve_start,
 206  * return number of non-UNDEFINDED points from there on,
 207  * if no more valid points are found, curve_start is set
 208  * to plot->p_count and 0 is returned
 209  */
 210
 211 static int
 212 next_curve(struct curve_points *plot, int *curve_start)
 213 {
 214     int curve_length;
 215
 216     /* Skip undefined points */
 217     while (*curve_start < plot->p_count
 218            && plot->points[*curve_start].type == UNDEFINED) {
 219         (*curve_start)++;
 220     };
 221     curve_length = 0;
 222     /* curve_length is first used as an offset, then the correct # points */
 223     while ((*curve_start) + curve_length < plot->p_count
 224            && plot->points[(*curve_start) + curve_length].type != UNDEFINED) {
 225         curve_length++;
 226     };
 227     return (curve_length);
 228 }
 229
 230
 231 /*
 232  * determine the number of curves in plot->points, separated by
 233  * UNDEFINED points
 234  */
 235
 236 static int
 237 num_curves(struct curve_points *plot)
 238 {
 239     int curves;
 240     int first_point;
 241     int num_points;
 242
 243     first_point = 0;
 244     curves = 0;
 245     while ((num_points = next_curve(plot, &first_point)) > 0) {
 246         curves++;
 247         first_point += num_points;
 248     }
 249     return (curves);
 250 }
 251
 252
 253
 254 /*
 255  * build up a cntr_struct list from curve_points
 256  * this funtion is only used for the alternate entry point to
 257  * Gershon's code and thus commented out
 258  ***deleted***
 259  */
 260
 261
 262 /* HBB 990205: rewrote the 'bezier' interpolation routine,
 263  * to prevent numerical overflow and other undesirable things happening
 264  * for large data files (num_data about 1000 or so), where binomial
 265  * coefficients would explode, and powers of 'sr' (0 < sr < 1) become
 266  * extremely small. Method used: compute logarithms of these
 267  * extremely large and small numbers, and only go back to the
 268  * real numbers once they've cancelled out each other, leaving
 269  * a reasonable-sized one. */
 270
 271 /*
 272  * cp_binomial() computes the binomial coefficients needed for BEZIER stuff
 273  *   and stores them into an array which is hooked to sdat.
 274  * (MGR 1992)
 275  */
 276 static double *
 277 cp_binomial(int points)
 278 {
 279     double *coeff;
 280     int n, k;
 281     int e;
 282
 283     e = points;                 /* well we're going from k=0 to k=p_count-1 */
 284     coeff = gp_alloc(e * sizeof(double), "bezier coefficients");
 285
 286     n = points - 1;
 287     e = n / 2;
 288     /* HBB 990205: calculate these in 'logarithmic space',
 289      * as they become _very_ large, with growing n (4^n) */
 290     coeff[0] = 0.0;
 291
 292     for (k = 0; k < e; k++) {
 293         coeff[k + 1] = coeff[k] + log(((double) (n - k)) / ((double) (k + 1)));
 294     }
 295
 296     for (k = n; k >= e; k--)
 297         coeff[k] = coeff[n - k];
 298
 299     return (coeff);
 300 }
 301
 302
 303 /* This is a subfunction of do_bezier() for BEZIER style computations.
 304  * It is passed the stepration (STEP/MAXSTEPS) and the addresses of
 305  * the double values holding the next x and y coordinates.
 306  * (MGR 1992)
 307  */
 308
 309 static void
 310 eval_bezier(
 311     struct curve_points *cp,
 312     int first_point,            /* where to start in plot->points (to find x-range) */
 313     int num_points,             /* to determine end in plot->points */
 314     double sr,                  /* position inside curve, range [0:1] */
 315     coordval *px,               /* OUTPUT: x and y */
 316     coordval *py,
 317     double *c)                  /* Bezier coefficient array */
 318 {
 319     unsigned int n = num_points - 1;
 320     struct coordinate GPHUGE *this_points;
 321
 322     this_points = (cp->points) + first_point;
 323
 324     if (sr == 0.0) {
 325         *px = this_points[0].x;
 326         *py = this_points[0].y;
 327     } else if (sr == 1.0) {
 328         *px = this_points[n].x;
 329         *py = this_points[n].y;
 330     } else {
 331         /* HBB 990205: do calculation in 'logarithmic space',
 332          * to avoid over/underflow errors, which would exactly cancel
 333          * out each other, anyway, in an exact calculation
 334          */
 335         unsigned int i;
 336         double lx = 0.0, ly = 0.0;
 337         double log_dsr_to_the_n = n * log(1 - sr);
 338         double log_sr_over_dsr = log(sr) - log(1 - sr);
 339
 340         for (i = 0; i <= n; i++) {
 341             double u = exp(c[i] + log_dsr_to_the_n + i * log_sr_over_dsr);
 342
 343             lx += this_points[i].x * u;
 344             ly += this_points[i].y * u;
 345         }
 346
 347         *px = lx;
 348         *py = ly;
 349     }
 350 }
 351
 352 /*
 353  * generate a new set of coordinates representing the bezier curve and
 354  * set it to the plot
 355  */
 356
 357 static void
 358 do_bezier(
 359     struct curve_points *cp,
 360     double *bc,                 /* Bezier coefficient array */
 361     int first_point,            /* where to start in plot->points */
 362     int num_points,             /* to determine end in plot->points */
 363     struct coordinate *dest)    /* where to put the interpolated data */
 364 {
 365     int i;
 366     coordval x, y;
 367
 368     /* min and max in internal (eg logged) co-ordinates. We update
 369      * these, then update the external extrema in user co-ordinates
 370      * at the end.
 371      */
 372
 373     double ixmin, ixmax, iymin, iymax;
 374     double sxmin, sxmax, symin, symax;  /* starting values of above */
 375
 376     x_axis = cp->x_axis;
 377     y_axis = cp->y_axis;
 378
 379     ixmin = sxmin = AXIS_LOG_VALUE(x_axis, X_AXIS.min);
 380     ixmax = sxmax = AXIS_LOG_VALUE(x_axis, X_AXIS.max);
 381     iymin = symin = AXIS_LOG_VALUE(y_axis, Y_AXIS.min);
 382     iymax = symax = AXIS_LOG_VALUE(y_axis, Y_AXIS.max);
 383
 384     for (i = 0; i < samples_1; i++) {
 385         eval_bezier(cp, first_point, num_points,
 386                     (double) i / (double) (samples_1 - 1),
 387                     &x, &y, bc);
 388
 389         /* now we have to store the points and adjust the ranges */
 390
 391         dest[i].type = INRANGE;
 392         STORE_AND_FIXUP_RANGE(dest[i].x, x, dest[i].type, ixmin, ixmax, X_AXIS.autoscale, NOOP, continue);
 393         STORE_AND_FIXUP_RANGE(dest[i].y, y, dest[i].type, iymin, iymax, Y_AXIS.autoscale, NOOP, NOOP);
 394
 395         dest[i].xlow = dest[i].xhigh = dest[i].x;
 396         dest[i].ylow = dest[i].yhigh = dest[i].y;
 397
 398         dest[i].z = -1;
 399     }
 400
 401     UPDATE_RANGE(ixmax > sxmax, X_AXIS.max, ixmax, x_axis);
 402     UPDATE_RANGE(ixmin < sxmin, X_AXIS.min, ixmin, x_axis);
 403     UPDATE_RANGE(iymax > symax, Y_AXIS.max, iymax, y_axis);
 404     UPDATE_RANGE(iymin < symin, Y_AXIS.min, iymin, y_axis);
 405 }
 406
 407 /*
 408  * call contouring routines -- main entry
 409  */
 410
 411 /*
 412  * it should be like this, but it doesn't run. If you find out why,
 413  * contact me: mgr@asgard.bo.open.de or Lars Hanke 2:243/4802.22@fidonet
 414  *
 415  * Well, all this had originally been inside contour.c, so maybe links
 416  * to functions and of contour.c are broken.
 417  * ***deleted***
 418  * end of unused entry point to Gershon's code
 419  *
 420  */
 421
 422 /*
 423  * Solve five diagonal linear system equation. The five diagonal matrix is
 424  * defined via matrix M, right side is r, and solution X i.e. M * X = R.
 425  * Size of system given in n. Return TRUE if solution exist.
 426  *  G. Engeln-Muellges/ F.Reutter:
 427  *  "Formelsammlung zur Numerischen Mathematik mit Standard-FORTRAN-Programmen"
 428  *  ISBN 3-411-01677-9
 429  *
 430  * /  m02 m03 m04   0   0   0   0    .       .       . \   /  x0  \    / r0  \
 431  * I  m11 m12 m13 m14   0   0   0    .       .       . I   I  x1  I   I  r1  I
 432  * I  m20 m21 m22 m23 m24   0   0    .       .       . I * I  x2  I = I  r2  I
 433  * I    0 m30 m31 m32 m33 m34   0    .       .       . I   I  x3  I   I  r3  I
 434  *      .   .   .   .   .   .   .    .       .       .        .        .
 435  * \                           m(n-3)0 m(n-2)1 m(n-1)2 /   \x(n-1)/   \r(n-1)/
 436  *
 437  */
 438 static int
 439 solve_five_diag(five_diag m[], double r[], double x[], int n)
 440 {
 441     int i;
 442     five_diag *hv;
 443
 444     hv = gp_alloc((n + 1) * sizeof(five_diag), "five_diag help vars");
 445
 446     hv[0][0] = m[0][2];
 447     if (hv[0][0] == 0) {
 448         free(hv);
 449         return FALSE;
 450     }
 451     hv[0][1] = m[0][3] / hv[0][0];
 452     hv[0][2] = m[0][4] / hv[0][0];
 453
 454     hv[1][3] = m[1][1];
 455     hv[1][0] = m[1][2] - hv[1][3] * hv[0][1];
 456     if (hv[1][0] == 0) {
 457         free(hv);
 458         return FALSE;
 459     }
 460     hv[1][1] = (m[1][3] - hv[1][3] * hv[0][2]) / hv[1][0];
 461     hv[1][2] = m[1][4] / hv[1][0];
 462
 463     for (i = 2; i < n; i++) {
 464         hv[i][3] = m[i][1] - m[i][0] * hv[i - 2][1];
 465         hv[i][0] = m[i][2] - m[i][0] * hv[i - 2][2] - hv[i][3] * hv[i - 1][1];
 466         if (hv[i][0] == 0) {
 467             free(hv);
 468             return FALSE;
 469         }
 470         hv[i][1] = (m[i][3] - hv[i][3] * hv[i - 1][2]) / hv[i][0];
 471         hv[i][2] = m[i][4] / hv[i][0];
 472     }
 473
 474     hv[0][4] = 0;
 475     hv[1][4] = r[0] / hv[0][0];
 476     for (i = 1; i < n; i++) {
 477         hv[i + 1][4] = (r[i] - m[i][0] * hv[i - 1][4] - hv[i][3] * hv[i][4]) / hv[i][0];
 478     }
 479
 480     x[n - 1] = hv[n][4];
 481     x[n - 2] = hv[n - 1][4] - hv[n - 2][1] * x[n - 1];
 482     for (i = n - 3; i >= 0; i--)
 483         x[i] = hv[i + 1][4] - hv[i][1] * x[i + 1] - hv[i][2] * x[i + 2];
 484
 485     free(hv);
 486     return TRUE;
 487 }
 488
 489
 490 /*
 491  * Calculation of approximation cubic splines
 492  * Input:  x[i], y[i], weights z[i]
 493  *
 494  * Returns matrix of spline coefficients
 495  */
 496 static spline_coeff *
 497 cp_approx_spline(
 498     struct curve_points *plot,
 499     int first_point,            /* where to start in plot->points */
 500     int num_points)             /* to determine end in plot->points */
 501 {
 502     spline_coeff *sc;
 503     five_diag *m;
 504     double *r, *x, *h, *xp, *yp;
 505     struct coordinate GPHUGE *this_points;
 506     int i;
 507
 508     x_axis = plot->x_axis;
 509     y_axis = plot->y_axis;
 510
 511     sc = gp_alloc((num_points) * sizeof(spline_coeff),
 512                                    "spline matrix");
 513
 514     if (num_points < 4)
 515         int_error(plot->token, "Can't calculate approximation splines, need at least 4 points");
 516
 517     this_points = (plot->points) + first_point;
 518
 519     for (i = 0; i < num_points; i++)
 520         if (this_points[i].z <= 0)
 521             int_error(plot->token, "Can't calculate approximation splines, all weights have to be > 0");
 522
 523     m = gp_alloc((num_points - 2) * sizeof(five_diag), "spline help matrix");
 524
 525     r = gp_alloc((num_points - 2) * sizeof(double), "spline right side");
 526     x = gp_alloc((num_points - 2) * sizeof(double), "spline solution vector");
 527     h = gp_alloc((num_points - 1) * sizeof(double), "spline help vector");
 528
 529     xp = gp_alloc((num_points) * sizeof(double), "x pos");
 530     yp = gp_alloc((num_points) * sizeof(double), "y pos");
 531
 532     /* KB 981107: With logarithmic axis first convert back to linear scale */
 533
 534     xp[0] = AXIS_DE_LOG_VALUE(x_axis, this_points[0].x);
 535     yp[0] = AXIS_DE_LOG_VALUE(y_axis, this_points[0].y);
 536     for (i = 1; i < num_points; i++) {
 537         xp[i] = AXIS_DE_LOG_VALUE(x_axis, this_points[i].x);
 538         yp[i] = AXIS_DE_LOG_VALUE(y_axis, this_points[i].y);
 539         h[i - 1] = xp[i] - xp[i - 1];
 540     }
 541
 542     /* set up the matrix and the vector */
 543
 544     for (i = 0; i <= num_points - 3; i++) {
 545         r[i] = 3 * ((yp[i + 2] - yp[i + 1]) / h[i + 1]
 546                     - (yp[i + 1] - yp[i]) / h[i]);
 547
 548         if (i < 2)
 549             m[i][0] = 0;
 550         else
 551             m[i][0] = 6 / this_points[i].z / h[i - 1] / h[i];
 552
 553         if (i < 1)
 554             m[i][1] = 0;
 555         else
 556             m[i][1] = h[i] - 6 / this_points[i].z / h[i] * (1 / h[i - 1] + 1 / h[i])
 557                 - 6 / this_points[i + 1].z / h[i] * (1 / h[i] + 1 / h[i + 1]);
 558
 559         m[i][2] = 2 * (h[i] + h[i + 1])
 560             + 6 / this_points[i].z / h[i] / h[i]
 561             + 6 / this_points[i + 1].z * (1 / h[i] + 1 / h[i + 1]) * (1 / h[i] + 1 / h[i + 1])
 562             + 6 / this_points[i + 2].z / h[i + 1] / h[i + 1];
 563
 564         if (i > num_points - 4)
 565             m[i][3] = 0;
 566         else
 567             m[i][3] = h[i + 1] - 6 / this_points[i + 1].z / h[i + 1] * (1 / h[i] + 1 / h[i + 1])
 568                 - 6 / this_points[i + 2].z / h[i + 1] * (1 / h[i + 1] + 1 / h[i + 2]);
 569
 570         if (i > num_points - 5)
 571             m[i][4] = 0;
 572         else
 573             m[i][4] = 6 / this_points[i + 2].z / h[i + 1] / h[i + 2];
 574     }
 575
 576     /* solve the matrix */
 577     if (!solve_five_diag(m, r, x, num_points - 2)) {
 578         free(h);
 579         free(x);
 580         free(r);
 581         free(m);
 582         free(xp);
 583         free(yp);
 584         int_error(plot->token, "Can't calculate approximation splines");
 585     }
 586     sc[0][2] = 0;
 587     for (i = 1; i <= num_points - 2; i++)
 588         sc[i][2] = x[i - 1];
 589     sc[num_points - 1][2] = 0;
 590
 591     sc[0][0] = yp[0] + 2 / this_points[0].z / h[0] * (sc[0][2] - sc[1][2]);
 592     for (i = 1; i <= num_points - 2; i++)
 593         sc[i][0] = yp[i] - 2 / this_points[i].z *
 594             (sc[i - 1][2] / h[i - 1]
 595              - sc[i][2] * (1 / h[i - 1] + 1 / h[i])
 596              + sc[i + 1][2] / h[i]);
 597     sc[num_points - 1][0] = yp[num_points - 1]
 598         - 2 / this_points[num_points - 1].z / h[num_points - 2]
 599         * (sc[num_points - 2][2] - sc[num_points - 1][2]);
 600
 601     for (i = 0; i <= num_points - 2; i++) {
 602         sc[i][1] = (sc[i + 1][0] - sc[i][0]) / h[i]
 603             - h[i] / 3 * (sc[i + 1][2] + 2 * sc[i][2]);
 604         sc[i][3] = (sc[i + 1][2] - sc[i][2]) / 3 / h[i];
 605     }
 606
 607     free(h);
 608     free(x);
 609     free(r);
 610     free(m);
 611     free(xp);
 612     free(yp);
 613
 614     return (sc);
 615 }
 616
 617
 618 /*
 619  * Calculation of cubic splines
 620  *
 621  * This can be treated as a special case of approximation cubic splines, with
 622  * all weights -> infinity.
 623  *
 624  * Returns matrix of spline coefficients
 625  */
 626 static spline_coeff *
 627 cp_tridiag(struct curve_points *plot, int first_point, int num_points)
 628 {
 629     spline_coeff *sc;
 630     tri_diag *m;
 631     double *r, *x, *h, *xp, *yp;
 632     struct coordinate GPHUGE *this_points;
 633     int i;
 634
 635     x_axis = plot->x_axis;
 636     y_axis = plot->y_axis;
 637     if (num_points < 3)
 638         int_error(plot->token, "Can't calculate splines, need at least 3 points");
 639
 640     this_points = (plot->points) + first_point;
 641
 642     sc = gp_alloc((num_points) * sizeof(spline_coeff), "spline matrix");
 643     m = gp_alloc((num_points - 2) * sizeof(tri_diag), "spline help matrix");
 644
 645     r = gp_alloc((num_points - 2) * sizeof(double), "spline right side");
 646     x = gp_alloc((num_points - 2) * sizeof(double), "spline solution vector");
 647     h = gp_alloc((num_points - 1) * sizeof(double), "spline help vector");
 648
 649     xp = gp_alloc((num_points) * sizeof(double), "x pos");
 650     yp = gp_alloc((num_points) * sizeof(double), "y pos");
 651
 652     /* KB 981107: With logarithmic axis first convert back to linear scale */
 653
 654     xp[0] = AXIS_DE_LOG_VALUE(x_axis,this_points[0].x);
 655     yp[0] = AXIS_DE_LOG_VALUE(y_axis,this_points[0].y);
 656     for (i = 1; i < num_points; i++) {
 657         xp[i] = AXIS_DE_LOG_VALUE(x_axis,this_points[i].x);
 658         yp[i] = AXIS_DE_LOG_VALUE(y_axis,this_points[i].y);
 659         h[i - 1] = xp[i] - xp[i - 1];
 660     }
 661
 662     /* set up the matrix and the vector */
 663
 664     for (i = 0; i <= num_points - 3; i++) {
 665         r[i] = 3 * ((yp[i + 2] - yp[i + 1]) / h[i + 1]
 666                     - (yp[i + 1] - yp[i]) / h[i]);
 667
 668         if (i < 1)
 669             m[i][0] = 0;
 670         else
 671             m[i][0] = h[i];
 672
 673         m[i][1] = 2 * (h[i] + h[i + 1]);
 674
 675         if (i > num_points - 4)
 676             m[i][2] = 0;
 677         else
 678             m[i][2] = h[i + 1];
 679     }
 680
 681     /* solve the matrix */
 682     if (!solve_tri_diag(m, r, x, num_points - 2)) {
 683         free(h);
 684         free(x);
 685         free(r);
 686         free(m);
 687         free(xp);
 688         free(yp);
 689         int_error(plot->token, "Can't calculate cubic splines");
 690     }
 691     sc[0][2] = 0;
 692     for (i = 1; i <= num_points - 2; i++)
 693         sc[i][2] = x[i - 1];
 694     sc[num_points - 1][2] = 0;
 695
 696     for (i = 0; i <= num_points - 1; i++)
 697         sc[i][0] = yp[i];
 698
 699     for (i = 0; i <= num_points - 2; i++) {
 700         sc[i][1] = (sc[i + 1][0] - sc[i][0]) / h[i]
 701             - h[i] / 3 * (sc[i + 1][2] + 2 * sc[i][2]);
 702         sc[i][3] = (sc[i + 1][2] - sc[i][2]) / 3 / h[i];
 703     }
 704
 705     free(h);
 706     free(x);
 707     free(r);
 708     free(m);
 709     free(xp);
 710     free(yp);
 711
 712     return (sc);
 713 }
 714
 715 static void
 716 do_cubic(
 717     struct curve_points *plot,  /* still containes old plot->points */
 718     spline_coeff *sc,           /* generated by cp_tridiag */
 719     int first_point,            /* where to start in plot->points */
 720     int num_points,             /* to determine end in plot->points */
 721     struct coordinate *dest)    /* where to put the interpolated data */
 722 {
 723     double xdiff, temp, x, y;
 724     double xstart, xend;        /* Endpoints of the sampled x range */
 725     int i, l;
 726     struct coordinate GPHUGE *this_points;
 727
 728     /* min and max in internal (eg logged) co-ordinates. We update
 729      * these, then update the external extrema in user co-ordinates
 730      * at the end.
 731      */
 732     double ixmin, ixmax, iymin, iymax;
 733     double sxmin, sxmax, symin, symax;  /* starting values of above */
 734
 735     x_axis = plot->x_axis;
 736     y_axis = plot->y_axis;
 737
 738     ixmin = sxmin = AXIS_LOG_VALUE(x_axis, X_AXIS.min);
 739     ixmax = sxmax = AXIS_LOG_VALUE(x_axis, X_AXIS.max);
 740     iymin = symin = AXIS_LOG_VALUE(y_axis, Y_AXIS.min);
 741     iymax = symax = AXIS_LOG_VALUE(y_axis, Y_AXIS.max);
 742
 743     this_points = (plot->points) + first_point;
 744
 745     l = 0;
 746
 747     /* HBB 20010727: Sample only across the actual x range, not the
 748      * full range of input data */
 749 #if SAMPLE_CSPLINES_TO_FULL_RANGE
 750     xstart = this_points[0].x;
 751     xend = this_points[num_points - 1].x;
 752 #else
 753     xstart = GPMAX(this_points[0].x, sxmin);
 754     xend = GPMIN(this_points[num_points - 1].x, sxmax);
 755
 756     if (xstart >= xend)
 757         int_error(plot->token,
 758                   "Cannot smooth: no data within fixed xrange!");
 759 #endif
 760     xdiff = (xend - xstart) / (samples_1 - 1);
 761
 762     for (i = 0; i < samples_1; i++) {
 763         x = xstart + i * xdiff;
 764
 765         /* Move forward to the spline interval this point is in */
 766         while ((x >= this_points[l + 1].x) && (l < num_points - 2))
 767             l++;
 768
 769         /* KB 981107: With logarithmic x axis the values were
 770          * converted back to linear scale before calculating the
 771          * coefficients. Use exponential for log x values. */
 772         temp = AXIS_DE_LOG_VALUE(x_axis, x)
 773             - AXIS_DE_LOG_VALUE(x_axis, this_points[l].x);
 774
 775         /* Evaluate cubic spline polynomial */
 776         y = ((sc[l][3] * temp + sc[l][2]) * temp + sc[l][1]) * temp + sc[l][0];
 777
 778         /* With logarithmic y axis, we need to convert from linear to
 779          * log scale now. */
 780         if (Y_AXIS.log) {
 781             if (y > 0.)
 782                 y = AXIS_DO_LOG(y_axis, y);
 783             else
 784                 y = symin - (symax - symin);
 785         }
 786
 787         dest[i].type = INRANGE;
 788         STORE_AND_FIXUP_RANGE(dest[i].x, x, dest[i].type, ixmin, ixmax, X_AXIS.autoscale, NOOP, continue);
 789         STORE_AND_FIXUP_RANGE(dest[i].y, y, dest[i].type, iymin, iymax, Y_AXIS.autoscale, NOOP, NOOP);
 790
 791         dest[i].xlow = dest[i].xhigh = dest[i].x;
 792         dest[i].ylow = dest[i].yhigh = dest[i].y;
 793
 794         dest[i].z = -1;
 795
 796     }
 797
 798     UPDATE_RANGE(ixmax > sxmax, X_AXIS.max, ixmax, x_axis);
 799     UPDATE_RANGE(ixmin < sxmin, X_AXIS.min, ixmin, x_axis);
 800     UPDATE_RANGE(iymax > symax, Y_AXIS.max, iymax, y_axis);
 801     UPDATE_RANGE(iymin < symin, Y_AXIS.min, iymin, y_axis);
 802
 803 }
 804
 805
 806 /*
 807  * do_freq() is like the other smoothers only in that it
 808  * needs to adjust the plot ranges. We don't have to copy
 809  * approximated curves or anything like that.
 810  */
 811
 812 static void
 813 do_freq(
 814     struct curve_points *plot,  /* still contains old plot->points */
 815     int first_point,            /* where to start in plot->points */
 816     int num_points)             /* to determine end in plot->points */
 817 {
 818     double x, y;
 819     int i;
 820     int x_axis = plot->x_axis;
 821     int y_axis = plot->y_axis;
 822     struct coordinate GPHUGE *this;
 823
 824     /* min and max in internal (eg logged) co-ordinates. We update
 825      * these, then update the external extrema in user co-ordinates
 826      * at the end.
 827      */
 828
 829     double ixmin, ixmax, iymin, iymax;
 830     double sxmin, sxmax, symin, symax;  /* starting values of above */
 831
 832     ixmin = sxmin = AXIS_LOG_VALUE(x_axis, X_AXIS.min);
 833     ixmax = sxmax = AXIS_LOG_VALUE(x_axis, X_AXIS.max);
 834     iymin = symin = AXIS_LOG_VALUE(y_axis, Y_AXIS.min);
 835     iymax = symax = AXIS_LOG_VALUE(y_axis, Y_AXIS.max);
 836
 837     this = (plot->points) + first_point;
 838
 839     for (i=0; i<num_points; i++) {
 840
 841         x = this[i].x;
 842         y = this[i].y;
 843
 844         this[i].type = INRANGE;
 845
 846         STORE_AND_FIXUP_RANGE(this[i].x, x, this[i].type, ixmin, ixmax, X_AXIS.autoscale, NOOP, continue);
 847         STORE_AND_FIXUP_RANGE(this[i].y, y, this[i].type, iymin, iymax, Y_AXIS.autoscale, NOOP, NOOP);
 848
 849         this[i].xlow = this[i].xhigh = this[i].x;
 850         this[i].ylow = this[i].yhigh = this[i].y;
 851         this[i].z = -1;
 852     }
 853
 854     UPDATE_RANGE(ixmax > sxmax, X_AXIS.max, ixmax, x_axis);
 855     UPDATE_RANGE(ixmin < sxmin, X_AXIS.min, ixmin, x_axis);
 856     UPDATE_RANGE(iymax > symax, Y_AXIS.max, iymax, y_axis);
 857     UPDATE_RANGE(iymin < symin, Y_AXIS.min, iymin, y_axis);
 858 }
 859
 860
 861 /*
 862  * Frequency plots have don't need new points allocated; we just need
 863  * to adjust the plot ranges. Wedging this into gen_interp() would
 864  * make that code even harder to read.
 865  */
 866
 867 void
 868 gen_interp_frequency(struct curve_points *plot)
 869 {
 870
 871     int i, curves;
 872     int first_point, num_points;
 873
 874     curves = num_curves(plot);
 875
 876     first_point = 0;
 877     for (i = 0; i < curves; i++) {
 878         num_points = next_curve(plot, &first_point);
 879         do_freq(plot, first_point, num_points);
 880         first_point += num_points + 1;
 881     }
 882     return;
 883 }
 884
 885 /*
 886  * This is the main entry point used for everything except frequencies.
 887  * As stated in the header, it is fine, but I'm not too happy with it.
 888  */
 889
 890 void
 891 gen_interp(struct curve_points *plot)
 892 {
 893
 894     spline_coeff *sc;
 895     double *bc;
 896     struct coordinate *new_points;
 897     int i, curves;
 898     int first_point, num_points;
 899
 900     curves = num_curves(plot);
 901     new_points = gp_alloc((samples_1 + 1) * curves * sizeof(struct coordinate),
 902                           "interpolation table");
 903
 904     first_point = 0;
 905     for (i = 0; i < curves; i++) {
 906         num_points = next_curve(plot, &first_point);
 907         switch (plot->plot_smooth) {
 908         case SMOOTH_CSPLINES:
 909             sc = cp_tridiag(plot, first_point, num_points);
 910             do_cubic(plot, sc, first_point, num_points,
 911                      new_points + i * (samples_1 + 1));
 912             free(sc);
 913             break;
 914         case SMOOTH_ACSPLINES:
 915             sc = cp_approx_spline(plot, first_point, num_points);
 916             do_cubic(plot, sc, first_point, num_points,
 917                      new_points + i * (samples_1 + 1));
 918             free(sc);
 919             break;
 920
 921         case SMOOTH_BEZIER:
 922         case SMOOTH_SBEZIER:
 923             bc = cp_binomial(num_points);
 924             do_bezier(plot, bc, first_point, num_points,
 925                       new_points + i * (samples_1 + 1));
 926             free((char *) bc);
 927             break;
 928         default:                /* keep gcc -Wall quiet */
 929             ;
 930         }
 931         new_points[(i + 1) * (samples_1 + 1) - 1].type = UNDEFINED;
 932         first_point += num_points;
 933     }
 934
 935     free(plot->points);
 936     plot->points = new_points;
 937     plot->p_max = curves * (samples_1 + 1);
 938     plot->p_count = plot->p_max - 1;
 939
 940     return;
 941 }
 942
 943 /*
 944  * sort_points
 945  *
 946  * sort data succession for further evaluation by plot_splines, etc.
 947  * This routine is mainly introduced for compilers *NOT* supporting the
 948  * UNIX qsort() routine. You can then easily replace it by more convenient
 949  * stuff for your compiler.
 950  * (MGR 1992)
 951  */
 952
 953 /* HBB 20010720: To avoid undefined behaviour that would be caused by
 954  * casting functions pointers around, changed arguments to what
 955  * qsort() *really* wants */
 956 /* HBB 20010720: removed 'static' to avoid HP-sUX gcc bug */
 957 int
 958 compare_points(SORTFUNC_ARGS arg1, SORTFUNC_ARGS arg2)
 959 {
 960     struct coordinate const *p1 = arg1;
 961     struct coordinate const *p2 = arg2;
 962
 963     if (p1->x > p2->x)
 964         return (1);
 965     if (p1->x < p2->x)
 966         return (-1);
 967     return (0);
 968 }
 969
 970 void
 971 sort_points(struct curve_points *plot)
 972 {
 973     int first_point, num_points;
 974
 975     first_point = 0;
 976     while ((num_points = next_curve(plot, &first_point)) > 0) {
 977         /* Sort this set of points, does qsort handle 1 point correctly? */
 978         /* HBB 20010720: removed casts -- they don't help a thing, but
 979          * may hide problems */
 980         qsort(plot->points + first_point, num_points,
 981               sizeof(struct coordinate), compare_points);
 982         first_point += num_points;
 983     }
 984     return;
 985 }
 986
 987
 988 /*
 989  * cp_implode() if averaging is selected this function computes the new
 990  *              entries and shortens the whole thing to the necessary
 991  *              size
 992  * MGR Addendum
 993  */
 994
 995 void
 996 cp_implode(struct curve_points *cp)
 997 {
 998     int first_point, num_points;
 999     int i, j, k;
1000     double x = 0., y = 0., sux = 0., slx = 0., suy = 0., sly = 0.;
1001     /* int x_axis = cp->x_axis; */
1002     /* int y_axis = cp->y_axis; */
1003     TBOOLEAN all_inrange = FALSE;
1004
1005     x_axis = cp->x_axis;
1006     y_axis = cp->y_axis;
1007     j = 0;
1008     first_point = 0;
1009     while ((num_points = next_curve(cp, &first_point)) > 0) {
1010         k = 0;
1011         for (i = first_point; i < first_point + num_points; i++) {
1012             /* HBB 20020801: don't try to use undefined datapoints */
1013             if (cp->points[i].type == UNDEFINED)
1014                 continue;
1015             if (!k) {
1016                 x = cp->points[i].x;
1017                 y = cp->points[i].y;
1018                 sux = cp->points[i].xhigh;
1019                 slx = cp->points[i].xlow;
1020                 suy = cp->points[i].yhigh;
1021                 sly = cp->points[i].ylow;
1022                 all_inrange = (cp->points[i].type == INRANGE);
1023                 k = 1;
1024             } else if (cp->points[i].x == x) {
1025                 y += cp->points[i].y;
1026                 sux += cp->points[i].xhigh;
1027                 slx += cp->points[i].xlow;
1028                 suy += cp->points[i].yhigh;
1029                 sly += cp->points[i].ylow;
1030                 if (cp->points[i].type != INRANGE)
1031                     all_inrange = FALSE;
1032                 k++;
1033             } else {
1034                 cp->points[j].x = x;
1035                 if ( cp->plot_smooth == SMOOTH_FREQUENCY )
1036                     k = 1;
1037                 cp->points[j].y = y /= (double) k;
1038                 cp->points[j].xhigh = sux / (double) k;
1039                 cp->points[j].xlow = slx / (double) k;
1040                 cp->points[j].yhigh = suy / (double) k;
1041                 cp->points[j].ylow = sly / (double) k;
1042                 /* HBB 20000405: I wanted to use STORE_AND_FIXUP_RANGE
1043                  * here, but won't: it assumes we want to modify the
1044                  * range, and that the range is given in 'input'
1045                  * coordinates.  For logarithmic axes, the overhead
1046                  * would be larger than the possible gain, so write it
1047                  * out explicitly, instead:
1048                  * */
1049                 cp->points[j].type = INRANGE;
1050                 if (! all_inrange) {
1051                     if (X_AXIS.log) {
1052                         if (x <= -VERYLARGE) {
1053                             cp->points[j].type = OUTRANGE;
1054                             goto is_outrange;
1055                         }
1056                         x = AXIS_UNDO_LOG(x_axis, x);
1057                     }
1058                     if (((x < X_AXIS.min) && !(X_AXIS.autoscale & AUTOSCALE_MIN))
1059                         || ((x > X_AXIS.max) && !(X_AXIS.autoscale & AUTOSCALE_MAX))) {
1060                         cp->points[j].type = OUTRANGE;
1061                         goto is_outrange;
1062                     }
1063                     if (Y_AXIS.log) {
1064                         if (y <= -VERYLARGE) {
1065                             cp->points[j].type = OUTRANGE;
1066                             goto is_outrange;
1067                         }
1068                         y = AXIS_UNDO_LOG(y_axis, y);
1069                     }
1070                     if (((y < Y_AXIS.min) && !(Y_AXIS.autoscale & AUTOSCALE_MIN))
1071                         || ((y > Y_AXIS.max) && !(Y_AXIS.autoscale & AUTOSCALE_MAX)))
1072                         cp->points[j].type = OUTRANGE;
1073                 is_outrange:
1074                     ;
1075                 } /* if(! all inrange) */
1076
1077                 j++;            /* next valid entry */
1078                 k = 0;          /* to read */
1079                 i--;            /* from this (-> last after for(;;)) entry */
1080             } /* else (same x position) */
1081         } /* for(points in curve) */
1082
1083         if (k) {
1084             cp->points[j].x = x;
1085             if ( cp->plot_smooth == SMOOTH_FREQUENCY )
1086                 k = 1;
1087             cp->points[j].y = y /= (double) k;
1088             cp->points[j].xhigh = sux / (double) k;
1089             cp->points[j].xlow = slx / (double) k;
1090             cp->points[j].yhigh = suy / (double) k;
1091             cp->points[j].ylow = sly / (double) k;
1092             cp->points[j].type = INRANGE;
1093             if (! all_inrange) {
1094                     if (X_AXIS.log) {
1095                         if (x <= -VERYLARGE) {
1096                             cp->points[j].type = OUTRANGE;
1097                             goto is_outrange2;
1098                         }
1099                         x = AXIS_UNDO_LOG(x_axis, x);
1100                     }
1101                     if (((x < X_AXIS.min) && !(X_AXIS.autoscale & AUTOSCALE_MIN))
1102                         || ((x > X_AXIS.max) && !(X_AXIS.autoscale & AUTOSCALE_MAX))) {
1103                         cp->points[j].type = OUTRANGE;
1104                         goto is_outrange2;
1105                     }
1106                     if (Y_AXIS.log) {
1107                         if (y <= -VERYLARGE) {
1108                             cp->points[j].type = OUTRANGE;
1109                             goto is_outrange2;
1110                         }
1111                         y = AXIS_UNDO_LOG(y_axis, y);
1112                     }
1113                     if (((y < Y_AXIS.min) && !(Y_AXIS.autoscale & AUTOSCALE_MIN))
1114                         || ((y > Y_AXIS.max) && !(Y_AXIS.autoscale & AUTOSCALE_MAX)))
1115                         cp->points[j].type = OUTRANGE;
1116                 is_outrange2:
1117                     ;
1118             }
1119             j++;                /* next valid entry */
1120         }
1121         /* insert invalid point to separate curves */
1122         if (j < cp->p_count) {
1123             cp->points[j].type = UNDEFINED;
1124             j++;
1125         }
1126         first_point += num_points;
1127     }                           /* end while */
1128     cp->p_count = j;
1129     cp_extend(cp, j);
1130 }