"Fossies" - the Fresh Open Source Software Archive

Member "grace-5.1.25/src/computils.c" (14 Feb 2015, 43273 Bytes) of package /linux/misc/grace-5.1.25.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "computils.c" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 5.1.24_vs_5.1.25.

    1 /*
    2  * Grace - GRaphing, Advanced Computation and Exploration of data
    3  * 
    4  * Home page: http://plasma-gate.weizmann.ac.il/Grace/
    5  * 
    6  * Copyright (c) 1991-1995 Paul J Turner, Portland, OR
    7  * Copyright (c) 1996-2000 Grace Development Team
    8  * 
    9  * Maintained by Evgeny Stambulchik
   10  * 
   11  * 
   12  *                           All Rights Reserved
   13  * 
   14  *    This program is free software; you can redistribute it and/or modify
   15  *    it under the terms of the GNU General Public License as published by
   16  *    the Free Software Foundation; either version 2 of the License, or
   17  *    (at your option) any later version.
   18  * 
   19  *    This program is distributed in the hope that it will be useful,
   20  *    but WITHOUT ANY WARRANTY; without even the implied warranty of
   21  *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   22  *    GNU General Public License for more details.
   23  * 
   24  *    You should have received a copy of the GNU General Public License
   25  *    along with this program; if not, write to the Free Software
   26  *    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
   27  */
   28 
   29 /*
   30  *
   31  * procedures for performing transformations from the command
   32  * line interpreter and the GUI.
   33  *
   34  */
   35 
   36 #include <config.h>
   37 #include <cmath.h>
   38 
   39 #include <stdio.h>
   40 #include <stdlib.h>
   41 #include <string.h>
   42 
   43 #include "globals.h"
   44 #include "utils.h"
   45 #include "draw.h"
   46 #include "ssdata.h"
   47 #include "graphs.h"
   48 #include "parser.h"
   49 #include "protos.h"
   50 
   51 static void forwarddiff(double *x, double *y, double *resx, double *resy, int n);
   52 static void backwarddiff(double *x, double *y, double *resx, double *resy, int n);
   53 static void centereddiff(double *x, double *y, double *resx, double *resy, int n);
   54 int get_points_inregion(int rno, int invr, int len, double *x, double *y, int *cnt, double **xt, double **yt);
   55 
   56 static char buf[256];
   57 
   58 
   59 void do_fourier_command(int gno, int setno, int ftype, int ltype)
   60 {
   61     switch (ftype) {
   62     case FFT_DFT:
   63     do_fourier(gno, setno, 0, 0, ltype, 0, 0, 0);
   64     break;
   65     case FFT_INVDFT       :
   66     do_fourier(gno, setno, 0, 0, ltype, 1, 0, 0);
   67     break;
   68     case FFT_FFT:
   69     do_fourier(gno, setno, 1, 0, ltype, 0, 0, 0);
   70     break;
   71     case FFT_INVFFT       :
   72     do_fourier(gno, setno, 1, 0, ltype, 1, 0, 0);
   73     break;
   74     }
   75 }
   76 
   77 
   78 /*
   79  * evaluate a formula
   80  */
   81 int do_compute(int gno, int setno, int graphto, int loadto, char *rarray, char *fstr)
   82 {
   83     if (is_set_active(gno, setno)) {
   84     if (gno != graphto || setno != loadto) {
   85         if (copysetdata(gno, setno, graphto, loadto) != RETURN_SUCCESS) {
   86             return RETURN_FAILURE;
   87             }
   88         }
   89     filter_set(graphto, loadto, rarray);
   90         set_parser_setno(graphto, loadto);
   91         if (scanner(fstr) != RETURN_SUCCESS) {
   92         if (graphto != gno || loadto != setno) {
   93         killset(graphto, loadto);
   94         }
   95         return RETURN_FAILURE;
   96     } else {
   97         set_dirtystate();
   98             return RETURN_SUCCESS;
   99         }
  100     } else {
  101         return RETURN_FAILURE;
  102     }
  103 }
  104 
  105 /*
  106  * forward, backward and centered differences
  107  */
  108 static void forwarddiff(double *x, double *y, double *resx, double *resy, int n)
  109 {
  110     int i, eflag = 0;
  111     double h;
  112 
  113     for (i = 1; i < n; i++) {
  114     resx[i - 1] = x[i - 1];
  115     h = x[i - 1] - x[i];
  116     if (h == 0.0) {
  117         resy[i - 1] = - MAXNUM;
  118         eflag = 1;
  119     } else {
  120         resy[i - 1] = (y[i - 1] - y[i]) / h;
  121     }
  122     }
  123     if (eflag) {
  124     errmsg("Warning: infinite slope, check set status before proceeding");
  125     }
  126 }
  127 
  128 static void backwarddiff(double *x, double *y, double *resx, double *resy, int n)
  129 {
  130     int i, eflag = 0;
  131     double h;
  132 
  133     for (i = 0; i < n - 1; i++) {
  134     resx[i] = x[i];
  135     h = x[i + 1] - x[i];
  136     if (h == 0.0) {
  137         resy[i] = - MAXNUM;
  138         eflag = 1;
  139     } else {
  140         resy[i] = (y[i + 1] - y[i]) / h;
  141     }
  142     }
  143     if (eflag) {
  144     errmsg("Warning: infinite slope, check set status before proceeding");
  145     }
  146 }
  147 
  148 static void centereddiff(double *x, double *y, double *resx, double *resy, int n)
  149 {
  150     int i, eflag = 0;
  151     double h1, h2;
  152 
  153     for (i = 1; i < n - 1; i++) {
  154     resx[i - 1] = x[i];
  155     h1 = x[i] - x[i - 1];
  156     h2 = x[i + 1] - x[i];
  157     if (h1 + h2 == 0.0) {
  158         resy[i - 1] = - MAXNUM;
  159         eflag = 1;
  160     } else {
  161         resy[i - 1] = (y[i + 1] - y[i - 1]) / (h1 + h2);
  162     }
  163     }
  164     if (eflag) {
  165     errmsg("Warning: infinite slope, check set status before proceeding");
  166     }
  167 }
  168 
  169 static void seasonaldiff(double *x, double *y,
  170              double *resx, double *resy, int n, int period)
  171 {
  172     int i;
  173 
  174     for (i = 0; i < n - period; i++) {
  175     resx[i] = x[i];
  176     resy[i] = y[i] - y[i + period];
  177     }
  178 }
  179 
  180 /*
  181  * trapezoidal rule
  182  */
  183 double trapint(double *x, double *y, double *resx, double *resy, int n)
  184 {
  185     int i;
  186     double sum = 0.0;
  187     double h;
  188 
  189     if (n < 2) {
  190         return 0.0;
  191     }
  192     
  193     if (resx != NULL) {
  194         resx[0] = x[0];
  195     }
  196     if (resy != NULL) {
  197         resy[0] = 0.0;
  198     }
  199     for (i = 1; i < n; i++) {
  200     h = (x[i] - x[i - 1]);
  201     if (resx != NULL) {
  202         resx[i] = x[i];
  203     }
  204     sum = sum + h * (y[i - 1] + y[i]) * 0.5;
  205     if (resy != NULL) {
  206         resy[i] = sum;
  207     }
  208     }
  209     return sum;
  210 }
  211 
  212 /*
  213  * apply a digital filter
  214  */
  215 void do_digfilter(int set1, int set2)
  216 {
  217     int digfiltset;
  218 
  219     if (!(is_set_active(get_cg(), set1) && is_set_active(get_cg(), set2))) {
  220     errmsg("Set not active");
  221     return;
  222     }
  223     if ((getsetlength(get_cg(), set1) < 3) || (getsetlength(get_cg(), set2) < 3)) {
  224     errmsg("Set length < 3");
  225     return;
  226     }
  227     digfiltset = nextset(get_cg());
  228     if (digfiltset != (-1)) {
  229     activateset(get_cg(), digfiltset);
  230     setlength(get_cg(), digfiltset, getsetlength(get_cg(), set1) - getsetlength(get_cg(), set2) + 1);
  231     sprintf(buf, "Digital filter from set %d applied to set %d", set2, set1);
  232     filterser(getsetlength(get_cg(), set1),
  233           getx(get_cg(), set1),
  234           gety(get_cg(), set1),
  235           getx(get_cg(), digfiltset),
  236           gety(get_cg(), digfiltset),
  237           gety(get_cg(), set2),
  238           getsetlength(get_cg(), set2));
  239     setcomment(get_cg(), digfiltset, buf);
  240     }
  241 }
  242 
  243 /*
  244  * linear convolution
  245  */
  246 void do_linearc(int gno1, int set1, int gno2, int set2)
  247 {
  248     int linearcset, i, itmp, cg = get_cg();
  249     double *xtmp;
  250 
  251     if (!(is_set_active(gno1, set1) && is_set_active(gno2, set2))) {
  252     errmsg("Set not active");
  253     return;
  254     }
  255     if ((getsetlength(gno1, set1) < 3) || (getsetlength(gno2, set2) < 3)) {
  256     errmsg("Set length < 3");
  257     return;
  258     }
  259     linearcset = nextset(cg);
  260     if (linearcset != (-1)) {
  261     activateset(cg, linearcset);
  262     setlength(cg, linearcset, (itmp = getsetlength(gno1, set1) + getsetlength(gno2, set2) - 1));
  263     linearconv(gety(gno2, set2), gety(gno1, set1), gety(cg, linearcset), getsetlength(gno2, set2), getsetlength(gno1, set1));
  264     xtmp = getx(cg, linearcset);
  265     for (i = 0; i < itmp; i++) {
  266         xtmp[i] = i;
  267     }
  268     sprintf(buf, "Linear convolution of set %d with set %d", set1, set2);
  269     setcomment(cg, linearcset, buf);
  270     }
  271 }
  272 
  273 /*
  274  * cross correlation/covariance
  275  */
  276 void do_xcor(int gno1, int set1, int gno2, int set2, int maxlag, int covar)
  277 {
  278     int xcorset, i, ierr, len, cg = get_cg();
  279     double *xtmp;
  280 
  281     if (!(is_set_active(gno1, set1) && is_set_active(gno2, set2))) {
  282     errmsg("Set not active");
  283     return;
  284     }
  285     
  286     len = getsetlength(gno1, set1);
  287     if (getsetlength(gno2, set2) != len) {
  288     errmsg("Sets must be of the same length");
  289     }
  290     if (len < 2) {
  291     errmsg("Set length < 2");
  292     return;
  293     }
  294     
  295     if (maxlag < 1 || maxlag > len) {
  296     errmsg("Lag incorrectly specified");
  297     return;
  298     }
  299     
  300     xcorset = nextset(cg);
  301     
  302     if (xcorset != (-1)) {
  303     char *fname;
  304         activateset(cg, xcorset);
  305     setlength(cg, xcorset, maxlag);
  306     if (covar) {
  307             fname = "covariance";
  308         } else {
  309             fname = "correlation";
  310         }
  311         if (set1 != set2) {
  312         sprintf(buf, "X-%s of G%d.S%d and G%d.S%d at maximum lag %d",
  313                     fname, gno1, set1, gno2, set2, maxlag);
  314     } else {
  315         sprintf(buf, "Auto-%s of G%d.S%d at maximum lag %d",
  316                     fname, gno1, set1, maxlag);
  317     }
  318     ierr = crosscorr(gety(gno1, set1), gety(gno2, set2), len,
  319                          maxlag, covar, gety(cg, xcorset));
  320     xtmp = getx(cg, xcorset);
  321     for (i = 0; i < maxlag; i++) {
  322         xtmp[i] = i;
  323     }
  324     setcomment(cg, xcorset, buf);
  325     }
  326 }
  327 
  328 
  329 /*
  330  * numerical integration
  331  */
  332 double do_int(int gno, int setno, int itype)
  333 {
  334     int intset;
  335     double sum = 0;
  336 
  337     if (!is_set_active(gno, setno)) {
  338     errmsg("Set not active");
  339     return 0.0;
  340     }
  341     if (getsetlength(gno, setno) < 3) {
  342     errmsg("Set length < 3");
  343     return 0.0;
  344     }
  345     if (itype == 0) {
  346     intset = nextset(gno);
  347     if (intset != (-1)) {
  348         activateset(gno, intset);
  349         setlength(gno, intset, getsetlength(gno, setno));
  350         sprintf(buf, "Cumulative sum of set %d", setno);
  351         sum = trapint(getx(gno, setno), gety(gno, setno), getx(gno, intset), gety(gno, intset), getsetlength(gno, setno));
  352         setcomment(gno, intset, buf);
  353     }
  354     } else {
  355     sum = trapint(getx(gno, setno), gety(gno, setno), NULL, NULL, getsetlength(gno, setno));
  356     }
  357     return sum;
  358 }
  359 
  360 /*
  361  * difference a set
  362  * itype means
  363  *  0 - forward
  364  *  1 - backward
  365  *  2 - centered difference
  366  */
  367 void do_differ(int gno, int setno, int itype)
  368 {
  369     int diffset;
  370 
  371     if (!is_set_active(gno, setno)) {
  372     errmsg("Set not active");
  373     return;
  374     }
  375     if (getsetlength(gno, setno) < 3) {
  376     errmsg("Set length < 3");
  377     return;
  378     }
  379     diffset = nextset(gno);
  380     if (diffset != (-1)) {
  381     activateset(gno, diffset);
  382     switch (itype) {
  383     case 0:
  384         sprintf(buf, "Forward difference of set %d", setno);
  385         setlength(gno, diffset, getsetlength(gno, setno) - 1);
  386         forwarddiff(getx(gno, setno), gety(gno, setno), getx(gno, diffset), gety(gno, diffset), getsetlength(gno, setno));
  387         break;
  388     case 1:
  389         sprintf(buf, "Backward difference of set %d", setno);
  390         setlength(gno, diffset, getsetlength(gno, setno) - 1);
  391         backwarddiff(getx(gno, setno), gety(gno, setno), getx(gno, diffset), gety(gno, diffset), getsetlength(gno, setno));
  392         break;
  393     case 2:
  394         sprintf(buf, "Centered difference of set %d", setno);
  395         setlength(gno, diffset, getsetlength(gno, setno) - 2);
  396         centereddiff(getx(gno, setno), gety(gno, setno), getx(gno, diffset), gety(gno, diffset), getsetlength(gno, setno));
  397         break;
  398     }
  399     setcomment(gno, diffset, buf);
  400     }
  401 }
  402 
  403 /*
  404  * seasonally difference a set
  405  */
  406 void do_seasonal_diff(int setno, int period)
  407 {
  408     int diffset;
  409 
  410     if (!is_set_active(get_cg(), setno)) {
  411     errmsg("Set not active");
  412     return;
  413     }
  414     if (getsetlength(get_cg(), setno) < 2) {
  415     errmsg("Set length < 2");
  416     return;
  417     }
  418     diffset = nextset(get_cg());
  419     if (diffset != (-1)) {
  420     activateset(get_cg(), diffset);
  421     setlength(get_cg(), diffset, getsetlength(get_cg(), setno) - period);
  422     seasonaldiff(getx(get_cg(), setno), gety(get_cg(), setno),
  423              getx(get_cg(), diffset), gety(get_cg(), diffset),
  424              getsetlength(get_cg(), setno), period);
  425     sprintf(buf, "Seasonal difference of set %d, period %d", setno, period);
  426     setcomment(get_cg(), diffset, buf);
  427     }
  428 }
  429 
  430 /*
  431  * regression with restrictions to region rno if rno >= 0
  432  */
  433 void do_regress(int gno, int setno, int ideg, int iresid, int rno, int invr, int fitset)
  434     /* 
  435      * gno, setno  - set to perform fit on
  436      * ideg   - degree of fit
  437      * irisid - 0 -> whole set, 1-> subset of setno
  438      * rno    - region number of subset
  439      * invr   - 1->invert region, 0 -> do nothing
  440      * fitset - set to which fitted function will be loaded
  441      *          Y values are computed at the x values in the set
  442      *          if -1 is specified, a set with the same x values as the
  443      *          one being fitted will be created and used
  444      */
  445 {
  446     int len, i, sdeg = ideg;
  447     int cnt = 0, fitlen = 0;
  448     double *x, *y, *xt = NULL, *yt = NULL, *xr = NULL, *yr = NULL;
  449     char buf[256];
  450     double c[20];   /* coefficients of fitted polynomial */
  451 
  452     if (!is_set_active(gno, setno)) {
  453         errmsg("Set not active");
  454         return;
  455     }
  456     len = getsetlength(gno, setno);
  457     x = getx(gno, setno);
  458     y = gety(gno, setno);
  459     if (rno == -1) {
  460         xt = x;
  461         yt = y;
  462     } else if (isactive_region(rno)) {
  463         if (!get_points_inregion(rno, invr, len, x, y, &cnt, &xt, &yt)) {
  464             if (cnt == 0) {
  465             errmsg("No points found in region, operation cancelled");
  466             }
  467             return;
  468         }
  469         len = cnt;
  470     } else {
  471         errmsg("Selected region is not active");
  472         return;
  473     }
  474     /*
  475      * first part for polynomials, second part for linear fits to transformed
  476      * data
  477      */
  478     if ((len < ideg && ideg <= 10) || (len < 2 && ideg > 10)) {
  479         errmsg("Too few points in set, operation cancelled");
  480         return;
  481     }
  482     /* determine is set provided or use abscissa from fitted set */
  483     if( fitset == -1 ) {                /* create set */
  484         if( (fitset = nextset(gno)) != -1 ) {
  485             activateset(gno, fitset);
  486             setlength(gno, fitset, len);
  487             fitlen = len;
  488             xr = getx(gno, fitset);
  489             for (i = 0; i < len; i++) {
  490                 xr[i] = xt[i];
  491             }   
  492             yr = gety(gno, fitset);
  493         }
  494     } else {                                    /* set has been provided */
  495         fitlen = getsetlength( gno, fitset );
  496         xr = getx(gno, fitset);
  497         yr = gety(gno, fitset);
  498     }
  499 
  500     /* transform data so that system has the form y' = A + B * x' */
  501     if (fitset != -1) {
  502     if (ideg == 12) {   /* y=A*x^B -> ln(y) = ln(A) + B * ln(x) */
  503         ideg = 1;
  504         for (i = 0; i < len; i++) {
  505             if (xt[i] <= 0.0) {
  506                 errmsg("One of X[i] <= 0.0");
  507                 return;
  508             }
  509             if (yt[i] <= 0.0) {
  510                 errmsg("One of Y[i] <= 0.0");
  511                 return;
  512             }
  513         }
  514         for (i = 0; i < len; i++) {
  515             xt[i] = log(xt[i]);
  516             yt[i] = log(yt[i]);
  517         }
  518         for( i=0; i<fitlen; i++ )
  519             if( xr[i] <= 0.0 ) {
  520                 errmsg("One of X[i] <= 0.0");
  521                 return;
  522             } 
  523         for( i=0; i<fitlen; i++ )
  524             xr[i] = log( xr[i] );
  525     } else if (ideg == 13) {   /*y=Aexp(B*x) -> ln(y) = ln(A) + B * x */
  526         ideg = 1;
  527         for (i = 0; i < len; i++) {
  528             if (yt[i] <= 0.0) {
  529                 errmsg("One of Y[i] <= 0.0");
  530                 return;
  531             }
  532         }
  533         for (i = 0; i < len; i++) {
  534             yt[i] = log(yt[i]);
  535         }
  536     } else if (ideg == 14) {    /* y = A + B * ln(x) */
  537         ideg = 1;
  538         for (i = 0; i < len; i++) {
  539             if (xt[i] <= 0.0) {
  540                 errmsg("One of X[i] <= 0.0");
  541                 return;
  542             }
  543         }
  544         for (i = 0; i < len; i++) {
  545             xt[i] = log(xt[i]); 
  546         }
  547         for( i=0; i<fitlen; i++ )
  548             if( xr[i] <= 0.0 ) {
  549                 errmsg("One of X[i] <= 0.0");
  550                 return;
  551             } 
  552         for( i=0; i<fitlen; i++ ){
  553             xr[i] = log( xr[i] );
  554         }     
  555     } else if (ideg == 15) {    /* y = 1/( A + B*x ) -> 1/y = a + B*x */
  556         ideg = 1;
  557         for (i = 0; i < len; i++) {
  558             if (yt[i] == 0.0) {
  559                 errmsg("One of Y[i] = 0.0");
  560                 return;
  561             }
  562         }
  563         for (i = 0; i < len; i++) {
  564             yt[i] = 1.0 / yt[i];
  565         }
  566     }
  567 
  568     if (fitcurve(xt, yt, len, ideg, c)) {
  569         killset(gno, fitset);
  570         goto bustout;
  571     }
  572 
  573     /* compute function at requested x ordinates */
  574     for( i=0; i<fitlen; i++ )
  575         yr[i] = leasev( c, ideg, xr[i] );
  576 
  577     sprintf(buf, "\nN.B. Statistics refer to the transformed model\n");
  578     /* apply inverse transform, output fitted function in human readable form */
  579     if( sdeg<11 ) {
  580         sprintf(buf, "\ny = %.5g %c %.5g * x", c[0], c[1]>=0?'+':'-', fabs(c[1]));
  581         for( i=2; i<=ideg; i++ )
  582             sprintf( buf+strlen(buf), " %c %.5g * x^%d", c[i]>=0?'+':'-', 
  583                                                             fabs(c[i]), i );
  584         strcat( buf, "\n" );
  585     } else if (sdeg == 12) {    /* ln(y) = ln(A) + b * ln(x) */
  586         sprintf( buf, "\ny = %.5g * x^%.5g\n", exp(c[0]), c[1] );
  587         for (i = 0; i < len; i++) {
  588             xt[i] = exp(xt[i]);
  589             yt[i] = exp(yt[i]);
  590         }
  591         for (i = 0; i < fitlen; i++){
  592             yr[i] = exp(yr[i]);
  593             xr[i] = exp(xr[i]);
  594         }
  595     } else if (sdeg == 13) {    /* ln(y) = ln(A) + B * x */
  596         sprintf( buf, "\ny = %.5g * exp( %.5g * x )\n", exp(c[0]), c[1] );
  597         for (i = 0; i < len; i++) {
  598             yt[i] = exp(yt[i]);
  599         }
  600         for (i = 0; i < fitlen; i++)
  601             yr[i] = exp(yr[i]);
  602     } else if (sdeg == 14) {    /* y = A + B * ln(x) */
  603         sprintf(buf, "\ny = %.5g %c %.5g * ln(x)\n", c[0], c[1]>=0?'+':'-',
  604                                             fabs(c[1]) );
  605         for (i = 0; i < len; i++) {
  606             xt[i] = exp(xt[i]);
  607         }
  608         for (i = 0; i < fitlen; i++)
  609             xr[i] = exp(xr[i]);
  610     } else if (sdeg == 15) {    /* y = 1/( A + B*x ) */
  611         sprintf( buf, "\ny = 1/(%.5g %c %.5g * x)\n", c[0], c[1]>=0?'+':'-',
  612                                                     fabs(c[1]) );
  613         for (i = 0; i < len; i++) {
  614             yt[i] = 1.0 / yt[i];
  615         }
  616         for (i = 0; i < fitlen; i++)
  617             yr[i] = 1.0 / yr[i];
  618     }
  619     stufftext(buf);
  620     sprintf(buf, "\nRegression of set %d results to set %d\n", setno, fitset);
  621     stufftext(buf);
  622     
  623     switch (iresid) {
  624     case 1:
  625         /* determine residual */
  626         for (i = 0; i < len; i++) {
  627             yr[i] = yt[i] - yr[i];
  628         }
  629         break;
  630     case 2:
  631         break;
  632     }
  633     sprintf(buf, "%d deg fit of set %d", ideg, setno);
  634     setcomment(gno, fitset, buf);
  635     }
  636     bustout:;
  637     if (rno >= 0 && cnt != 0) { /* had a region and allocated memory there */
  638     xfree(xt);
  639     xfree(yt);
  640     }
  641 }
  642 
  643 /*
  644  * running averages, medians, min, max, std. deviation
  645  */
  646 void do_runavg(int gno, int setno, int runlen, int runtype, int rno, int invr)
  647 {
  648     int runset;
  649     int len, cnt = 0;
  650     double *x, *y, *xt = NULL, *yt = NULL, *xr, *yr;
  651 
  652     if (!is_set_active(gno, setno)) {
  653     errmsg("Set not active");
  654     return;
  655     }
  656     if (runlen < 2) {
  657     errmsg("Length of running average < 2");
  658     return;
  659     }
  660     len = getsetlength(gno, setno);
  661     x = getx(gno, setno);
  662     y = gety(gno, setno);
  663     if (rno == -1) {
  664     xt = x;
  665     yt = y;
  666     } else if (isactive_region(rno)) {
  667     if (!get_points_inregion(rno, invr, len, x, y, &cnt, &xt, &yt)) {
  668         if (cnt == 0) {
  669         errmsg("No points found in region, operation cancelled");
  670         }
  671         return;
  672     }
  673     len = cnt;
  674     } else {
  675     errmsg("Selected region is not active");
  676     return;
  677     }
  678     if (runlen >= len) {
  679     errmsg("Length of running average > set length");
  680     goto bustout;
  681     }
  682     runset = nextset(gno);
  683     if (runset != (-1)) {
  684     activateset(gno, runset);
  685     setlength(gno, runset, len - runlen + 1);
  686     xr = getx(gno, runset);
  687     yr = gety(gno, runset);
  688     switch (runtype) {
  689     case 0:
  690         runavg(xt, yt, xr, yr, len, runlen);
  691         sprintf(buf, "%d-pt. avg. on set %d ", runlen, setno);
  692         break;
  693     case 1:
  694         runmedian(xt, yt, xr, yr, len, runlen);
  695         sprintf(buf, "%d-pt. median on set %d ", runlen, setno);
  696         break;
  697     case 2:
  698         runminmax(xt, yt, xr, yr, len, runlen, 0);
  699         sprintf(buf, "%d-pt. min on set %d ", runlen, setno);
  700         break;
  701     case 3:
  702         runminmax(xt, yt, xr, yr, len, runlen, 1);
  703         sprintf(buf, "%d-pt. max on set %d ", runlen, setno);
  704         break;
  705     case 4:
  706         runstddev(xt, yt, xr, yr, len, runlen);
  707         sprintf(buf, "%d-pt. std dev., set %d ", runlen, setno);
  708         break;
  709     }
  710     setcomment(gno, runset, buf);
  711     }
  712   bustout:;
  713     if (rno >= 0 && cnt != 0) { /* had a region and allocated memory there */
  714     xfree(xt);
  715     xfree(yt);
  716     }
  717 }
  718 
  719 /*
  720  * DFT by FFT or definition
  721  */
  722 void do_fourier(int gno, int setno, int fftflag, int load, int loadx, int invflag, int type, int wind)
  723 {
  724     int i, ilen;
  725     double *x, *y, *xx, *yy, delt, T;
  726     int i2 = 0, specset;
  727 
  728     if (!is_set_active(get_cg(), setno)) {
  729     errmsg("Set not active");
  730     return;
  731     }
  732     ilen = getsetlength(get_cg(), setno);
  733     if (ilen < 2) {
  734     errmsg("Set length < 2");
  735     return;
  736     }
  737     if (fftflag) {
  738     if ((i2 = ilog2(ilen)) <= 0) {
  739         errmsg("Set length not a power of 2");
  740         return;
  741     }
  742     }
  743     specset = nextset(get_cg());
  744     if (specset != -1) {
  745     activateset(get_cg(), specset);
  746     setlength(get_cg(), specset, ilen);
  747     xx = getx(get_cg(), specset);
  748     yy = gety(get_cg(), specset);
  749     x = getx(get_cg(), setno);
  750     y = gety(get_cg(), setno);
  751     copyx(get_cg(), setno, specset);
  752     copyy(get_cg(), setno, specset);
  753     if (wind != 0) {    /* apply data window if needed */
  754         apply_window(xx, yy, ilen, type, wind);
  755     }
  756     if (type == 0) {    /* real data */
  757         for (i = 0; i < ilen; i++) {
  758         xx[i] = yy[i];
  759         yy[i] = 0.0;
  760         }
  761     }
  762     if (fftflag) {
  763         fft(xx, yy, ilen, i2, invflag);
  764     } else {
  765         dft(xx, yy, ilen, invflag);
  766     }
  767     switch (load) {
  768     case 0:
  769         delt = (x[ilen-1] - x[0])/(ilen -1.0);
  770         T = (x[ilen - 1] - x[0]);
  771         xx = getx(get_cg(), specset);
  772         yy = gety(get_cg(), specset);
  773         for (i = 0; i < ilen / 2; i++) {
  774           /* carefully get amplitude of complex xform: 
  775          use abs(a[i]) + abs(a[-i]) except for zero term */
  776           if(i) yy[i] = hypot(xx[i], yy[i])+hypot(xx[ilen-i], yy[ilen-i]);
  777           else yy[i]=2*fabs(xx[i]);
  778         switch (loadx) {
  779         case 0:
  780             xx[i] = i;
  781             break;
  782         case 1:
  783             /* xx[i] = 2.0 * M_PI * i / ilen; */
  784             xx[i] = i / T;
  785             break;
  786         case 2:
  787             if (i == 0) {
  788             xx[i] = T + delt;   /* the mean */
  789             } else {
  790             /* xx[i] = (double) ilen / (double) i; */
  791             xx[i] = T / i;
  792             }
  793             break;
  794         }
  795         }
  796         setlength(get_cg(), specset, ilen / 2);
  797         break;
  798     case 1:
  799         delt = (x[ilen-1] - x[0])/(ilen -1.0);
  800         T = (x[ilen - 1] - x[0]);
  801         setlength(get_cg(), specset, ilen / 2);
  802         xx = getx(get_cg(), specset);
  803         yy = gety(get_cg(), specset);
  804         for (i = 0; i < ilen / 2; i++) {
  805         yy[i] = -atan2(yy[i], xx[i]);
  806         switch (loadx) {
  807         case 0:
  808             xx[i] = i;
  809             break;
  810         case 1:
  811             /* xx[i] = 2.0 * M_PI * i / ilen; */
  812             xx[i] = i / T;
  813             break;
  814         case 2:
  815             if (i == 0) {
  816             xx[i] = T + delt;
  817             } else {
  818             /* xx[i] = (double) ilen / (double) i; */
  819             xx[i] = T / i;
  820             }
  821             break;
  822         }
  823         }
  824         break;
  825     }
  826     if (fftflag) {
  827         sprintf(buf, "FFT of set %d", setno);
  828     } else {
  829         sprintf(buf, "DFT of set %d", setno);
  830     }
  831     setcomment(get_cg(), specset, buf);
  832     }
  833 }
  834 
  835 /*
  836  * Apply a window to a set, result goes to a new set.
  837  */
  838 void do_window(int setno, int type, int wind)
  839 {
  840     int ilen;
  841     double *xx, *yy;
  842     int specset;
  843 
  844     if (!is_set_active(get_cg(), setno)) {
  845     errmsg("Set not active");
  846     return;
  847     }
  848     ilen = getsetlength(get_cg(), setno);
  849     if (ilen < 2) {
  850     errmsg("Set length < 2");
  851     return;
  852     }
  853     specset = nextset(get_cg());
  854     if (specset != -1) {
  855     char *wtype[6];
  856     wtype[0] = "Triangular";
  857     wtype[1] = "Hanning";
  858     wtype[2] = "Welch";
  859     wtype[3] = "Hamming";
  860     wtype[4] = "Blackman";
  861     wtype[5] = "Parzen";
  862 
  863     activateset(get_cg(), specset);
  864     setlength(get_cg(), specset, ilen);
  865     xx = getx(get_cg(), specset);
  866     yy = gety(get_cg(), specset);
  867     copyx(get_cg(), setno, specset);
  868     copyy(get_cg(), setno, specset);
  869     if (wind != 0) {
  870         apply_window(xx, yy, ilen, type, wind);
  871         sprintf(buf, "%s windowed set %d", wtype[wind - 1], setno);
  872     } else {        /* shouldn't happen */
  873     }
  874     setcomment(get_cg(), specset, buf);
  875     }
  876 }
  877 
  878 void apply_window(double *xx, double *yy, int ilen, int type, int wind)
  879 {
  880     int i;
  881 
  882     for (i = 0; i < ilen; i++) {
  883     switch (wind) {
  884     case 1:     /* triangular */
  885         if (type != 0) {
  886         xx[i] *= 1.0 - fabs((i - 0.5 * (ilen - 1.0)) / (0.5 * (ilen - 1.0)));
  887         }
  888         yy[i] *= 1.0 - fabs((i - 0.5 * (ilen - 1.0)) / (0.5 * (ilen - 1.0)));
  889 
  890         break;
  891     case 2:     /* Hanning */
  892         if (type != 0) {
  893         xx[i] = xx[i] * (0.5 - 0.5 * cos(2.0 * M_PI * i / (ilen - 1.0)));
  894         }
  895         yy[i] = yy[i] * (0.5 - 0.5 * cos(2.0 * M_PI * i / (ilen - 1.0)));
  896         break;
  897     case 3:     /* Welch (from Numerical Recipes) */
  898         if (type != 0) {
  899         xx[i] *= 1.0 - pow((i - 0.5 * (ilen - 1.0)) / (0.5 * (ilen + 1.0)), 2.0);
  900         }
  901         yy[i] *= 1.0 - pow((i - 0.5 * (ilen - 1.0)) / (0.5 * (ilen + 1.0)), 2.0);
  902         break;
  903     case 4:     /* Hamming */
  904         if (type != 0) {
  905         xx[i] = xx[i] * (0.54 - 0.46 * cos(2.0 * M_PI * i / (ilen - 1.0)));
  906         }
  907         yy[i] = yy[i] * (0.54 - 0.46 * cos(2.0 * M_PI * i / (ilen - 1.0)));
  908         break;
  909     case 5:     /* Blackman */
  910         if (type != 0) {
  911         xx[i] = xx[i] * (0.42 - 0.5 * cos(2.0 * M_PI * i / (ilen - 1.0)) + 0.08 * cos(4.0 * M_PI * i / (ilen - 1.0)));
  912         }
  913         yy[i] = yy[i] * (0.42 - 0.5 * cos(2.0 * M_PI * i / (ilen - 1.0)) + 0.08 * cos(4.0 * M_PI * i / (ilen - 1.0)));
  914         break;
  915     case 6:     /* Parzen (from Numerical Recipes) */
  916         if (type != 0) {
  917         xx[i] *= 1.0 - fabs((i - 0.5 * (ilen - 1)) / (0.5 * (ilen + 1)));
  918         }
  919         yy[i] *= 1.0 - fabs((i - 0.5 * (ilen - 1)) / (0.5 * (ilen + 1)));
  920         break;
  921     }
  922     }
  923 }
  924 
  925 
  926 /*
  927  * histograms
  928  */
  929 int do_histo(int fromgraph, int fromset, int tograph, int toset,
  930           double *bins, int nbins, int cumulative, int normalize)
  931 {
  932     int i, ndata;
  933     int *hist;
  934     double *x, *y, *data;
  935     plotarr p;
  936     
  937     if (!is_set_active(fromgraph, fromset)) {
  938     errmsg("Set not active");
  939     return RETURN_FAILURE;
  940     }
  941     if (nbins <= 0) {
  942     errmsg("Number of bins <= 0");
  943     return RETURN_FAILURE;
  944     }
  945     if (toset == SET_SELECT_NEXT) {
  946     toset = nextset(tograph);
  947     }
  948     if (!is_valid_setno(tograph, toset)) {
  949     errmsg("Can't activate destination set");
  950         return RETURN_FAILURE;
  951     }
  952     
  953     ndata = getsetlength(fromgraph, fromset);
  954     data = gety(fromgraph, fromset);
  955     
  956     hist = xmalloc(nbins*SIZEOF_INT);
  957     if (hist == NULL) {
  958         errmsg("xmalloc failed in do_histo()");
  959         return RETURN_FAILURE;
  960     }
  961 
  962     if (histogram(ndata, data, nbins, bins, hist) == RETURN_FAILURE){
  963         xfree(hist);
  964         return RETURN_FAILURE;
  965     }
  966     
  967     activateset(tograph, toset);
  968     setlength(tograph, toset, nbins + 1);
  969     x = getx(tograph, toset);
  970     y = gety(tograph, toset);
  971     
  972     x[0] = bins[0];
  973     y[0] = 0.0;
  974     for (i = 1; i < nbins + 1; i++) {
  975         x[i] = bins[i];
  976         y[i] = hist[i - 1];
  977         if (cumulative) {
  978             y[i] += y[i - 1];
  979         }
  980     }
  981     
  982     if (normalize) {
  983         for (i = 1; i < nbins + 1; i++) {
  984             double factor;
  985             if (cumulative) {
  986                 factor = 1.0/ndata;
  987             } else {
  988                 factor = 1.0/((bins[i] - bins[i - 1])*ndata);
  989             }
  990             y[i] *= factor;
  991         }
  992     }
  993     
  994     xfree(hist);
  995 
  996     get_graph_plotarr(tograph, toset, &p);
  997     p.sym = SYM_NONE;
  998     p.linet = LINE_TYPE_LEFTSTAIR;
  999     p.dropline = TRUE;
 1000     p.baseline = FALSE;
 1001     p.baseline_type = BASELINE_TYPE_0;
 1002     p.lines = 1;
 1003     p.symlines = 1;
 1004     sprintf(p.comments, "Histogram from G%d.S%d", fromgraph, fromset);
 1005     set_graph_plotarr(tograph, toset, &p);
 1006 
 1007     return RETURN_SUCCESS;
 1008 }
 1009 
 1010 int histogram(int ndata, double *data, int nbins, double *bins, int *hist)
 1011 {
 1012     int i, j, bsign;
 1013     double minval, maxval;
 1014     
 1015     if (nbins < 1) {
 1016         errmsg("Number of bins < 1");
 1017         return RETURN_FAILURE;
 1018     }
 1019     
 1020     bsign = monotonicity(bins, nbins + 1, TRUE);
 1021     if (bsign == 0) {
 1022         errmsg("Non-monotonic bins");
 1023         return RETURN_FAILURE;
 1024     }
 1025     
 1026     for (i = 0; i < nbins; i++) {
 1027         hist[i] = 0;
 1028     }
 1029     
 1030     /* TODO: binary search */
 1031     for (i = 0; i < ndata; i++) {
 1032         for (j = 0; j < nbins; j++) {
 1033             if (bsign > 0) {
 1034                 minval = bins[j];
 1035                 maxval = bins[j + 1];
 1036             } else {
 1037                 minval = bins[j + 1];
 1038                 maxval = bins[j];
 1039             }
 1040             if (data[i] >= minval && data[i] <= maxval) {
 1041                 hist[j] += 1;
 1042                 break;
 1043             }
 1044         }
 1045     }
 1046     return RETURN_SUCCESS;
 1047 }
 1048 
 1049 
 1050 /*
 1051  * sample a set, by start/step or logical expression
 1052  */
 1053 void do_sample(int setno, int typeno, char *exprstr, int startno, int stepno)
 1054 {
 1055     int len, npts = 0, i, resset;
 1056     double *x, *y;
 1057     int reslen;
 1058     double *result;
 1059     int gno = get_cg();
 1060 
 1061     if (!is_set_active(gno, setno)) {
 1062     errmsg("Set not active");
 1063     return;
 1064     }
 1065 
 1066     len = getsetlength(gno, setno);
 1067 
 1068     resset = nextset(gno);
 1069     if (resset < 0) {
 1070     return;
 1071     }
 1072 
 1073     x = getx(gno, setno);
 1074     y = gety(gno, setno);
 1075 
 1076     if (typeno == 0) {
 1077     if (len <= 2) {
 1078         errmsg("Set has <= 2 points");
 1079         return;
 1080     }
 1081     if (startno < 1) {
 1082         errmsg("Start point < 1 (locations in sets are numbered starting from 1)");
 1083         return;
 1084     }
 1085     if (stepno < 1) {
 1086         errmsg("Step < 1");
 1087         return;
 1088     }
 1089     for (i = startno - 1; i < len; i += stepno) {
 1090         add_point(gno, resset, x[i], y[i]);
 1091         npts++;
 1092     }
 1093     sprintf(buf, "Sample, %d, %d set #%d", startno, stepno, setno);
 1094     } else {
 1095         if (set_parser_setno(gno, setno) != RETURN_SUCCESS) {
 1096         errmsg("Bad set");
 1097             killset(gno, resset);
 1098         return;
 1099         }
 1100         if (v_scanner(exprstr, &reslen, &result) != RETURN_SUCCESS) {
 1101             killset(gno, resset);
 1102         return;
 1103         }
 1104         if (reslen != len) {
 1105         errmsg("Internal error");
 1106             killset(gno, resset);
 1107         return;
 1108         }
 1109 
 1110         npts = 0;
 1111     sprintf(buf, "Sample from %d, using '%s'", setno, exprstr);
 1112     for (i = 0; i < len; i++) {
 1113         if ((int) rint(result[i])) {
 1114         add_point(gno, resset, x[i], y[i]);
 1115         npts++;
 1116         }
 1117     }
 1118         xfree(result);
 1119     }
 1120     if (npts > 0) {
 1121     setcomment(gno, resset, buf);
 1122     }
 1123 }
 1124 
 1125 #define prune_xconv(res,x,xtype)    \
 1126     switch (deltatypeno) {      \
 1127     case PRUNE_VIEWPORT:        \
 1128     res = xy_xconv(x);          \
 1129     break;              \
 1130     case PRUNE_WORLD:           \
 1131     switch (xtype) {        \
 1132     case PRUNE_LIN:         \
 1133         res = x;            \
 1134         break;          \
 1135     case PRUNE_LOG:         \
 1136         res = log(x);       \
 1137         break;          \
 1138     }               \
 1139     }
 1140 
 1141 #define prune_yconv(res,y,ytype)    \
 1142     switch (deltatypeno) {      \
 1143     case PRUNE_VIEWPORT:        \
 1144     res = xy_yconv(y);          \
 1145     break;              \
 1146     case PRUNE_WORLD:           \
 1147     switch (ytype) {        \
 1148     case PRUNE_LIN:         \
 1149         res = y;            \
 1150         break;          \
 1151     case PRUNE_LOG:         \
 1152         res = log(y);       \
 1153         break;          \
 1154     }               \
 1155     }
 1156 
 1157 /*
 1158  * Prune data
 1159  */
 1160 void do_prune(int setno, int typeno, int deltatypeno, float deltax, float deltay, int dxtype, int dytype)
 1161 {
 1162     int len, npts = 0, d, i, j, k, drop, resset;
 1163     double *x, *y, *resx, *resy, xtmp, ytmp, ddx = 0.0, ddy = 0.0;
 1164     double xj = 0.0, xjm = 0.0, xjp = 0.0, yj = 0.0, yjm = 0.0, yjp = 0.0;
 1165 
 1166     if (!is_set_active(get_cg(), setno)) {
 1167         errmsg("Set not active");
 1168         return;
 1169     }
 1170     len = getsetlength(get_cg(), setno);
 1171     if (len <= 2) {
 1172     errmsg("Set has <= 2 points");
 1173     return;
 1174     }
 1175     x = getx(get_cg(), setno);
 1176     y = gety(get_cg(), setno);
 1177     switch (typeno) {
 1178     case PRUNE_CIRCLE:
 1179     case PRUNE_ELLIPSE:
 1180     case PRUNE_RECTANGLE:
 1181     deltax = fabs(deltax);
 1182     if (deltax == 0)
 1183         return;
 1184     break;
 1185     }
 1186     switch (typeno) {
 1187     case PRUNE_CIRCLE:
 1188     deltay = deltax;
 1189     break;
 1190     case PRUNE_ELLIPSE:
 1191     case PRUNE_RECTANGLE:
 1192     case PRUNE_INTERPOLATION:
 1193     deltay = fabs(deltay);
 1194     if (deltay == 0)
 1195         return;
 1196     break;
 1197     }
 1198     if (deltatypeno == PRUNE_WORLD) {
 1199     if (dxtype == PRUNE_LOG && deltax < 1.0) {
 1200         deltax = 1.0 / deltax;
 1201     }
 1202     if (dytype == PRUNE_LOG && deltay < 1.0) {
 1203         deltay = 1.0 / deltay;
 1204     }
 1205     }
 1206     resset = nextset(get_cg());
 1207     if (resset < 0) {
 1208         return;
 1209     }
 1210     add_point(get_cg(), resset, x[0], y[0]);
 1211     npts++;
 1212     resx = getx(get_cg(), resset);
 1213     resy = gety(get_cg(), resset);
 1214     switch (typeno) {
 1215     case PRUNE_CIRCLE:
 1216     case PRUNE_ELLIPSE:
 1217     for (i = 1; i < len; i++) {
 1218         xtmp = x[i];
 1219         ytmp = y[i];
 1220         drop = FALSE;
 1221         for (j = npts - 1; j >= 0 && drop == FALSE; j--) {
 1222         switch (deltatypeno) {
 1223         case PRUNE_VIEWPORT:
 1224             ddx = (xy_xconv(xtmp) - xy_xconv(resx[j])) / deltax;
 1225             if (fabs(ddx) < 1.0) {
 1226             ddy = (xy_yconv(ytmp) - xy_yconv(resy[j])) / deltay;
 1227             if (ddx * ddx + ddy * ddy < 1.0) {
 1228                 drop = TRUE;
 1229             }
 1230             }
 1231             break;
 1232         case PRUNE_WORLD:
 1233             switch (dxtype) {
 1234             case PRUNE_LIN:
 1235             ddx = (xtmp - resx[j]) / deltax;
 1236             break;
 1237             case PRUNE_LOG:
 1238             ddx = (xtmp / resx[j]);
 1239             if (ddx < 1.0) {
 1240                 ddx = 1.0 / ddx;
 1241             }
 1242             ddx /= deltax;
 1243             break;
 1244             }
 1245             if (fabs(ddx) < 1.0) {
 1246             switch (dytype) {
 1247             case PRUNE_LIN:
 1248                 ddy = (ytmp - resy[j]) / deltay;
 1249                 break;
 1250             case PRUNE_LOG:
 1251                 ddy = (ytmp / resy[j]);
 1252                 if (ddy < 1.0) {
 1253                 ddy = 1.0 / ddy;
 1254                 }
 1255                 ddy /= deltay;
 1256                 break;
 1257             }
 1258             if (ddx * ddx + ddy * ddy < 1.0) {
 1259                 drop = TRUE;
 1260             }
 1261             }
 1262             break;
 1263         }
 1264         }
 1265         if (drop == FALSE) {
 1266         add_point(get_cg(), resset, xtmp, ytmp);
 1267         npts++;
 1268         resx = getx(get_cg(), resset);
 1269         resy = gety(get_cg(), resset);
 1270         }
 1271     }
 1272     sprintf(buf, "Prune from %d, %s dx = %g dy = %g", setno, 
 1273         (typeno == 0) ? "Circle" : "Ellipse", deltax, deltay);
 1274     break;
 1275     case PRUNE_RECTANGLE:
 1276     for (i = 1; i < len; i++) {
 1277         xtmp = x[i];
 1278         ytmp = y[i];
 1279         drop = FALSE;
 1280         for (j = npts - 1; j >= 0 && drop == FALSE; j--) {
 1281         switch (deltatypeno) {
 1282         case PRUNE_VIEWPORT:
 1283             ddx = fabs(xy_xconv(xtmp) - xy_xconv(resx[j]));
 1284             if (ddx < deltax) {
 1285             ddy = fabs(xy_yconv(ytmp) - xy_yconv(resy[j]));
 1286             if (ddy < deltay) {
 1287                 drop = TRUE;
 1288             }
 1289             }
 1290             break;
 1291         case PRUNE_WORLD:
 1292             switch (dxtype) {
 1293             case PRUNE_LIN:
 1294             ddx = fabs(xtmp - resx[j]);
 1295             break;
 1296             case PRUNE_LOG:
 1297             ddx = (xtmp / resx[j]);
 1298             if (ddx < 1.0) {
 1299                 ddx = 1.0 / ddx;
 1300             }
 1301             break;
 1302             }
 1303             if (ddx < deltax) {
 1304             switch (dytype) {
 1305             case PRUNE_LIN:
 1306                 ddy = fabs(ytmp - resy[j]);
 1307                 break;
 1308             case PRUNE_LOG:
 1309                 ddy = (ytmp / resy[j]);
 1310                 if (ddy < 1.0) {
 1311                 ddy = 1.0 / ddy;
 1312                 }
 1313                 break;
 1314             }
 1315             if (ddy < deltay) {
 1316                 drop = TRUE;
 1317             }
 1318             }
 1319             break;
 1320         }
 1321         }
 1322         if (drop == FALSE) {
 1323         add_point(get_cg(), resset, xtmp, ytmp);
 1324         npts++;
 1325         resx = getx(get_cg(), resset);
 1326         resy = gety(get_cg(), resset);
 1327         }
 1328     }
 1329     sprintf(buf, "Prune from %d, %s dx = %g dy = %g", setno, 
 1330         "Rectangle", deltax, deltay);
 1331     break;
 1332     case PRUNE_INTERPOLATION:
 1333     k = 0;
 1334     prune_xconv(xjm, x[0], dxtype);
 1335     prune_yconv(yjm, y[0], dytype);
 1336     while (k < len - 2) {
 1337         d = 1;
 1338         i = k + d + 1;
 1339         drop = TRUE;
 1340         while (TRUE) {
 1341         prune_xconv(xjp, x[i], dxtype);
 1342         prune_yconv(yjp, y[i], dytype);
 1343         for (j = k + 1; j < i && drop == TRUE; j++) {
 1344             prune_xconv(xj, x[j], dxtype);
 1345             prune_yconv(yj, y[j], dytype);
 1346             if (xjp == xjm) {
 1347             ytmp = 0.5 * (yjp + yjm);
 1348             } else {
 1349             ytmp = (yjp*(xj-xjm)+yjm*(xjp-xj))/(xjp-xjm);
 1350             }
 1351             switch (deltatypeno) {
 1352             case PRUNE_VIEWPORT:
 1353             ddy = fabs(ytmp - yj);
 1354             break;
 1355             case PRUNE_WORLD:
 1356             switch (dytype) {
 1357             case PRUNE_LIN:
 1358                 ddy = fabs(ytmp - yj);
 1359                 break;
 1360             case PRUNE_LOG:
 1361                 ddy = exp(fabs(ytmp - yj));
 1362                 break;
 1363             }
 1364             }
 1365             if (ddy > deltay) {
 1366             drop = FALSE;
 1367             }
 1368         }
 1369         if (drop == FALSE || i == len - 1) {
 1370             break;
 1371         }
 1372         d *= 2;
 1373         i = k + d + 1;
 1374         if (i >= len) {
 1375             i = len - 1;
 1376         }
 1377         }
 1378         if (drop == FALSE) {
 1379         i = k + 1;
 1380         drop = TRUE;
 1381         while (d > 1) {
 1382             d /= 2;
 1383             i += d;
 1384             prune_xconv(xjp, x[i], dxtype);
 1385             prune_yconv(yjp, y[i], dytype);
 1386             drop = TRUE;
 1387             for (j = k + 1; j < i && drop == TRUE; j++) {
 1388             prune_xconv(xj, x[j], dxtype);
 1389             prune_yconv(yj, y[j], dytype);
 1390             ytmp = (yjp*(xj-xjm)+yjm*(xjp-xj))/(xjp-xjm);
 1391             switch (deltatypeno) {
 1392             case PRUNE_VIEWPORT:
 1393                 ddy = fabs(ytmp - yj);
 1394                 break;
 1395             case PRUNE_WORLD:
 1396                 switch (dytype) {
 1397                 case PRUNE_LIN:
 1398                 ddy = fabs(ytmp - yj);
 1399                 break;
 1400                 case PRUNE_LOG:
 1401                 ddy = exp(fabs(ytmp - yj));
 1402                 break;
 1403                 }
 1404             }
 1405             if (ddy > deltay) {
 1406                 drop = FALSE;
 1407             }
 1408             }
 1409             if (drop == FALSE) {
 1410             i -= d;
 1411             }
 1412         }
 1413         }
 1414         k = i;
 1415         prune_xconv(xjm, x[k], dxtype);
 1416         prune_yconv(yjm, y[k], dytype);
 1417         add_point(get_cg(), resset, x[k], y[k]);
 1418         npts++;
 1419         resx = getx(get_cg(), resset);
 1420         resy = gety(get_cg(), resset);
 1421     }
 1422     if (k == len - 2) {
 1423         add_point(get_cg(), resset, x[len-1], y[len-1]);
 1424         npts++;
 1425     }
 1426     sprintf(buf, "Prune from %d, %s dy = %g", setno, 
 1427         "Interpolation", deltay);
 1428     break;
 1429     }
 1430     setcomment(get_cg(), resset, buf);
 1431 }
 1432 
 1433 int get_points_inregion(int rno, int invr, int len, double *x, double *y, int *cnt, double **xt, double **yt)
 1434 {
 1435     int i, clen = 0;
 1436     double *xtmp, *ytmp;
 1437     *cnt = 0;
 1438     if (isactive_region(rno)) {
 1439     for (i = 0; i < len; i++) {
 1440         if (invr) {
 1441         if (!inregion(rno, x[i], y[i])) {
 1442             clen++;
 1443         }
 1444         } else {
 1445         if (inregion(rno, x[i], y[i])) {
 1446             clen++;
 1447         }
 1448         }
 1449     }
 1450     if (clen == 0) {
 1451         return 0;
 1452     }
 1453     xtmp = (double *) xcalloc(clen, SIZEOF_DOUBLE);
 1454     if (xtmp == NULL) {
 1455         return 0;
 1456     }
 1457     ytmp = (double *) xcalloc(clen, SIZEOF_DOUBLE);
 1458     if (ytmp == NULL) {
 1459         xfree(xtmp);
 1460         return 0;
 1461     }
 1462     clen = 0;
 1463     for (i = 0; i < len; i++) {
 1464         if (invr) {
 1465         if (!inregion(rno, x[i], y[i])) {
 1466             xtmp[clen] = x[i];
 1467             ytmp[clen] = y[i];
 1468             clen++;
 1469         }
 1470         } else {
 1471         if (inregion(rno, x[i], y[i])) {
 1472             xtmp[clen] = x[i];
 1473             ytmp[clen] = y[i];
 1474             clen++;
 1475         }
 1476         }
 1477     }
 1478     } else {
 1479     return 0;
 1480     }
 1481     *cnt = clen;
 1482     *xt = xtmp;
 1483     *yt = ytmp;
 1484     return 1;
 1485 }
 1486 
 1487 int interpolate(double *mesh, double *yint, int meshlen,
 1488     double *x, double *y, int len, int method)
 1489 {
 1490     double *b, *c, *d;
 1491     double dx;
 1492     int i, ifound;
 1493     int m;
 1494 
 1495     /* For linear interpolation, non-strict monotonicity is fine */
 1496     m = monotonicity(x, len, method == INTERP_LINEAR ? FALSE:TRUE);
 1497     if (m == 0) {
 1498         errmsg("Can't interpolate a set with non-monotonic abscissas");
 1499         return RETURN_FAILURE;
 1500     }
 1501 
 1502     switch (method) {
 1503     case INTERP_SPLINE:
 1504     case INTERP_ASPLINE:
 1505         b = xcalloc(len, SIZEOF_DOUBLE);
 1506         c = xcalloc(len, SIZEOF_DOUBLE);
 1507         d = xcalloc(len, SIZEOF_DOUBLE);
 1508         if (b == NULL || c == NULL || d == NULL) {
 1509             xfree(b);
 1510             xfree(c);
 1511             xfree(d);
 1512             return RETURN_FAILURE;
 1513         }
 1514         if (method == INTERP_ASPLINE){
 1515             /* Akima spline */
 1516             aspline(len, x, y, b, c, d);
 1517         } else {
 1518             /* Plain cubic spline */
 1519             spline(len, x, y, b, c, d);
 1520         }
 1521 
 1522     seval(mesh, yint, meshlen, x, y, b, c, d, len);
 1523 
 1524         xfree(b);
 1525         xfree(c);
 1526         xfree(d);
 1527         break;
 1528     default:
 1529         /* linear interpolation */
 1530 
 1531         for (i = 0; i < meshlen; i++) {
 1532             ifound = find_span_index(x, len, m, mesh[i]);
 1533             if (ifound < 0) {
 1534                 ifound = 0;
 1535             } else if (ifound > len - 2) {
 1536                 ifound = len - 2;
 1537             }
 1538             dx = x[ifound + 1] - x[ifound];
 1539             if (dx != 0.0) {
 1540                 yint[i] = y[ifound] + (mesh[i] - x[ifound])*
 1541                     ((y[ifound + 1] - y[ifound])/dx);
 1542             } else {
 1543                 yint[i] = (y[ifound] + y[ifound + 1])/2;
 1544             }
 1545         }
 1546         break;
 1547     }
 1548     
 1549     return RETURN_SUCCESS;
 1550 }
 1551 
 1552 /* interpolate a set at abscissas from mesh
 1553  * method - type of spline (or linear interpolation)
 1554  * if strict is set, perform interpolation only within source set bounds
 1555  * (i.e., no extrapolation)
 1556  */
 1557 int do_interp(int gno_src, int setno_src, int gno_dest, int setno_dest,
 1558     double *mesh, int meshlen, int method, int strict)
 1559 {
 1560     int len, n, ncols;
 1561     double *x, *xint;
 1562     char *s;
 1563     
 1564     if (!is_valid_setno(gno_src, setno_src)) {
 1565     errmsg("Interpolated set not active");
 1566     return RETURN_FAILURE;
 1567     }
 1568     if (mesh == NULL || meshlen < 1) {
 1569         errmsg("NULL sampling mesh");
 1570         return RETURN_FAILURE;
 1571     }
 1572     
 1573     len = getsetlength(gno_src, setno_src);
 1574     ncols = dataset_cols(gno_src, setno_src);
 1575 
 1576     if (setno_dest == SET_SELECT_NEXT) {
 1577     setno_dest = nextset(gno_dest);
 1578     }
 1579     if (!is_valid_setno(gno_dest, setno_dest)) {
 1580     errmsg("Can't activate destination set");
 1581     return RETURN_FAILURE;
 1582     }
 1583 
 1584     if (dataset_cols(gno_dest, setno_dest) != ncols) {
 1585         copyset(gno_src, setno_src, gno_dest, setno_dest);
 1586     }
 1587     
 1588     setlength(gno_dest, setno_dest, meshlen);
 1589     activateset(gno_dest, setno_dest);
 1590     
 1591     x = getcol(gno_src, setno_src, DATA_X);
 1592     for (n = 1; n < ncols; n++) {
 1593         int res;
 1594         double *y, *yint;
 1595         
 1596         y    = getcol(gno_src, setno_src, n);
 1597         yint = getcol(gno_dest, setno_dest, n);
 1598         
 1599         res = interpolate(mesh, yint, meshlen, x, y, len, method);
 1600         if (res != RETURN_SUCCESS) {
 1601             killset(gno_dest, setno_dest);
 1602             return RETURN_FAILURE;
 1603         }
 1604     }
 1605 
 1606     xint = getcol(gno_dest, setno_dest, DATA_X);
 1607     memcpy(xint, mesh, meshlen*SIZEOF_DOUBLE);
 1608 
 1609     if (strict) {
 1610         double xmin, xmax;
 1611         int i, imin, imax;
 1612         minmax(x, len, &xmin, &xmax, &imin, &imax);
 1613         for (i = meshlen - 1; i >= 0; i--) {
 1614             if (xint[i] < xmin || xint[i] > xmax) {
 1615                 del_point(gno_dest, setno_dest, i);
 1616             }
 1617         }
 1618     }
 1619     
 1620     switch (method) {
 1621     case INTERP_SPLINE:
 1622         s = "cubic spline";
 1623         break;
 1624     case INTERP_ASPLINE:
 1625         s = "Akima spline";
 1626         break;
 1627     default:
 1628         s = "linear interpolation";
 1629         break;
 1630     }
 1631     sprintf(buf, "Interpolated from G%d.S%d using %s", gno_src, setno_src, s);
 1632     setcomment(gno_dest, setno_dest, buf);
 1633     
 1634     return RETURN_SUCCESS;
 1635 }
 1636 
 1637 int get_restriction_array(int gno, int setno,
 1638     int rtype, int negate, char **rarray)
 1639 {
 1640     int i, n, regno;
 1641     double *x, *y;
 1642     world w;
 1643     WPoint wp;
 1644     
 1645     if (rtype == RESTRICT_NONE) {
 1646         *rarray = NULL;
 1647         return RETURN_SUCCESS;
 1648     }
 1649     
 1650     n = getsetlength(gno, setno);
 1651     if (n <= 0) {
 1652         *rarray = NULL;
 1653         return RETURN_FAILURE;
 1654     }
 1655     
 1656     *rarray = xmalloc(n*SIZEOF_CHAR);
 1657     if (*rarray == NULL) {
 1658         return RETURN_FAILURE;
 1659     }
 1660     
 1661     x = getcol(gno, setno, DATA_X);
 1662     y = getcol(gno, setno, DATA_Y);
 1663     
 1664     switch (rtype) {
 1665     case RESTRICT_REG0:
 1666     case RESTRICT_REG1:
 1667     case RESTRICT_REG2:
 1668     case RESTRICT_REG3:
 1669     case RESTRICT_REG4:
 1670         regno = rtype - RESTRICT_REG0;
 1671         for (i = 0; i < n; i++) {
 1672             (*rarray)[i] = inregion(regno, x[i], y[i]) ? !negate : negate;
 1673         }
 1674         break;
 1675     case RESTRICT_WORLD:
 1676         get_graph_world(gno, &w);
 1677         for (i = 0; i < n; i++) {
 1678             wp.x = x[i];
 1679             wp.y = y[i];
 1680             (*rarray)[i] = is_wpoint_inside(&wp, &w) ? !negate : negate;
 1681         }
 1682         break;
 1683     default:
 1684         errmsg("Internal error in get_restriction_array()");
 1685         XCFREE(*rarray);
 1686         return RETURN_FAILURE;
 1687         break;
 1688     }
 1689     return RETURN_SUCCESS;
 1690 }
 1691 
 1692 int monotonicity(double *array, int len, int strict)
 1693 {
 1694     int i;
 1695     int s0, s1;
 1696     
 1697     if (len < 2) {
 1698         errmsg("Monotonicity of an array of length < 2 is meaningless");
 1699         return 0;
 1700     }
 1701     
 1702     s0 = sign(array[1] - array[0]);
 1703     for (i = 2; i < len; i++) {
 1704         s1 = sign(array[i] - array[i - 1]);
 1705         if (s1 != s0) {
 1706             if (strict) {
 1707                 return 0;
 1708             } else if (s0 == 0) {
 1709                 s0 = s1;
 1710             } else if (s1 != 0) {
 1711                 return 0;
 1712             }
 1713         }
 1714     }
 1715     
 1716     return s0;
 1717 }
 1718 
 1719 int find_span_index(double *array, int len, int m, double x)
 1720 {
 1721     int ind, low = 0, high = len - 1;
 1722     
 1723     if (len < 2 || m == 0) {
 1724         errmsg("find_span_index() called with a non-monotonic array");
 1725         return -2;
 1726     } else if (m > 0) {
 1727         /* ascending order */
 1728         if (x < array[0]) {
 1729             return -1;
 1730         } else if (x > array[len - 1]) {
 1731             return len - 1;
 1732         } else {
 1733             while (low <= high) {
 1734             ind = (low + high) / 2;
 1735             if (x < array[ind]) {
 1736                 high = ind - 1;
 1737             } else if (x > array[ind + 1]) {
 1738                 low = ind + 1;
 1739             } else {
 1740                 return ind;
 1741             }
 1742             }
 1743         }
 1744     } else {
 1745         /* descending order */
 1746         if (x > array[0]) {
 1747             return -1;
 1748         } else if (x < array[len - 1]) {
 1749             return len - 1;
 1750         } else {
 1751             while (low <= high) {
 1752             ind = (low + high) / 2;
 1753             if (x > array[ind]) {
 1754                 high = ind - 1;
 1755             } else if (x < array[ind + 1]) {
 1756                 low = ind + 1;
 1757             } else {
 1758                 return ind;
 1759             }
 1760             }
 1761         }
 1762     }
 1763 
 1764     /* should never happen */
 1765     errmsg("internal error in find_span_index()");
 1766     return -2;
 1767 }