Annotation of OpenXM_contrib/gnuplot/interpol.c, Revision 1.1
1.1 ! maekawa 1: #ifndef lint
! 2: static char *RCSid = "$Id: interpol.c,v 1.29 1998/04/14 00:15:45 drd Exp $";
! 3: #endif
! 4:
! 5: /* GNUPLOT - interpol.c */
! 6:
! 7: /*[
! 8: * Copyright 1986 - 1993, 1998 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, is_log_x, base_log_x, xmin, xmax, autoscale_lx
! 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 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 "plot.h"
! 118: #include "setshow.h"
! 119:
! 120:
! 121: /* in order to support multiple axes, and to simplify ranging in
! 122: * parametric plots, we use arrays to store some things. For 2d plots,
! 123: * elements are z=0,y1=1,x1=2,z2=4,y2=5,x2=6 these are given symbolic
! 124: * names in plot.h
! 125: */
! 126:
! 127: extern double min_array[AXIS_ARRAY_SIZE];
! 128: extern double max_array[AXIS_ARRAY_SIZE];
! 129: extern int auto_array[AXIS_ARRAY_SIZE];
! 130: extern TBOOLEAN log_array[AXIS_ARRAY_SIZE];
! 131: extern double base_array[AXIS_ARRAY_SIZE];
! 132: extern double log_base_array[AXIS_ARRAY_SIZE];
! 133:
! 134:
! 135: #define Inc_c_token if (++c_token >= num_tokens) \
! 136: int_error ("Syntax error", c_token);
! 137:
! 138:
! 139: /*
! 140: * IMHO, code is getting too cluttered with repeated chunks of
! 141: * code. Some macros to simplify, I hope.
! 142: *
! 143: * do { } while(0) is comp.lang.c recommendation for complex macros
! 144: * also means that break can be specified as an action, and it will
! 145: *
! 146: */
! 147:
! 148:
! 149: /* store VALUE or log(VALUE) in STORE, set TYPE as appropriate
! 150: * Do OUT_ACTION or UNDEF_ACTION as appropriate
! 151: * adjust range provided type is INRANGE (ie dont adjust y if x is outrange
! 152: * VALUE must not be same as STORE
! 153: */
! 154:
! 155: #define STORE_AND_FIXUP_RANGE(STORE, VALUE, TYPE, MIN, MAX, AUTO, OUT_ACTION, UNDEF_ACTION)\
! 156: do { STORE=VALUE; \
! 157: if (TYPE != INRANGE) break; /* dont set y range if x is outrange, for example */ \
! 158: if ( VALUE<MIN ) { \
! 159: if (AUTO & 1) MIN=VALUE; else { TYPE=OUTRANGE; OUT_ACTION; break; } \
! 160: } \
! 161: if ( VALUE>MAX ) {\
! 162: if (AUTO & 2) MAX=VALUE; else { TYPE=OUTRANGE; OUT_ACTION; } \
! 163: } \
! 164: } while(0)
! 165:
! 166: #define UPDATE_RANGE(TEST,OLD,NEW,AXIS) \
! 167: do { if (TEST) { \
! 168: if (log_array[AXIS]) OLD = pow(base_array[AXIS], NEW); else OLD = NEW; \
! 169: } \
! 170: } while(0)
! 171:
! 172: /* use this instead empty macro arguments to work around NeXT cpp bug */
! 173: /* if this fails on any system, we might use ((void)0) */
! 174:
! 175: #define NOOP /* */
! 176:
! 177: #define spline_coeff_size 4
! 178: typedef double spline_coeff[spline_coeff_size];
! 179: typedef double five_diag[5];
! 180:
! 181: static int next_curve __PROTO((struct curve_points * plot, int *curve_start));
! 182: static int num_curves __PROTO((struct curve_points * plot));
! 183: static double *cp_binomial __PROTO((int points));
! 184: GP_INLINE static double s_pow __PROTO((double base, unsigned int exponent));
! 185: static void eval_bezier __PROTO((struct curve_points * cp, int first_point, int num_points, double sr, coordval * px, coordval * py, double *c));
! 186: static void do_bezier __PROTO((struct curve_points * cp, double *bc, int first_point, int num_points, struct coordinate * dest));
! 187: static int solve_five_diag __PROTO((five_diag m[], double r[], double x[], int n));
! 188: static spline_coeff *cp_approx_spline __PROTO((struct curve_points * plot, int first_point, int num_points));
! 189: static spline_coeff *cp_tridiag __PROTO((struct curve_points * plot, int first_point, int num_points));
! 190: static void do_cubic __PROTO((struct curve_points * plot, spline_coeff * sc, int first_point, int num_points, struct coordinate * dest));
! 191: static int compare_points __PROTO((struct coordinate * p1, struct coordinate * p2));
! 192:
! 193:
! 194: /*
! 195: * position curve_start to index the next non-UNDEFINDED point,
! 196: * start search at initial curve_start,
! 197: * return number of non-UNDEFINDED points from there on,
! 198: * if no more valid points are found, curve_start is set
! 199: * to plot->p_count and 0 is returned
! 200: */
! 201:
! 202: static int next_curve(plot, curve_start)
! 203: struct curve_points *plot;
! 204: int *curve_start;
! 205: {
! 206: int curve_length;
! 207:
! 208: /* Skip undefined points */
! 209: while (*curve_start < plot->p_count
! 210: && plot->points[*curve_start].type == UNDEFINED) {
! 211: (*curve_start)++;
! 212: };
! 213: curve_length = 0;
! 214: /* curve_length is first used as an offset, then the correkt # points */
! 215: while ((*curve_start) + curve_length < plot->p_count
! 216: && plot->points[(*curve_start) + curve_length].type != UNDEFINED) {
! 217: curve_length++;
! 218: };
! 219: return (curve_length);
! 220: }
! 221:
! 222:
! 223: /*
! 224: * determine the number of curves in plot->points, separated by
! 225: * UNDEFINED points
! 226: */
! 227:
! 228: static int num_curves(plot)
! 229: struct curve_points *plot;
! 230: {
! 231: int curves;
! 232: int first_point;
! 233: int num_points;
! 234:
! 235: first_point = 0;
! 236: curves = 0;
! 237: while ((num_points = next_curve(plot, &first_point)) > 0) {
! 238: curves++;
! 239: first_point += num_points;
! 240: }
! 241: return (curves);
! 242: }
! 243:
! 244:
! 245:
! 246: /*
! 247: * build up a cntr_struct list from curve_points
! 248: * this funtion is only used for the alternate entry point to
! 249: * Gershon's code and thus commented out
! 250: ***deleted***
! 251: */
! 252:
! 253:
! 254: /*
! 255: * cp_binomial() computes the binomial coefficients needed for BEZIER stuff
! 256: * and stores them into an array which is hooked to sdat.
! 257: * (MGR 1992)
! 258: */
! 259: static double *cp_binomial(points)
! 260: int points;
! 261: {
! 262: register double *coeff;
! 263: register int n, k;
! 264: int e;
! 265:
! 266: e = points; /* well we're going from k=0 to k=p_count-1 */
! 267: coeff = (double *) gp_alloc(e * sizeof(double), "bezier coefficients");
! 268:
! 269: n = points - 1;
! 270: e = n / 2;
! 271: coeff[0] = 1.0;
! 272:
! 273: for (k = 0; k < e; k++) {
! 274: coeff[k + 1] = coeff[k] * ((double) (n - k)) / ((double) (k + 1));
! 275: }
! 276:
! 277: for (k = n; k >= e; k--)
! 278: coeff[k] = coeff[n - k];
! 279:
! 280: return (coeff);
! 281: }
! 282:
! 283:
! 284: /* This is a subfunction of do_bezier() for BEZIER style computations.
! 285: * It is passed the stepration (STEP/MAXSTEPS) and the addresses of
! 286: * the double values holding the next x and y coordinates.
! 287: * (MGR 1992)
! 288: */
! 289:
! 290: /*
! 291: * well, this routine runs faster with the 68040 striptease FPU
! 292: * and it handles zeroes correctly - there had been some trouble with TC
! 293: */
! 294:
! 295: GP_INLINE static double s_pow(base, exponent)
! 296: double base;
! 297: unsigned int exponent;
! 298: {
! 299: double y;
! 300:
! 301: if (exponent == 0)
! 302: return (1.0);
! 303: if (base == 0.0)
! 304: return (0.0);
! 305:
! 306: /* consider i in binary = abcd
! 307: * x^i = x^(8a+4b+2c+d) = x^8a * x^4b * x^2b * x^d
! 308: * = a?x^2^2^2:1 * b?x^2^2:1 + ...
! 309: * so for each bit in exponent, square x, multiplying it into accumulator
! 310: *
! 311: */
! 312:
! 313: y = 1;
! 314: while (exponent) {
! 315: if (exponent & 1)
! 316: y *= base;
! 317: base *= base;
! 318: /* if exponent was signed, this could be trouble ! */
! 319: exponent >>= 1;
! 320: }
! 321:
! 322: return (y);
! 323: }
! 324:
! 325: static void eval_bezier(cp, first_point, num_points, sr, px, py, c)
! 326: struct curve_points *cp;
! 327: int first_point; /* where to start in plot->points (to find x-range) */
! 328: int num_points; /* to determine end in plot->points */
! 329: double sr;
! 330: coordval *px;
! 331: coordval *py;
! 332: double *c;
! 333: {
! 334: unsigned int n = num_points - 1;
! 335: /* HBB 980308: added 'GPHUGE' tag for DOS */
! 336: struct coordinate GPHUGE *this_points;
! 337:
! 338: this_points = (cp->points) + first_point;
! 339:
! 340: if (sr == 0.0) {
! 341: *px = this_points[0].x;
! 342: *py = this_points[0].y;
! 343: } else if (sr == 1.0) {
! 344: *px = this_points[n].x;
! 345: *py = this_points[n].y;
! 346: } else {
! 347: unsigned int i;
! 348: double lx = 0.0, ly = 0.0;
! 349: double sr_to_the_i = 1.0;
! 350: double dsr_to_the_di = s_pow(1 - sr, n);
! 351: double reciprocal_dsr = 1.0 / (1 - sr);
! 352:
! 353: for (i = 0; i <= n; i++) {
! 354: double u = c[i] * sr_to_the_i * s_pow(1 - sr, n - i);
! 355: lx += this_points[i].x * u;
! 356: ly += this_points[i].y * u;
! 357: sr_to_the_i *= sr;
! 358: dsr_to_the_di *= reciprocal_dsr;
! 359: }
! 360:
! 361: *px = lx;
! 362: *py = ly;
! 363: }
! 364: }
! 365:
! 366: /*
! 367: * generate a new set of coordinates representing the bezier curve and
! 368: * set it to the plot
! 369: */
! 370:
! 371: static void do_bezier(cp, bc, first_point, num_points, dest)
! 372: struct curve_points *cp;
! 373: double *bc;
! 374: int first_point; /* where to start in plot->points */
! 375: int num_points; /* to determine end in plot->points */
! 376: struct coordinate *dest; /* where to put the interpolated data */
! 377: {
! 378: int i;
! 379: coordval x, y;
! 380: int xaxis = cp->x_axis;
! 381: int yaxis = cp->y_axis;
! 382:
! 383: /* min and max in internal (eg logged) co-ordinates. We update
! 384: * these, then update the external extrema in user co-ordinates
! 385: * at the end.
! 386: */
! 387:
! 388: double ixmin, ixmax, iymin, iymax;
! 389: double sxmin, sxmax, symin, symax; /* starting values of above */
! 390:
! 391: if (log_array[xaxis]) {
! 392: ixmin = sxmin = log(min_array[xaxis]) / log_base_array[xaxis];
! 393: ixmax = sxmax = log(max_array[xaxis]) / log_base_array[xaxis];
! 394: } else {
! 395: ixmin = sxmin = min_array[xaxis];
! 396: ixmax = sxmax = max_array[xaxis];
! 397: }
! 398:
! 399: if (log_array[yaxis]) {
! 400: iymin = symin = log(min_array[yaxis]) / log_base_array[yaxis];
! 401: iymax = symax = log(max_array[yaxis]) / log_base_array[yaxis];
! 402: } else {
! 403: iymin = symin = min_array[yaxis];
! 404: iymax = symax = max_array[yaxis];
! 405: }
! 406:
! 407: for (i = 0; i < samples; i++) {
! 408: eval_bezier(cp, first_point, num_points, (double) i / (double) (samples - 1), &x, &y, bc);
! 409:
! 410: /* now we have to store the points and adjust the ranges */
! 411:
! 412: dest[i].type = INRANGE;
! 413: STORE_AND_FIXUP_RANGE(dest[i].x, x, dest[i].type, ixmin, ixmax, auto_array[xaxis], NOOP, continue);
! 414: STORE_AND_FIXUP_RANGE(dest[i].y, y, dest[i].type, iymin, iymax, auto_array[yaxis], NOOP, NOOP);
! 415:
! 416: dest[i].xlow = dest[i].xhigh = dest[i].x;
! 417: dest[i].ylow = dest[i].yhigh = dest[i].y;
! 418:
! 419: dest[i].z = -1;
! 420: }
! 421:
! 422: UPDATE_RANGE(ixmax > sxmax, max_array[xaxis], ixmax, xaxis);
! 423: UPDATE_RANGE(ixmin < sxmin, min_array[xaxis], ixmin, xaxis);
! 424: UPDATE_RANGE(iymax > symax, max_array[yaxis], iymax, yaxis);
! 425: UPDATE_RANGE(iymin < symin, min_array[yaxis], iymin, yaxis);
! 426: }
! 427:
! 428: /*
! 429: * call contouring routines -- main entry
! 430: */
! 431:
! 432: /*
! 433: * it should be like this, but it doesn't run. If you find out why,
! 434: * contact me: mgr@asgard.bo.open.de or Lars Hanke 2:243/4802.22@fidonet
! 435: *
! 436: * Well, all this had originally been inside contour.c, so maybe links
! 437: * to functions and of contour.c are broken.
! 438: * ***deleted***
! 439: * end of unused entry point to Gershon's code
! 440: *
! 441: */
! 442:
! 443: /*
! 444: * Solve five diagonal linear system equation. The five diagonal matrix is
! 445: * defined via matrix M, right side is r, and solution X i.e. M * X = R.
! 446: * Size of system given in n. Return TRUE if solution exist.
! 447: * G. Engeln-Muellges/ F.Reutter:
! 448: * "Formelsammlung zur Numerischen Mathematik mit Standard-FORTRAN-Programmen"
! 449: * ISBN 3-411-01677-9
! 450: *
! 451: * / m02 m03 m04 0 0 0 0 . . . \ / x0 \ / r0 \
! 452: * I m11 m12 m13 m14 0 0 0 . . . I I x1 I I r1 I
! 453: * I m20 m21 m22 m23 m24 0 0 . . . I * I x2 I = I r2 I
! 454: * I 0 m30 m31 m32 m33 m34 0 . . . I I x3 I I r3 I
! 455: * . . . . . . . . . . . .
! 456: * \ m(n-3)0 m(n-2)1 m(n-1)2 / \x(n-1)/ \r(n-1)/
! 457: *
! 458: */
! 459: static int solve_five_diag(m, r, x, n)
! 460: five_diag m[];
! 461: double r[], x[];
! 462: int n;
! 463: {
! 464: int i;
! 465: five_diag *hv;
! 466:
! 467: hv = (five_diag *) gp_alloc((n + 1) * sizeof(five_diag), "five_diag help vars");
! 468:
! 469: hv[0][0] = m[0][2];
! 470: if (hv[0][0] == 0) {
! 471: free(hv);
! 472: return FALSE;
! 473: }
! 474: hv[0][1] = m[0][3] / hv[0][0];
! 475: hv[0][2] = m[0][4] / hv[0][0];
! 476:
! 477: hv[1][3] = m[1][1];
! 478: hv[1][0] = m[1][2] - hv[1][3] * hv[0][1];
! 479: if (hv[1][0] == 0) {
! 480: free(hv);
! 481: return FALSE;
! 482: }
! 483: hv[1][1] = (m[1][3] - hv[1][3] * hv[0][2]) / hv[1][0];
! 484: hv[1][2] = m[1][4] / hv[1][0];
! 485:
! 486: for (i = 2; i <= n - 1; i++) {
! 487: hv[i][3] = m[i][1] - m[i][0] * hv[i - 2][1];
! 488: hv[i][0] = m[i][2] - m[i][0] * hv[i - 2][2] - hv[i][3] * hv[i - 1][1];
! 489: if (hv[i][0] == 0) {
! 490: free(hv);
! 491: return FALSE;
! 492: }
! 493: hv[i][1] = (m[i][3] - hv[i][3] * hv[i - 1][2]) / hv[i][0];
! 494: hv[i][2] = m[i][4] / hv[i][0];
! 495: }
! 496:
! 497: hv[0][4] = 0;
! 498: hv[1][4] = r[0] / hv[0][0];
! 499: for (i = 1; i <= n - 1; i++) {
! 500: hv[i + 1][4] = (r[i] - m[i][0] * hv[i - 1][4] - hv[i][3] * hv[i][4]) / hv[i][0];
! 501: }
! 502:
! 503: x[n - 1] = hv[n][4];
! 504: x[n - 2] = hv[n - 1][4] - hv[n - 2][1] * x[n - 1];
! 505: for (i = n - 3; i >= 0; i--)
! 506: x[i] = hv[i + 1][4] - hv[i][1] * x[i + 1] - hv[i][2] * x[i + 2];
! 507:
! 508: free(hv);
! 509: return TRUE;
! 510: }
! 511:
! 512:
! 513: /*
! 514: * Calculation of approximation cubic splines
! 515: * Input: x[i], y[i], weights z[i]
! 516: *
! 517: * Returns matrix of spline coefficients
! 518: */
! 519: static spline_coeff *cp_approx_spline(plot, first_point, num_points)
! 520: struct curve_points *plot;
! 521: int first_point; /* where to start in plot->points */
! 522: int num_points; /* to determine end in plot->points */
! 523: {
! 524: spline_coeff *sc;
! 525: five_diag *m;
! 526: int xaxis = plot->x_axis;
! 527: int yaxis = plot->y_axis;
! 528: double *r, *x, *h, *xp, *yp;
! 529: /* HBB 980308: added 'GPHUGE' tag */
! 530: struct coordinate GPHUGE *this_points;
! 531: int i;
! 532:
! 533: sc = (spline_coeff *) gp_alloc((num_points) * sizeof(spline_coeff),
! 534: "spline matrix");
! 535:
! 536: if (num_points < 4)
! 537: int_error("Can't calculate approximation splines, need at least 4 points", NO_CARET);
! 538:
! 539: this_points = (plot->points) + first_point;
! 540:
! 541: for (i = 0; i <= num_points - 1; i++)
! 542: if (this_points[i].z <= 0)
! 543: int_error("Can't calculate approximation splines, all weights have to be > 0", NO_CARET);
! 544:
! 545: m = (five_diag *) gp_alloc((num_points - 2) * sizeof(five_diag), "spline help matrix");
! 546:
! 547: r = (double *) gp_alloc((num_points - 2) * sizeof(double), "spline right side");
! 548: x = (double *) gp_alloc((num_points - 2) * sizeof(double), "spline solution vector");
! 549: h = (double *) gp_alloc((num_points - 1) * sizeof(double), "spline help vector");
! 550:
! 551: xp = (double *) gp_alloc((num_points) * sizeof(double), "x pos");
! 552: yp = (double *) gp_alloc((num_points) * sizeof(double), "y pos");
! 553:
! 554: /* KB 981107: With logarithmic axis first convert back to linear scale */
! 555:
! 556: if (log_array[xaxis]) {
! 557: for (i = 0; i <= num_points - 1; i++)
! 558: xp[i] = exp(this_points[i].x * log_base_array[xaxis]);
! 559: } else {
! 560: for (i = 0; i <= num_points - 1; i++)
! 561: xp[i] = this_points[i].x;
! 562: }
! 563: if (log_array[yaxis]) {
! 564: for (i = 0; i <= num_points - 1; i++)
! 565: yp[i] = exp(this_points[i].y * log_base_array[yaxis]);
! 566: } else {
! 567: for (i = 0; i <= num_points - 1; i++)
! 568: yp[i] = this_points[i].y;
! 569: }
! 570:
! 571: for (i = 0; i <= num_points - 2; i++)
! 572: h[i] = xp[i + 1] - xp[i];
! 573:
! 574: /* set up the matrix and the vector */
! 575:
! 576: for (i = 0; i <= num_points - 3; i++) {
! 577: r[i] = 3 * ((yp[i + 2] - yp[i + 1]) / h[i + 1]
! 578: - (yp[i + 1] - yp[i]) / h[i]);
! 579:
! 580: if (i < 2)
! 581: m[i][0] = 0;
! 582: else
! 583: m[i][0] = 6 / this_points[i].z / h[i - 1] / h[i];
! 584:
! 585: if (i < 1)
! 586: m[i][1] = 0;
! 587: else
! 588: m[i][1] = h[i] - 6 / this_points[i].z / h[i] * (1 / h[i - 1] + 1 / h[i])
! 589: - 6 / this_points[i + 1].z / h[i] * (1 / h[i] + 1 / h[i + 1]);
! 590:
! 591: m[i][2] = 2 * (h[i] + h[i + 1])
! 592: + 6 / this_points[i].z / h[i] / h[i]
! 593: + 6 / this_points[i + 1].z * (1 / h[i] + 1 / h[i + 1]) * (1 / h[i] + 1 / h[i + 1])
! 594: + 6 / this_points[i + 2].z / h[i + 1] / h[i + 1];
! 595:
! 596: if (i > num_points - 4)
! 597: m[i][3] = 0;
! 598: else
! 599: m[i][3] = h[i + 1] - 6 / this_points[i + 1].z / h[i + 1] * (1 / h[i] + 1 / h[i + 1])
! 600: - 6 / this_points[i + 2].z / h[i + 1] * (1 / h[i + 1] + 1 / h[i + 2]);
! 601:
! 602: if (i > num_points - 5)
! 603: m[i][4] = 0;
! 604: else
! 605: m[i][4] = 6 / this_points[i + 2].z / h[i + 1] / h[i + 2];
! 606: }
! 607:
! 608: /* solve the matrix */
! 609: if (!solve_five_diag(m, r, x, num_points - 2)) {
! 610: free(h);
! 611: free(x);
! 612: free(r);
! 613: free(m);
! 614: free(xp);
! 615: free(yp);
! 616: int_error("Can't calculate approximation splines", NO_CARET);
! 617: }
! 618: sc[0][2] = 0;
! 619: for (i = 1; i <= num_points - 2; i++)
! 620: sc[i][2] = x[i - 1];
! 621: sc[num_points - 1][2] = 0;
! 622:
! 623: sc[0][0] = yp[0] + 2 / this_points[0].z / h[0] * (sc[0][2] - sc[1][2]);
! 624: for (i = 1; i <= num_points - 2; i++)
! 625: sc[i][0] = yp[i] - 2 / this_points[i].z *
! 626: (sc[i - 1][2] / h[i - 1]
! 627: - sc[i][2] * (1 / h[i - 1] + 1 / h[i])
! 628: + sc[i + 1][2] / h[i]);
! 629: sc[num_points - 1][0] = yp[num_points - 1]
! 630: - 2 / this_points[num_points - 1].z / h[num_points - 2]
! 631: * (sc[num_points - 2][2] - sc[num_points - 1][2]);
! 632:
! 633: for (i = 0; i <= num_points - 2; i++) {
! 634: sc[i][1] = (sc[i + 1][0] - sc[i][0]) / h[i]
! 635: - h[i] / 3 * (sc[i + 1][2] + 2 * sc[i][2]);
! 636: sc[i][3] = (sc[i + 1][2] - sc[i][2]) / 3 / h[i];
! 637: }
! 638:
! 639: free(h);
! 640: free(x);
! 641: free(r);
! 642: free(m);
! 643: free(xp);
! 644: free(yp);
! 645:
! 646: return (sc);
! 647: }
! 648:
! 649:
! 650: /*
! 651: * Calculation of cubic splines
! 652: *
! 653: * This can be treated as a special case of approximation cubic splines, with
! 654: * all weights -> infinity.
! 655: *
! 656: * Returns matrix of spline coefficients
! 657: */
! 658: static spline_coeff *cp_tridiag(plot, first_point, num_points)
! 659: struct curve_points *plot;
! 660: int first_point, num_points;
! 661: {
! 662: spline_coeff *sc;
! 663: tri_diag *m;
! 664: int xaxis = plot->x_axis;
! 665: int yaxis = plot->y_axis;
! 666: double *r, *x, *h, *xp, *yp;
! 667: /* HBB 980308: added 'GPHUGE' tag */
! 668: struct coordinate GPHUGE *this_points;
! 669: int i;
! 670:
! 671: if (num_points < 3)
! 672: int_error("Can't calculate splines, need at least 3 points", NO_CARET);
! 673:
! 674: this_points = (plot->points) + first_point;
! 675:
! 676: sc = (spline_coeff *) gp_alloc((num_points) * sizeof(spline_coeff), "spline matrix");
! 677: m = (tri_diag *) gp_alloc((num_points - 2) * sizeof(tri_diag), "spline help matrix");
! 678:
! 679: r = (double *) gp_alloc((num_points - 2) * sizeof(double), "spline right side");
! 680: x = (double *) gp_alloc((num_points - 2) * sizeof(double), "spline solution vector");
! 681: h = (double *) gp_alloc((num_points - 1) * sizeof(double), "spline help vector");
! 682:
! 683: xp = (double *) gp_alloc((num_points) * sizeof(double), "x pos");
! 684: yp = (double *) gp_alloc((num_points) * sizeof(double), "y pos");
! 685:
! 686: /* KB 981107: With logarithmic axis first convert back to linear scale */
! 687:
! 688: if (log_array[xaxis]) {
! 689: for (i = 0; i <= num_points - 1; i++)
! 690: xp[i] = exp(this_points[i].x * log_base_array[xaxis]);
! 691: } else {
! 692: for (i = 0; i <= num_points - 1; i++)
! 693: xp[i] = this_points[i].x;
! 694: }
! 695: if (log_array[yaxis]) {
! 696: for (i = 0; i <= num_points - 1; i++)
! 697: yp[i] = exp(this_points[i].y * log_base_array[yaxis]);
! 698: } else {
! 699: for (i = 0; i <= num_points - 1; i++)
! 700: yp[i] = this_points[i].y;
! 701: }
! 702:
! 703: for (i = 0; i <= num_points - 2; i++)
! 704: h[i] = xp[i + 1] - xp[i];
! 705:
! 706: /* set up the matrix and the vector */
! 707:
! 708: for (i = 0; i <= num_points - 3; i++) {
! 709: r[i] = 3 * ((yp[i + 2] - yp[i + 1]) / h[i + 1]
! 710: - (yp[i + 1] - yp[i]) / h[i]);
! 711:
! 712: if (i < 1)
! 713: m[i][0] = 0;
! 714: else
! 715: m[i][0] = h[i];
! 716:
! 717: m[i][1] = 2 * (h[i] + h[i + 1]);
! 718:
! 719: if (i > num_points - 4)
! 720: m[i][2] = 0;
! 721: else
! 722: m[i][2] = h[i + 1];
! 723: }
! 724:
! 725: /* solve the matrix */
! 726: if (!solve_tri_diag(m, r, x, num_points - 2)) {
! 727: free(h);
! 728: free(x);
! 729: free(r);
! 730: free(m);
! 731: free(xp);
! 732: free(yp);
! 733: int_error("Can't calculate cubic splines", NO_CARET);
! 734: }
! 735: sc[0][2] = 0;
! 736: for (i = 1; i <= num_points - 2; i++)
! 737: sc[i][2] = x[i - 1];
! 738: sc[num_points - 1][2] = 0;
! 739:
! 740: for (i = 0; i <= num_points - 1; i++)
! 741: sc[i][0] = yp[i];
! 742:
! 743: for (i = 0; i <= num_points - 2; i++) {
! 744: sc[i][1] = (sc[i + 1][0] - sc[i][0]) / h[i]
! 745: - h[i] / 3 * (sc[i + 1][2] + 2 * sc[i][2]);
! 746: sc[i][3] = (sc[i + 1][2] - sc[i][2]) / 3 / h[i];
! 747: }
! 748:
! 749: free(h);
! 750: free(x);
! 751: free(r);
! 752: free(m);
! 753: free(xp);
! 754: free(yp);
! 755:
! 756: return (sc);
! 757: }
! 758:
! 759: static void do_cubic(plot, sc, first_point, num_points, dest)
! 760: struct curve_points *plot; /* still containes old plot->points */
! 761: spline_coeff *sc; /* generated by cp_tridiag */
! 762: int first_point; /* where to start in plot->points */
! 763: int num_points; /* to determine end in plot->points */
! 764: struct coordinate *dest; /* where to put the interpolated data */
! 765: {
! 766: double xdiff, temp, x, y;
! 767: int i, l;
! 768: int xaxis = plot->x_axis;
! 769: int yaxis = plot->y_axis;
! 770: /* HBB 980308: added 'GPHUGE' tag */
! 771: struct coordinate GPHUGE *this_points;
! 772:
! 773: /* min and max in internal (eg logged) co-ordinates. We update
! 774: * these, then update the external extrema in user co-ordinates
! 775: * at the end.
! 776: */
! 777:
! 778: double ixmin, ixmax, iymin, iymax;
! 779: double sxmin, sxmax, symin, symax; /* starting values of above */
! 780:
! 781: if (log_array[xaxis]) {
! 782: ixmin = sxmin = log(min_array[xaxis]) / log_base_array[xaxis];
! 783: ixmax = sxmax = log(max_array[xaxis]) / log_base_array[xaxis];
! 784: } else {
! 785: ixmin = sxmin = min_array[xaxis];
! 786: ixmax = sxmax = max_array[xaxis];
! 787: }
! 788:
! 789: if (log_array[yaxis]) {
! 790: iymin = symin = log(min_array[yaxis]) / log_base_array[yaxis];
! 791: iymax = symax = log(max_array[yaxis]) / log_base_array[yaxis];
! 792: } else {
! 793: iymin = symin = min_array[yaxis];
! 794: iymax = symax = max_array[yaxis];
! 795: }
! 796:
! 797:
! 798: this_points = (plot->points) + first_point;
! 799:
! 800: l = 0;
! 801:
! 802: xdiff = (this_points[num_points - 1].x - this_points[0].x) / (samples - 1);
! 803:
! 804: for (i = 0; i < samples; i++) {
! 805: x = this_points[0].x + i * xdiff;
! 806:
! 807: while ((x >= this_points[l + 1].x) && (l < num_points - 2))
! 808: l++;
! 809:
! 810: /* KB 981107: With logarithmic x axis the values were converted back to linear */
! 811: /* scale before calculating the coefficients. Use exponential for log x values. */
! 812:
! 813: if (log_array[xaxis]) {
! 814: temp = exp(x * log_base_array[xaxis]) - exp(this_points[l].x * log_base_array[xaxis]);
! 815: y = ((sc[l][3] * temp + sc[l][2]) * temp + sc[l][1]) * temp + sc[l][0];
! 816: } else {
! 817: temp = x - this_points[l].x;
! 818: y = ((sc[l][3] * temp + sc[l][2]) * temp + sc[l][1]) * temp + sc[l][0];
! 819: }
! 820: /* With logarithmic y axis, we need to convert from linear to log scale now. */
! 821: if (log_array[yaxis]) {
! 822: if (y > 0.)
! 823: y = log(y) / log_base_array[yaxis];
! 824: else
! 825: y = symin - (symax - symin);
! 826: }
! 827: dest[i].type = INRANGE;
! 828: STORE_AND_FIXUP_RANGE(dest[i].x, x, dest[i].type, ixmin, ixmax, auto_array[xaxis], NOOP, continue);
! 829: STORE_AND_FIXUP_RANGE(dest[i].y, y, dest[i].type, iymin, iymax, auto_array[yaxis], NOOP, NOOP);
! 830:
! 831: dest[i].xlow = dest[i].xhigh = dest[i].x;
! 832: dest[i].ylow = dest[i].yhigh = dest[i].y;
! 833:
! 834: dest[i].z = -1;
! 835:
! 836: }
! 837:
! 838: UPDATE_RANGE(ixmax > sxmax, max_array[xaxis], ixmax, xaxis);
! 839: UPDATE_RANGE(ixmin < sxmin, min_array[xaxis], ixmin, xaxis);
! 840: UPDATE_RANGE(iymax > symax, max_array[yaxis], iymax, yaxis);
! 841: UPDATE_RANGE(iymin < symin, min_array[yaxis], iymin, yaxis);
! 842:
! 843: }
! 844:
! 845: /*
! 846: * This is the main entry point used. As stated in the header, it is fine,
! 847: * but I'm not too happy with it.
! 848: */
! 849:
! 850: void gen_interp(plot)
! 851: struct curve_points *plot;
! 852: {
! 853:
! 854: spline_coeff *sc;
! 855: double *bc;
! 856: struct coordinate *new_points;
! 857: int i, curves;
! 858: int first_point, num_points;
! 859:
! 860: curves = num_curves(plot);
! 861: new_points = (struct coordinate *) gp_alloc((samples + 1) * curves * sizeof(struct coordinate), "interpolation table");
! 862:
! 863: first_point = 0;
! 864: for (i = 0; i < curves; i++) {
! 865: num_points = next_curve(plot, &first_point);
! 866: switch (plot->plot_smooth) {
! 867: case CSPLINES:
! 868: sc = cp_tridiag(plot, first_point, num_points);
! 869: do_cubic(plot, sc, first_point, num_points,
! 870: new_points + i * (samples + 1));
! 871: free(sc);
! 872: break;
! 873: case ACSPLINES:
! 874: sc = cp_approx_spline(plot, first_point, num_points);
! 875: do_cubic(plot, sc, first_point, num_points,
! 876: new_points + i * (samples + 1));
! 877: free(sc);
! 878: break;
! 879:
! 880: case BEZIER:
! 881: case SBEZIER:
! 882: bc = cp_binomial(num_points);
! 883: do_bezier(plot, bc, first_point, num_points,
! 884: new_points + i * (samples + 1));
! 885: free((char *) bc);
! 886: break;
! 887: default: /* keep gcc -Wall quiet */
! 888: ;
! 889: }
! 890: new_points[(i + 1) * (samples + 1) - 1].type = UNDEFINED;
! 891: first_point += num_points;
! 892: }
! 893:
! 894: free(plot->points);
! 895: plot->points = new_points;
! 896: plot->p_max = curves * (samples + 1);
! 897: plot->p_count = plot->p_max - 1;
! 898:
! 899: return;
! 900: }
! 901:
! 902: /*
! 903: * sort_points
! 904: *
! 905: * sort data succession for further evaluation by plot_splines, etc.
! 906: * This routine is mainly introduced for compilers *NOT* supporting the
! 907: * UNIX qsort() routine. You can then easily replace it by more convenient
! 908: * stuff for your compiler.
! 909: * (MGR 1992)
! 910: */
! 911:
! 912: static int compare_points(p1, p2)
! 913: struct coordinate *p1;
! 914: struct coordinate *p2;
! 915: {
! 916: if (p1->x > p2->x)
! 917: return (1);
! 918: if (p1->x < p2->x)
! 919: return (-1);
! 920: return (0);
! 921: }
! 922:
! 923: void sort_points(plot)
! 924: struct curve_points *plot;
! 925: {
! 926: int first_point, num_points;
! 927:
! 928: first_point = 0;
! 929: while ((num_points = next_curve(plot, &first_point)) > 0) {
! 930: /* Sort this set of points, does qsort handle 1 point correctly? */
! 931: qsort((char *) (plot->points + first_point), num_points,
! 932: sizeof(struct coordinate), (sortfunc) compare_points);
! 933: first_point += num_points;
! 934: }
! 935: return;
! 936: }
! 937:
! 938:
! 939: /*
! 940: * cp_implode() if averaging is selected this function computes the new
! 941: * entries and shortens the whole thing to the necessary
! 942: * size
! 943: * MGR Addendum
! 944: */
! 945:
! 946: void cp_implode(cp)
! 947: struct curve_points *cp;
! 948: {
! 949: int first_point, num_points;
! 950: int i, j, k;
! 951: double x, y, sux, slx, suy, sly;
! 952: enum coord_type dot;
! 953:
! 954:
! 955: j = 0;
! 956: first_point = 0;
! 957: while ((num_points = next_curve(cp, &first_point)) > 0) {
! 958: k = 0;
! 959: for (i = first_point; i < first_point + num_points; i++) {
! 960: if (!k) {
! 961: x = cp->points[i].x;
! 962: y = cp->points[i].y;
! 963: sux = cp->points[i].xhigh;
! 964: slx = cp->points[i].xlow;
! 965: suy = cp->points[i].yhigh;
! 966: sly = cp->points[i].ylow;
! 967: dot = INRANGE;
! 968: if (cp->points[i].type != INRANGE)
! 969: dot = UNDEFINED; /* This means somthing other than usual *//* just signal to check if INRANGE */
! 970: k = 1;
! 971: } else if (cp->points[i].x == x) {
! 972: y += cp->points[i].y;
! 973: sux += cp->points[i].xhigh;
! 974: slx += cp->points[i].xlow;
! 975: suy += cp->points[i].yhigh;
! 976: sly += cp->points[i].ylow;
! 977: if (cp->points[i].type != INRANGE)
! 978: dot = UNDEFINED;
! 979: k++;
! 980: } else {
! 981: cp->points[j].x = x;
! 982: cp->points[j].y = y / (double) k;
! 983: cp->points[j].xhigh = sux / (double) k;
! 984: cp->points[j].xlow = slx / (double) k;
! 985: cp->points[j].yhigh = suy / (double) k;
! 986: cp->points[j].ylow = sly / (double) k;
! 987: cp->points[j].type = INRANGE;
! 988: if (dot != INRANGE) {
! 989: if ((cp->points[j].x > xmax) || (cp->points[j].x < xmin))
! 990: cp->points[j].type = OUTRANGE;
! 991: else if (!autoscale_ly) {
! 992: if ((cp->points[j].y > ymax) || (cp->points[j].y < ymin))
! 993: cp->points[j].type = OUTRANGE;
! 994: } else
! 995: cp->points[j].type = INRANGE;
! 996: }
! 997: j++; /* next valid entry */
! 998: k = 0; /* to read */
! 999: i--; /* from this (-> last after for(;;)) entry */
! 1000: }
! 1001: }
! 1002: if (k) {
! 1003: cp->points[j].x = x;
! 1004: cp->points[j].y = y / (double) k;
! 1005: cp->points[j].xhigh = sux / (double) k;
! 1006: cp->points[j].xlow = slx / (double) k;
! 1007: cp->points[j].yhigh = suy / (double) k;
! 1008: cp->points[j].ylow = sly / (double) k;
! 1009: cp->points[j].type = INRANGE;
! 1010: if (dot != INRANGE) {
! 1011: if ((cp->points[j].x > xmax) || (cp->points[j].x < xmin))
! 1012: cp->points[j].type = OUTRANGE;
! 1013: else if (!autoscale_ly) {
! 1014: if ((cp->points[j].y > ymax) || (cp->points[j].y < ymin))
! 1015: cp->points[j].type = OUTRANGE;
! 1016: } else
! 1017: cp->points[j].type = INRANGE;
! 1018: }
! 1019: j++; /* next valid entry */
! 1020: }
! 1021: /* insert invalid point to separate curves */
! 1022: if (j < cp->p_count) {
! 1023: cp->points[j].type = UNDEFINED;
! 1024: j++;
! 1025: }
! 1026: first_point += num_points;
! 1027: } /* end while */
! 1028: cp->p_count = j;
! 1029: cp_extend(cp, j);
! 1030: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to interpol.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gnuplot