home *** CD-ROM | disk | FTP | other *** search
/ System Booster / System Booster.iso / Archives / GNU / GNUPLOTsrc.lha / interpol.c < prev    next >
Encoding:
C/C++ Source or Header  |  1996-01-22  |  24.5 KB  |  909 lines
  1. #ifndef lint
  2. static char *RCSid="$Id: interpol.c,v 1.17 1995/12/14 21:08:56 drd Exp $";
  3. #endif
  4.  
  5. /*
  6.  * C-Source file identification Header
  7.  *
  8.  * This file belongs to a project which is:
  9.  *
  10.  * done 1993 by MGR-Software, Asgard  (Lars Hanke)
  11.  * written by Lars Hanke
  12.  *
  13.  * Contact me via:
  14.  *
  15.  *  InterNet: mgr@asgard.bo.open.de
  16.  *      FIDO: Lars Hanke @ 2:243/4802.22   (as long as they keep addresses)
  17.  *
  18.  **************************************************************************
  19.  *
  20.  *   Project: gnuplot
  21.  *    Module:
  22.  *      File: interpol.c
  23.  *
  24.  *   Revisor: Lars Hanke
  25.  *   Revised: 26/09/93
  26.  *  Revision: 1.0
  27.  *
  28.  **************************************************************************
  29.  *
  30.  * LEGAL
  31.  *  This module is part of gnuplot and distributed under whatever terms
  32.  *  gnuplot is or will be published, unless exclusive rights are claimed.
  33.  *
  34.  * DESCRIPTION
  35.  *  Supplies 2-D data interpolation and approximation routines
  36.  *
  37.  * IMPORTS
  38.  *  plot.h
  39.  *    - cp_extend()
  40.  *    - structs: curve_points, coordval, coordinate
  41.  *
  42.  *  setshow.h
  43.  *    - samples, is_log_x, base_log_x, xmin, xmax, autoscale_lx
  44.  *    - plottypes
  45.  *
  46.  *  proto.h
  47.  *    - solve_tri_diag()
  48.  *    - typedef tri_diag
  49.  *
  50.  * EXPORTS
  51.  *  gen_interp()
  52.  *  sort_points()
  53.  *  cp_implode()
  54.  *
  55.  * BUGS and TODO
  56.  *  I would really have liked to use Gershon Elbers contouring code for
  57.  *  all the stuff done here, but I failed. So I used my own code.
  58.  *  If somebody is able to consolidate Gershon's code for this purpose
  59.  *  a lot of gnuplot users would be very happy - due to memory problems.
  60.  *
  61.  **************************************************************************
  62.  *
  63.  * HISTORY
  64.  * Changes:
  65.  *  Nov 24, 1995  Markus Schuh (M.Schuh@meteo.uni-koeln.de):
  66.  *      changed the algorithm for csplines
  67.  *      added algorithm for approximation csplines
  68.  *      copied point storage and range fix from plot2d.c
  69.  *
  70.  *  Dec 12, 1995 David Denholm
  71.  *      oops - at the time this is called, stored co-ords are
  72.  *      internal (ie maybe log of data) but min/max are in
  73.  *      user co-ordinates. 
  74.  *      Work with min and max of internal co-ords, and
  75.  *      check at the end whether external min and max need to
  76.  *      be increased. (since samples is typically 100 ; we
  77.  *      dont want to take more logs than necessary)
  78.  *      Also, need to take into account which axes are active
  79.  *
  80.  */
  81.  
  82. #include <math.h>
  83. #include "plot.h"
  84. #include "setshow.h"
  85.  
  86.  
  87. /* in order to support multiple axes, and to
  88.  * simplify ranging in parametric plots, we use
  89.  * arrays to store some things. For 2d plots,
  90.  * elements are  z=0,y1=1,x1=2,z2=4,y2=5,x2=6
  91.  * these are given symbolic names in plot.h
  92.  */
  93.  
  94. extern double          min_array[AXIS_ARRAY_SIZE], max_array[AXIS_ARRAY_SIZE];
  95. extern int             auto_array[AXIS_ARRAY_SIZE];
  96. extern TBOOLEAN        log_array[AXIS_ARRAY_SIZE];
  97. extern double          base_array[AXIS_ARRAY_SIZE];
  98. extern double          log_base_array[AXIS_ARRAY_SIZE];
  99.  
  100.  
  101. /* both min() and MIN() are defined by some compilers - grr ! */
  102. #define GPMIN(a,b)    ((a)>(b) ? (b) : (a))
  103.  
  104. #define Inc_c_token if (++c_token >= num_tokens)    \
  105.                         int_error ("Syntax error", c_token);
  106.  
  107.  
  108. /*
  109.  * IMHO, code is getting too cluttered with repeated chunks of
  110.  * code. Some macros to simplify, I hope.
  111.  *
  112.  * do { } while(0) is comp.lang.c recommendation for complex macros
  113.  * also means that break can be specified as an action, and it will
  114.  * 
  115.  */
  116.  
  117.  
  118.  
  119. /* store VALUE or log(VALUE) in STORE, set TYPE as appropriate
  120.  * Do OUT_ACTION or UNDEF_ACTION as appropriate
  121.  * adjust range provided type is INRANGE (ie dont adjust y if x is outrange
  122.  * VALUE must not be same as STORE
  123.  */
  124.  
  125. #define STORE_AND_FIXUP_RANGE(STORE, VALUE, TYPE, MIN, MAX, AUTO, OUT_ACTION, UNDEF_ACTION)\
  126. do { STORE=VALUE; \
  127.      if (TYPE != INRANGE) break;  /* dont set y range if x is outrange, for example */ \
  128.      if ( VALUE<MIN ) \
  129.       if (AUTO & 1) MIN=VALUE; else { TYPE=OUTRANGE; OUT_ACTION; break; }  \
  130.      if ( VALUE>MAX ) \
  131.      if (AUTO & 2) MAX=VALUE; else { TYPE=OUTRANGE; OUT_ACTION; }   \
  132. } while(0)
  133.  
  134. #define UPDATE_RANGE(TEST,OLD,NEW,AXIS) \
  135. do { if (TEST) \
  136.     if (log_array[AXIS]) OLD = pow(base_array[AXIS], NEW); \
  137.    else OLD = NEW; }while(0)
  138.      
  139. /* use this instead empty macro arguments to work around NeXT cpp bug */
  140. /* if this fails on any system, we might use ((void)0) */
  141.  
  142. #define NOOP /* */
  143.  
  144.  
  145. #define inrange(z,min,max) ((min<max) ? ((z>=min)&&(z<=max)) : ((z>=max)&&(z<=min)) )
  146.  
  147. #define spline_coeff_size 4
  148. typedef double spline_coeff[spline_coeff_size];
  149. typedef double five_diag[5];
  150.  
  151. static double *cp_binomial __P((struct curve_points *cp));
  152. static double s_pow __P((double base, unsigned int exponent));
  153. static void eval_bezier __P((struct curve_points *cp, double sr, coordval *px, coordval *py, double *c));
  154. static void do_bezier __P((struct curve_points *cp, double *bc));
  155. static spline_coeff *cp_tridiag __P((struct curve_points *plot));
  156. static int solve_five_diag __P((five_diag m[], double r[], double x[], int n));
  157. static spline_coeff *cp_approx_spline __P((struct curve_points *plot));
  158. static void do_cubic __P((struct curve_points *plot, spline_coeff *sc));
  159. static int compare_points __P((struct coordinate *p1, struct coordinate *p2));
  160.  
  161. /*
  162.  * build up a cntr_struct list from curve_points
  163.  * this funtion is only used for the alternate entry point to
  164.  * Gershon's code and thus commented out
  165.  
  166. static struct cntr_struct *conv2cntr(struct curve_points *plot)
  167. {
  168.   int i;
  169.   struct cntr_struct *c,*oc;
  170.  
  171.   if(!plot) return(NULL);
  172.   oc=NULL;
  173.   for (i=plot->p_count-1;i >= 0; i--){
  174.      c=(struct cntr_struct *)alloc(sizeof(struct cntr_struct),"Spline list");
  175.      c->next=oc;
  176.      c->X=plot->points[i].x;
  177.      c->Y=plot->points[i].y;
  178.      oc=c;
  179.   }
  180.  
  181.   return(c);
  182. }
  183. */
  184.  
  185. /*
  186.  * cp_binomial() computes the binomial coefficients needed for BEZIER stuff
  187.  *   and stores them into an array which is hooked to sdat.
  188.  * (MGR 1992)
  189.  */
  190. static double *cp_binomial(cp)
  191. struct curve_points *cp;
  192. {
  193.    register double *coeff;
  194.    register int n,k;
  195.    int e;
  196.  
  197.    e=cp->p_count;        /* well we're going from k=0 to k=p_count-1 */
  198.    coeff=(double *)alloc(e*sizeof(double),"bezier coefficients");
  199.  
  200.    n=cp->p_count-1;
  201.    e=n/2;
  202.    coeff[0]=1.0;
  203.  
  204.    for(k=0;k<e;k++)
  205.       coeff[k+1] = coeff[k] * ((double)(n-k))/((double)(k+1));
  206.  
  207.    for(k=n;k>=e;k--) coeff[k] = coeff[n-k];
  208.  
  209.    return(coeff);
  210. }
  211.  
  212. /* This is a subfunction of do_bezier() for BEZIER style computations.
  213.  * It is passed the stepration (STEP/MAXSTEPS) and the addresses of
  214.  * the double values holding the next x and y coordinates.
  215.  * (MGR 1992)
  216.  */
  217.  
  218. /*
  219.  * well, this routine runs faster with the 68040 striptease FPU
  220.  * and it handles zeroes correctly - there had been some trouble with TC
  221.  */
  222.  
  223. #ifdef AMIGA_SC_6_1
  224. __inline
  225. #endif
  226. static double s_pow(base, exponent)
  227. double base;
  228. unsigned int exponent;
  229. {
  230.   int i;
  231.   double x,y;
  232.  
  233.   static int spoken=0;
  234.  
  235.   if(exponent==0) return(1.0);
  236.   if(base==0.0) return(0.0);
  237.   x=base;
  238.   for(i=1;i<exponent;i++) x*=base;
  239.  
  240.   /* consider i in binary = abcd
  241.    * x^i = x^(8a+4b+2c+d) = x^8a * x^4b * x^2b * x^d
  242.    *                      = a?x^2^2^2:1 * b?x^2^2:1 + ...
  243.    * so for each bit in exponent, square x, multiplying it into accumulator
  244.      *
  245.     * - probably not clear - that's why I'm comparing the two versions !
  246.     */
  247.  
  248.   if (!spoken) {
  249.       fprintf(stderr, "Please note that s_pow in interpol.c is comparing two algorithms.\nThis message will be taken out (I hope) before release\n");
  250.       spoken=1;
  251.   }
  252.  
  253.     y=1;
  254.     while (exponent)
  255.     {
  256.           if (exponent&1)
  257.               y *= base;
  258.           base *= base;
  259.           exponent >>= 1;  /* if exponent was signed, this could be trouble ! */
  260.     }
  261.  
  262.     if (fabs(x-y) > 0.001)
  263.         fprintf(stderr, "uh oh - div got s_pow wrong in interpol.c : %f != %f\n", y, x);
  264.       
  265.   return(x);
  266. }
  267.  
  268. static void eval_bezier(cp,sr,px,py,c)
  269. struct curve_points *cp;
  270. double sr;
  271. coordval *px;
  272. coordval *py;
  273. double *c;
  274. {
  275.    int i;
  276.    int n=cp->p_count-1;
  277.    double lx=0.0, ly=0.0;
  278.     double sr_to_the_i = 1.0;
  279.     double dsr_to_the_di = s_pow(1-sr, n);
  280.     double reciprocal_dsr = 1.0/(1-sr);
  281.  
  282.    for(i=0;i<=n;i++){
  283.       double u = c[i] * sr_to_the_i * s_pow(1-sr,n-i);
  284.       lx+=cp->points[i].x*u;
  285.       ly+=cp->points[i].y*u;
  286.       sr_to_the_i *= sr;
  287.       dsr_to_the_di *= reciprocal_dsr;
  288.    }
  289.  
  290.    *px = lx;
  291.    *py = ly;
  292. }
  293.  
  294. /*
  295.  * generate a new set of coordinates representing the bezier curve and
  296.  * set it to the plot
  297.  */
  298.  
  299. static void do_bezier(cp,bc)
  300. struct curve_points *cp;
  301. double *bc;
  302. {
  303.   struct coordinate *s_tbl;
  304.   int i;
  305.   double x,y;
  306.   int xaxis = cp->x_axis;
  307.   int yaxis = cp->y_axis;
  308.  
  309.   /* min and max in internal (eg logged) co-ordinates. We update
  310.    * these, then update the external extrema in user co-ordinates
  311.    * at the end.
  312.    */
  313.    
  314.   double ixmin,ixmax,iymin,iymax;
  315.   double sxmin,sxmax,symin,symax; /* starting values of above */
  316.  
  317.   if (log_array[xaxis])
  318.   {
  319.           ixmin = sxmin = log(min_array[xaxis]) / log_base_array[xaxis];
  320.           ixmax = sxmax = log(max_array[xaxis]) / log_base_array[xaxis];
  321.   }
  322.   else
  323.   {
  324.           ixmin = sxmin = min_array[xaxis];
  325.           ixmax = sxmax = max_array[xaxis];
  326.   }
  327.  
  328.   if (log_array[yaxis])
  329.   {
  330.           iymin = symin = log(min_array[yaxis]) / log_base_array[yaxis];
  331.           iymax = symax = log(max_array[yaxis]) / log_base_array[yaxis];
  332.   }
  333.   else
  334.   {
  335.           iymin = symin = min_array[yaxis];
  336.           iymax = symax = max_array[yaxis];
  337.   }
  338.  
  339.   s_tbl=(struct coordinate *)alloc(samples*sizeof(struct coordinate),"bezier table");
  340.  
  341.   for(i=0;i<samples;i++){
  342.      eval_bezier(cp,(double)i/(double)(samples-1),&x,&y,bc);
  343.  
  344.      /* now we have to store the points and adjust the ranges */
  345.  
  346.      s_tbl[i].type=INRANGE;
  347.      STORE_AND_FIXUP_RANGE(s_tbl[i].x, x, s_tbl[i].type, ixmin, ixmax, auto_array[xaxis], NOOP, continue );
  348.      STORE_AND_FIXUP_RANGE(s_tbl[i].y, y, s_tbl[i].type, iymin, iymax, auto_array[yaxis], NOOP, NOOP );
  349.  
  350.      s_tbl[i].xlow = s_tbl[i].xhigh = s_tbl[i].x;
  351.      s_tbl[i].ylow = s_tbl[i].yhigh = s_tbl[i].y;
  352.      
  353.      s_tbl[i].z = -1;
  354.   }
  355.   free((char *)cp->points);
  356.   cp->points=s_tbl;
  357.   cp->p_count=samples;
  358.   cp->p_max=samples;
  359.  
  360.  
  361.     UPDATE_RANGE(ixmax > sxmax, max_array[xaxis], ixmax, xaxis);
  362.     UPDATE_RANGE(ixmin < sxmin, min_array[xaxis], ixmin, xaxis);
  363.     UPDATE_RANGE(iymax > symax, max_array[yaxis], iymax, yaxis);
  364.     UPDATE_RANGE(iymin < symin, min_array[yaxis], iymin, yaxis);
  365. }
  366.  
  367. /*
  368.  * call contouring routines -- main entry
  369.  */
  370.  
  371. /*
  372.  * it should be like this, but it doesn't run. If you find out why,
  373.  * contact me: mgr@asgard.bo.open.de or Lars Hanke 2:243/4802.22@fidonet
  374.  *
  375.  * Well, all this had originally been inside contour.c, so maybe links
  376.  * to functions and of contour.c are broken.
  377.  *
  378.  
  379. --> this denotes comment nests
  380.  
  381. void gen_interp(struct curve_points *plot)
  382. {
  383.   struct cntr_struct *my_cntr = NULL;
  384.   struct gnuplot_contours *save_list = contour_list;
  385.  
  386.   int i;
  387.   double x;
  388.   double *bc;
  389.   int num_pts = num_approx_pts;
  390.  
  391.   hermit_table=NULL;
  392.   contour_list=NULL;
  393.   my_cntr=conv2cntr(plot);
  394.   if(!my_cntr) return;
  395.  
  396.   switch(plot->plot_smooth){
  397.      case CSPLINES: calc_hermit_table();
  398.                     interp_kind=INTERP_CUBIC;
  399.                     break;
  400.      case BEZIER:
  401.      case SBEZIER:  bc=cp_binomial(plot);
  402.                     do_bezier(plot,bc);
  403.                     free((char *)bc);
  404. -->                  no break
  405.      default:       free_contour(my_cntr);
  406.                     return;
  407.   }
  408.  
  409.   num_approx_pts=samples/(plot->p_count);
  410.  
  411. --> could not find any place where x_min to y_max were used, so I set
  412. --> them to zero
  413.  
  414.   put_contour(my_cntr,0,0,0,0,0,OPEN_CONTOUR);
  415.  
  416.   free(plot->points);
  417.  
  418.   plot->p_max=contour_list->num_pts;
  419.   plot->p_count=contour_list->num_pts;
  420.   plot->points=contour_list->coords;
  421.   free(contour_list);
  422.   contour_list=save_list;
  423.   num_approx_pts=num_pts;
  424.  
  425.   if(hermit_table){
  426.      free(hermit_table);
  427.      hermit_table=NULL;
  428.   }
  429.   free_contour(my_cntr);
  430.  
  431.   for(i=0;i<plot->p_count;i++){
  432.      x=(is_log_x)? pow(base_log_x,plot->points[i].x) : plot->points[i].x;
  433.      if((x >= xmin) && (x <= xmax)) plot->points[i].type=INRANGE;
  434.      else plot->points[i].type=OUTRANGE;
  435.   }
  436.  
  437. }
  438.  *
  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.     five_diag *hv;
  465.     int i;
  466.  
  467.     hv=(five_diag *)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.     x[n-1] = hv[n][4];
  503.     x[n-2] = hv[n-1][4] - hv[n-2][1]*x[n-1];
  504.     for( i=n-3; i>=0; i--)
  505.       x[i] = hv[i+1][4] - hv[i][1]*x[i+1] - hv[i][2]*x[i+2];
  506.  
  507.     free(hv);
  508.     return TRUE;
  509. }
  510.  
  511. /*
  512.  * Calculation of approximation cubic splines
  513.  * Input:  x[i], y[i], weights z[i]
  514.  *         
  515.  * Returns matrix of spline coefficients
  516.  */
  517. static spline_coeff *cp_approx_spline(plot)
  518. struct curve_points *plot;
  519. {
  520.   spline_coeff  *sc;
  521.   five_diag     *m;
  522.   double        *r, *x, *h;
  523.   
  524.   int i;
  525.   
  526.   sc=(spline_coeff *)alloc((plot->p_count)*sizeof(spline_coeff),"spline matrix");
  527.  
  528.   if (plot->p_count < 3)
  529.     int_error("Can't calculate approximation splines, need more than 3 points", NO_CARET);
  530.  
  531.   for( i=0; i<=plot->p_count-1; i++)
  532.     if (plot->points[i].z <=0)
  533.       int_error("Can't calculate approximation splines, all weights have to be > 0", NO_CARET);
  534.   
  535.   m = (five_diag *)alloc((plot->p_count-2)*sizeof(five_diag),"spline help matrix");
  536.  
  537.   r = (double *)alloc((plot->p_count-2)*sizeof(double),"spline right side");
  538.   x = (double *)alloc((plot->p_count-2)*sizeof(double),"spline solution vector");
  539.   h = (double *)alloc((plot->p_count-1)*sizeof(double),"spline help vector");
  540.   
  541.   for( i=0; i<=plot->p_count-2; i++)
  542.     h[i] = plot->points[i+1].x - plot->points[i].x;
  543.   
  544.   /* set up the matrix and the vector*/
  545.   
  546.   for( i=0; i<=plot->p_count-3; i++){
  547.     r[i] = 3*((plot->points[i+2].y-plot->points[i+1].y)/h[i+1]
  548.           -(plot->points[i+1].y-plot->points[i].y)/h[i]); 
  549.     
  550.     if ( i < 2 )
  551.       m[i][0] = 0;
  552.     else
  553.       m[i][0] = 6 / plot->points[i].z / h[i-1] / h[i];
  554.     
  555.     if ( i < 1 )
  556.       m[i][1] = 0;
  557.     else
  558.       m[i][1] = h[i] - 6 / plot->points[i].z / h[i] * ( 1/h[i-1] + 1/h[i])
  559.     - 6 / plot->points[i+1].z / h[i] * ( 1/h[i] + 1/h[i+1]);
  560.     
  561.     m[i][2] = 2 * ( h[i] + h[i+1] )
  562.       + 6 / plot->points[i].z / h[i] / h[i]
  563.     + 6 / plot->points[i+1].z * ( 1/h[i] + 1/h[i+1] ) * ( 1/h[i] + 1/h[i+1] )
  564.       + 6 / plot->points[i+2].z / h[i+1] / h[i+1];
  565.     
  566.     if ( i > plot->p_count-4 )
  567.       m[i][3] = 0;
  568.     else
  569.       m[i][3] = h[i+1] - 6 / plot->points[i+1].z / h[i+1] * ( 1/h[i] + 1/h[i+1] )
  570.     - 6 / plot->points[i+2].z / h[i+1] * ( 1/h[i+1] + 1/h[i+2] );
  571.     
  572.     if ( i > plot->p_count-5 )
  573.       m[i][4] = 0;
  574.     else
  575.       m[i][4] = 6 / plot->points[i+2].z / h[i+1] / h[i+2];
  576.   }
  577.   
  578.   /* solve the matrix */
  579.   if (!solve_five_diag(m, r, x, plot->p_count-2)){
  580.     free(h);
  581.     free(x);
  582.     free(r);
  583.     free(m);
  584.     int_error("Can't calculate approximation splines", NO_CARET);
  585.   }
  586.   
  587.   sc[0][2] = 0;
  588.   for( i=1; i<=plot->p_count-2; i++ )
  589.     sc[i][2] = x[i-1];
  590.   sc[plot->p_count-1][2] = 0;
  591.   
  592.   sc[0][0] = plot->points[0].y + 2 / plot->points[0].z / h[0] * ( sc[0][2] - sc[1][2] );
  593.   for( i=1; i<=plot->p_count-2; i++ )
  594.     sc[i][0] = plot->points[i].y - 2 / plot->points[i].z *
  595.       ( sc[i-1][2] / h[i-1] 
  596.        - sc[i][2] * ( 1/h[i-1] + 1/h[i] )
  597.        + sc[i+1][2] / h[i] );
  598.   sc[plot->p_count-1][0] = plot->points[plot->p_count-1].y 
  599.     - 2 / plot->points[plot->p_count-1].z / h[plot->p_count-2] 
  600.       * ( sc[plot->p_count-2][2] - sc[plot->p_count-1][2] );
  601.   
  602.   for( i=0; i<=plot->p_count-2; i++ ){
  603.     sc[i][1] = ( sc[i+1][0] - sc[i][0] ) / h[i]
  604.       - h[i] / 3 * ( sc[i+1][2] + 2 * sc[i][2] );
  605.     sc[i][3] = ( sc[i+1][2] - sc[i][2] ) / 3 / h[i];
  606.   }
  607.   
  608.   free(h);
  609.   free(x);
  610.   free(r);
  611.   free(m);
  612.   
  613.   return(sc);
  614. }
  615.  
  616. /*
  617.  * Calculation of cubic splines
  618.  * Input:  x[i], y[i]
  619.  * This can be treated as a special case of approximation cubic splines, with
  620.  * all weights -> infinity.
  621.  *         
  622.  * Returns matrix of spline coefficients
  623.  */
  624. static spline_coeff *cp_tridiag(plot)
  625. struct curve_points *plot;
  626. {
  627.   spline_coeff *sc;
  628.   tri_diag     *m;
  629.   double       *r, *x, *h;
  630.   
  631.   int i;
  632.   
  633.   sc=(spline_coeff *)alloc((plot->p_count)*sizeof(spline_coeff),"spline matrix");
  634.   
  635.   m = (tri_diag *)alloc((plot->p_count-2)*sizeof(tri_diag),"spline help matrix");
  636.   
  637.   r = (double *)alloc((plot->p_count-2)*sizeof(double),"spline right side");
  638.   x = (double *)alloc((plot->p_count-2)*sizeof(double),"spline solution vector");
  639.   h = (double *)alloc((plot->p_count-1)*sizeof(double),"spline help vector");
  640.   
  641.   for( i=0; i<=plot->p_count-2; i++)
  642.     h[i] = plot->points[i+1].x - plot->points[i].x;
  643.   
  644.   /* set up the matrix and the vector*/
  645.   
  646.   for( i=0; i<=plot->p_count-3; i++){
  647.     r[i] = 3*((plot->points[i+2].y-plot->points[i+1].y)/h[i+1]
  648.           -(plot->points[i+1].y-plot->points[i].y)/h[i]); 
  649.     
  650.     if ( i < 1 )
  651.       m[i][0] = 0;
  652.     else
  653.       m[i][0] = h[i];
  654.     
  655.     m[i][1] = 2 * ( h[i] + h[i+1] );
  656.     
  657.     if ( i > plot->p_count-4 )
  658.       m[i][2] = 0;
  659.     else
  660.       m[i][2] = h[i+1];
  661.   }
  662.   
  663.   /* solve the matrix */
  664.   if (!solve_tri_diag(m, r, x, plot->p_count-2)){
  665.     free(h);
  666.     free(x);
  667.     free(r);
  668.     free(m);
  669.     int_error("Can't calculate cubic splines", NO_CARET);
  670.   }
  671.   
  672.   sc[0][2] = 0;
  673.   for( i=1; i<=plot->p_count-2; i++ )
  674.     sc[i][2] = x[i-1];
  675.   sc[plot->p_count-1][2] = 0;
  676.   
  677.   for( i=0; i<=plot->p_count-1; i++ )
  678.     sc[i][0] = plot->points[i].y;
  679.   
  680.   for( i=0; i<=plot->p_count-2; i++ ){
  681.     sc[i][1] = ( sc[i+1][0] - sc[i][0] ) / h[i]
  682.       - h[i] / 3 * ( sc[i+1][2] + 2 * sc[i][2] );
  683.     sc[i][3] = ( sc[i+1][2] - sc[i][2] ) / 3 / h[i];
  684.   }
  685.   
  686.   free(h);
  687.   free(x);
  688.   free(r);
  689.   free(m);
  690.  
  691.   return(sc);
  692. }
  693.  
  694. static void do_cubic(plot,sc)
  695. struct curve_points *plot;
  696. spline_coeff *sc;
  697. {
  698.   struct coordinate *s_tbl;
  699.   double xdiff,temp,x,y;
  700.   int i,l;
  701.   int xaxis = plot->x_axis;
  702.   int yaxis = plot->y_axis;
  703.  
  704.   /* min and max in internal (eg logged) co-ordinates. We update
  705.    * these, then update the external extrema in user co-ordinates
  706.    * at the end.
  707.    */
  708.    
  709.   double ixmin,ixmax,iymin,iymax;
  710.   double sxmin,sxmax,symin,symax; /* starting values of above */
  711.  
  712.   if (log_array[xaxis])
  713.   {
  714.           ixmin = sxmin = log(min_array[xaxis]) / log_base_array[xaxis];
  715.           ixmax = sxmax = log(max_array[xaxis]) / log_base_array[xaxis];
  716.   }
  717.   else
  718.   {
  719.           ixmin = sxmin = min_array[xaxis];
  720.           ixmax = sxmax = max_array[xaxis];
  721.   }
  722.  
  723.   if (log_array[yaxis])
  724.   {
  725.           iymin = symin = log(min_array[yaxis]) / log_base_array[yaxis];
  726.           iymax = symax = log(max_array[yaxis]) / log_base_array[yaxis];
  727.   }
  728.   else
  729.   {
  730.           iymin = symin = min_array[yaxis];
  731.           iymax = symax = max_array[yaxis];
  732.   }
  733.  
  734.   s_tbl=(struct coordinate *)alloc(samples*sizeof(struct coordinate),"spline table");
  735.  
  736.   l = 0;
  737.  
  738.   xdiff = (plot->points[plot->p_count-1].x -plot->points[0].x)/(samples-1);
  739.  
  740.   for(i = 0; i < samples; i++) {
  741.      x = plot->points[0].x + i*xdiff;
  742.  
  743.      while((x >= plot->points[l+1].x) && (l < plot->p_count-2)) l++;
  744.      
  745.      temp = x - plot->points[l].x;
  746.      
  747.      y = ((sc[l][3]*temp+sc[l][2])*temp+sc[l][1])*temp+sc[l][0];
  748.  
  749.      s_tbl[i].type=INRANGE;
  750.      STORE_AND_FIXUP_RANGE(s_tbl[i].x, x, s_tbl[i].type, ixmin, ixmax, auto_array[xaxis], NOOP, continue );
  751.      STORE_AND_FIXUP_RANGE(s_tbl[i].y, y, s_tbl[i].type, iymin, iymax, auto_array[yaxis], NOOP, NOOP );
  752.  
  753.      s_tbl[i].xlow = s_tbl[i].xhigh = s_tbl[i].x;
  754.      s_tbl[i].ylow = s_tbl[i].yhigh = s_tbl[i].y;
  755.  
  756.      s_tbl[i].z = -1;
  757.  
  758.   }
  759.  
  760.   free(plot->points);
  761.   plot->points = s_tbl;
  762.   plot->p_count = samples;
  763.   plot->p_max = samples;
  764.  
  765.     UPDATE_RANGE(ixmax > sxmax, max_array[xaxis], ixmax, xaxis);
  766.     UPDATE_RANGE(ixmin < sxmin, min_array[xaxis], ixmin, xaxis);
  767.     UPDATE_RANGE(iymax > symax, max_array[yaxis], iymax, yaxis);
  768.     UPDATE_RANGE(iymin < symin, min_array[yaxis], iymin, yaxis);
  769. }
  770.  
  771. /*
  772.  * This is the main entry point used. As stated in the header, it is fine,
  773.  * but I'm not too happy with it.
  774.  */
  775.  
  776. void gen_interp(plot)
  777. struct curve_points *plot;
  778. {
  779.  
  780.   spline_coeff *sc;
  781.   double *bc;
  782.  
  783.   switch(plot->plot_smooth){
  784.      case CSPLINES:
  785.                     sc=cp_tridiag(plot);
  786.                     do_cubic(plot,sc);
  787.                     free(sc);
  788.                     break;
  789.      case ACSPLINES:
  790.                     sc=cp_approx_spline(plot);
  791.                     do_cubic(plot,sc);
  792.                     free(sc);
  793.                     break;
  794.  
  795.      case BEZIER:
  796.      case SBEZIER:  bc=cp_binomial(plot);
  797.                     do_bezier(plot,bc);
  798.                     free((char *)bc);
  799.             break;
  800.      default: /* keep gcc -Wall quiet */
  801.             ;
  802.   }
  803.  
  804.   return;
  805.  
  806. }
  807.  
  808. /* sort_points
  809.  * sort data succession for further evaluation by plot_splines, etc.
  810.  * This routine is mainly introduced for compilers *NOT* supporting the
  811.  * UNIX qsort() routine. You can then easily replace it by more convenient
  812.  * stuff for your compiler.
  813.  * (MGR 1992)
  814.  */
  815.  
  816. static int compare_points(p1,p2)
  817. struct coordinate *p1;
  818. struct coordinate *p2;
  819. {
  820.    if(p1->x > p2->x) return(1);
  821.    if(p1->x < p2->x) return(-1);
  822.    return(0);
  823. }
  824.  
  825. void sort_points(plot)
  826. struct curve_points *plot;
  827. {
  828.    qsort((char *)(plot->points), plot->p_count, sizeof(struct coordinate),
  829.          (sortfunc)compare_points);
  830. }
  831.  
  832. /*
  833.  * cp_implode() if averaging is selected this function computes the new
  834.  *              entries and shortens the whole thing to the necessary
  835.  *              size
  836.  * MGR Addendum
  837.  */
  838.  
  839. void cp_implode(cp)
  840. struct curve_points *cp;
  841. {
  842.    int i,j,k;
  843.    double x,y,sux,slx,suy,sly;
  844.    enum coord_type dot;
  845.  
  846.    k=0;
  847.    j=0;
  848.    for(i=0;i<cp->p_count;i++){
  849.       if(!k){
  850.          x=cp->points[i].x;
  851.          y=cp->points[i].y;
  852.          sux=cp->points[i].xhigh;
  853.          slx=cp->points[i].xlow;
  854.          suy=cp->points[i].yhigh;
  855.          sly=cp->points[i].ylow;
  856.          dot=INRANGE;
  857.          if(cp->points[i].type!=INRANGE) dot=UNDEFINED;
  858.          k=1;
  859.       } else if(cp->points[i].x==x){
  860.          y+=cp->points[i].y;
  861.          sux+=cp->points[i].xhigh;
  862.          slx+=cp->points[i].xlow;
  863.          suy+=cp->points[i].yhigh;
  864.          sly+=cp->points[i].ylow;
  865.          if(cp->points[i].type!=INRANGE) dot=UNDEFINED;
  866.          k++;
  867.       } else {
  868.          cp->points[j].x=x;
  869.          cp->points[j].y=y/(double)k;
  870.          cp->points[j].xhigh=sux/(double)k;
  871.          cp->points[j].xlow=slx/(double)k;
  872.          cp->points[j].yhigh=suy/(double)k;
  873.          cp->points[j].ylow=sly/(double)k;
  874.          cp->points[j].type=INRANGE;
  875.          if(dot!=INRANGE){
  876.             if((cp->points[j].x>xmax) || (cp->points[j].x<xmin))
  877.                cp->points[j].type=OUTRANGE;
  878.             else if(!autoscale_ly){
  879.                if((cp->points[j].y>ymax) || (cp->points[j].y<ymin))
  880.                   cp->points[j].type=OUTRANGE;
  881.             } else cp->points[j].type=INRANGE;
  882.          }
  883.          j++;  /* next valid entry */
  884.          k=0;  /* to read */
  885.          i--;  /* from this (-> last after for(;;)) entry */
  886.       }
  887.    }
  888.    if(k){
  889.       cp->points[j].x=x;
  890.       cp->points[j].y=y/(double)k;
  891.       cp->points[j].xhigh=sux/(double)k;
  892.       cp->points[j].xlow=slx/(double)k;
  893.       cp->points[j].yhigh=suy/(double)k;
  894.       cp->points[j].ylow=sly/(double)k;
  895.       cp->points[j].type=INRANGE;
  896.       if(dot!=INRANGE){
  897.          if((cp->points[j].x>xmax) || (cp->points[j].x<xmin))
  898.             cp->points[j].type=OUTRANGE;
  899.          else if(!autoscale_ly){
  900.             if((cp->points[j].y>ymax) || (cp->points[j].y<ymin))
  901.                cp->points[j].type=OUTRANGE;
  902.          } else cp->points[j].type=INRANGE;
  903.       }
  904.       j++;  /* next valid entry */
  905.    }
  906.    cp->p_count=j;
  907.    cp_extend(cp,j);
  908. }
  909.