grace_fit.c

Go to the documentation of this file.

/* 
Copyright 2021 Arizona Board of Regents on behalf of the University of Arizona.

For commercial uses, to obtain a license to sell and or sublicense copies of
the Software please contact the University of Arizona at Tech Launch Arizona:
lewish@tla.arizona.edu.

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction for non-commercial use educational and
research purposes, including without limitation the rights to use, copy,
modify, merge, publish, distribute, and to permit persons to whom the Software
is furnished to do so, subject to the following conditions:

The entire above copyright notice and this permission notice shall be included
in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.

1. Additional Required Provisions

1.1.  Arbitration. The parties agree that if a dispute arises between them
concerning this Agreement, the parties may be required to submit the matter to
arbitration pursuant to Arizona law.

1.2.  Applicable Law and Venue. This Agreement shall be interpreted pursuant
to the laws of the State of Arizona. Any arbitration or litigation between the
Parties shall be conducted in Pima County, ARIZONA, and LICENSEE hereby
submits to venue and jurisdiction in Pima County, ARIZONA.

1.3.  Non-Discrimination. The Parties agree to be bound by state and federal
laws and regulations governing equal opportunity and non-discrimination and
immigration.

1.4.   Appropriation of Funds. The Parties recognize that performance by
ARIZONA may depend upon appropriation of funds by the State Legislature of
ARIZONA. If the Legislature fails to appropriate the necessary funds, or if
ARIZONA’S appropriation is reduced during the fiscal year, ARIZONA may cancel
this Agreement without further duty or obligation. ARIZONA will notify
LICENSEE as soon as reasonably possible after it knows of the loss of funds.

1.5.   Conflict of Interest. This Agreement is subject to the provisions of
A.R.S. 38-511 and other conflict of interest regulations. Within three years
of the EFFECTIVE DATE, ARIZONA may cancel this Agreement if any person
significantly involved in initiating, negotiating, drafting, securing, or
creating this Agreement for or on behalf of ARIZONA becomes an employee or
consultant in any capacity of LICENSEE with respect to the subject matter of
this Agreement.

 */


#include "extern.h"

/*
 * Grace - Graphics for Exploratory Data Analysis
 * 
 * Home page: http://plasma-gate.weizmann.ac.il/Grace/
 * 
 * Copyright (c) 1991-95 Paul J Turner, Portland, OR
 * Copyright (c) 1996-98 GRACE Development Team
 * 
 * Maintained by Evgeny Stambulchik <fnevgeny@plasma-gate.weizmann.ac.il>
 * 
 * 
 *                           All Rights Reserved
 * 
 *    This program is free software; you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation; either version 2 of the License, or
 *    (at your option) any later version.
 * 
 *    This program is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 * 
 *    You should have received a copy of the GNU General Public License
 *    along with this program; if not, write to the Free Software
 *    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/*
 *
 * curve fitting, and other numerical routines used in compose.
 *
 * Contents:
 *
 * void   stasum()            - compute mean and variance
 * double leasev()            - evaluate least squares polynomial
 * int    fitcurve()          - compute coefficients for a polynomial fit of degree >1
 * void   runavg()            - compute a running average
 * void   runstddev()         - compute a running standard deviation
 * void   runmedian()         - compute a running median
 * void   runminmax()         - compute a running minimum or maximum
 * void   filterser()         - apply a digital filter
 * void   linearconv()        - convolve one set with another
 * int    crosscorr()         - cross/auto correlation
 * int    linear_regression() - linear regression
 * void   spline()            - compute a spline fit
 * double seval()             - evaluate the spline computed in spline()
 */

#include <../calc/config.h>
#include <math.h>

#include "defines.h"
#include "utils.h"
#include "channels.h"

int dofitcurve(int cnt, double *xd, double *yd, int nd, double *c);

/*
 * free and check for NULL pointer
 */
void xfree(void *ptr)
{
    if (ptr != NULL) {
        free(ptr);
    }
}


/*
 * sort two arrays
 */
void sort_xy(double *tmp1, double *tmp2, int up, int sorton, int stype)
{

    int d, i, j;
    int lo = 0;
    double t1, t2;

    if (sorton == 1) {
        double *ttmp;

        ttmp = tmp1;
        tmp1 = tmp2;
        tmp2 = ttmp;
    }
    up--;

    for (d = up - lo + 1; d > 1;) {
        if (d < 5)
            d = 1;
        else
            d = (5 * d - 1) / 11;
        for (i = up - d; i >= lo; i--) {
            t1 = tmp1[i];
            t2 = tmp2[i];
            if (!stype) {
                for (j = i + d; j <= up && (t1 > tmp1[j]); j += d) {
                    tmp1[j - d] = tmp1[j];
                    tmp2[j - d] = tmp2[j];
                }
                tmp1[j - d] = t1;
                tmp2[j - d] = t2;
            } else {
                for (j = i + d; j <= up && (t1 < tmp1[j]); j += d) {
                    tmp1[j - d] = tmp1[j];
                    tmp2[j - d] = tmp2[j];
                }
                tmp1[j - d] = t1;
                tmp2[j - d] = t2;
            }
        }
    }
}


/*
        compute mean and standard dev
*/
void stasum(double *x, int n, double *xbar, double *sd)
{
    int i;

    *xbar = 0;
    *sd = 0;
    
    if (x == NULL) {
        return;
    }
    
    if (n < 1) {
        return;
    }
    
    for (i = 0; i < n; i++) {
        *xbar = (*xbar) + x[i];
    }
    *xbar = (*xbar)/n;
    
    for (i = 0; i < n; i++) {
        *sd = (*sd) + (x[i] - *xbar) * (x[i] - *xbar);
    }
    *sd = sqrt(*sd/n);
}


/*
        evaluate least squares polynomial
*/
double leasev(double *c, int degree, double x)
{
  double temp;
  int i;

    /*
     * Commented out - refer to "Numerical Methods in C" for the
     * reason, section 5.2! The revolution has begun.
    temp = 0.0;
    for (i = 0; i <= degree; i++) {
        if ((i == 0) && (x == 0.0)) {
            temp = temp + c[i];
        } else {
            temp = temp + c[i] * pow(x, (double) (i));
        }
    } */

  temp = c[degree];

  for(i=degree-1;i>=0;i--)
  {
    temp = c[i] + temp*x;
  }

  return (temp);
}


/* polynomial curve fitting */
int fitcurve(double *x, double *y, int n, int ideg, double *coeff)
/*
 * x - x values
 * y - y values
 * n - number of points
 * ideg - degree of fit
 * coeff - coefficients of fit
 */
{
  int i, ifail;

  ifail = 1;

  if(ideg > 1)
  {
    dofitcurve(n, x, y, ideg, coeff);
  }
  else
  {
    ifail = linear_regression(n, x, y, coeff);

    if(ifail == 1)
    {
      st_fail("Linear_regression entered with N < 2");
      return ifail;
    } 
    else 
    if(ifail == 2)
    {
      st_fail("Linear_regression - all values of x or y are the same");
      return ifail;
    }
  }
        /* check coefficients */
  for(i=0;i<=ideg;i++)
  {
    if(!gsl_finite(coeff[i]))
    {
      st_fail("Linear_regression - all values of x or y are the same");
      ifail = 3;
      return ifail;
    }
  }

  ifail = 0;

  return(ifail);
}


/* compute a running average */
void runavg(double *x, double *y, double *ax, double *ay, int n, int ilen)
{
    int i;
    double sumy = 0.0;
    double sumx = 0.0;

    for (i = 0; i < ilen; i++) {
        sumx = sumx + x[i];
        sumy = sumy + y[i];
    }
    ax[0] = sumx / ilen;
    ay[0] = sumy / ilen;
    for (i = 1; i < (n - ilen + 1); i++) {
        sumx = x[i + ilen - 1] - x[i - 1] + sumx;
        ax[i] = sumx / ilen;
        sumy = y[i + ilen - 1] - y[i - 1] + sumy;
        ay[i] = sumy / ilen;
    }
}

/*
        compute a running standard deviation
*/
void runstddev(double *x, double *y, double *ax, double *ay, int n, int ilen)
{
    int i;
    double ybar, ysd;
    double sumx = 0.0;

    for (i = 0; i < ilen; i++) {
        sumx = sumx + x[i];
    }
    ax[0] = sumx / ilen;
    stasum(y, ilen, &ybar, &ysd);
    ay[0] = ysd;
    for (i = 1; i < (n - ilen + 1); i++) {
        stasum(y + i, ilen, &ybar, &ysd);
        sumx = x[i + ilen - 1] - x[i - 1] + sumx;
        ax[i] = sumx / ilen;
        ay[i] = ysd;
    }
}

/*
        compute a running median
*/
void runmedian(double *x, double *y, double *ax, double *ay, int n, int ilen)
{
    int i, j, nlen = n - ilen + 1;
    double *tmpx, *tmpy;

    tmpx = (double *) calloc(ilen, sizeof(double));
    if (tmpx == NULL) {
        st_fail("Can't calloc tmpx in runmedian");
        return;
    }
    tmpy = (double *) calloc(ilen, sizeof(double));
    if (tmpy == NULL) {
        st_fail("Can't calloc tmpy in runmedian");
        cxfree(tmpx);
        return;
    }
    for (i = 0; i < nlen; i++) {
        for (j = 0; j < ilen; j++) {
            tmpx[j] = x[j + i];
            tmpy[j] = y[j + i];
        }
        sort_xy(tmpx, tmpy, ilen, 1, 0);

        if (ilen % 2) {
            ax[i] = x[i + (ilen / 2)];
            ay[i] = tmpy[ilen / 2];
        } else {
            ax[i] = (x[i + ilen / 2] + x[i + (ilen - 1) / 2]) * 0.5;
            ay[i] = (tmpy[ilen / 2] + tmpy[(ilen - 1) / 2]) * 0.5;
        }
    }
    cxfree(tmpx);
    cxfree(tmpy);
}

/*
        compute a running minimum or maximum
*/
void runminmax(double *x, double *y, double *ax, double *ay, int n, int ilen, int type)
{
    int i, j;
    double min, max;
    double sumx = 0.0;

    min = max = y[0];
    for (i = 0; i < ilen; i++) {
        sumx = sumx + x[i];
        if (min > y[i])
            min = y[i];
        if (max < y[i])
            max = y[i];
    }
    ax[0] = sumx / ilen;
    if (type == 0) {
        ay[0] = min;
    } else if (type == 1) {
        ay[0] = max;
    } else {
        st_fail("Unknown type in runminmax, setting type = min");
        type = 0;
    }
    for (i = 1; i < (n - ilen + 1); i++) {
        sumx = x[i + ilen - 1] - x[i - 1] + sumx;
        ax[i] = sumx / ilen;
        min = y[i];
        max = y[i];
        for (j = 0; j < ilen; j++) {
            if (min > y[i + j])
                min = y[i + j];
            if (max < y[i + j])
                max = y[i + j];
        }
        if (type == 0) {
            ay[i] = min;
        } else if (type == 1) {
            ay[i] = max;
        }
    }
}

/*
        Apply a digital filter of length len to a set in x, y,
        of length n with the results going to resx, resy.
        the length of the result is set by the caller
*/
void filterser(int n, double *x, double *y, double *resx, double *resy, double *h, int len)
{
    int i, j, outlen, eo, ld2;
    double sum;

    outlen = n - len + 1;
    eo = len % 2;
    ld2 = len / 2;
    for (i = 0; i < outlen; i++) {
        sum = 0.0;
        for (j = 0; j < len; j++) {
            sum = sum + h[j] * y[j + i];
        }
        resy[i] = sum;
        if (eo)
            resx[i] = x[i + ld2];
        else
            resx[i] = (x[i + ld2] + x[i + ld2 - 1]) / 2.0;
    }
}

/*
        linear convolution of set x (length n) with h (length m) and
        result to y. the length of y is set by the caller
*/
void linearconv(double *x, double *h, double *y, int n, int m)
{
    int i, j, itmp;

    for (i = 0; i < n + m - 1; i++) {
        for (j = 0; j < m; j++) {
            itmp = i - j;
            if ((itmp >= 0) && (itmp < n)) {
                y[i] = y[i] + h[j] * x[itmp];
            }
        }
    }
}

/*
 * cross correlation
 */
int crosscorr(double *x, double *y, int n, int maxlag, double *xcov, double *xcor)
{
    double xbar, xsd;
    double ybar, ysd;
    int i, j;

    if (maxlag + 2 > n)
        return 1;
    for (i = 0; i <= maxlag; i++) {
        stasum(&x[i], n - i, &xbar, &xsd);
        if (xsd == 0.0)
            return 2;
        stasum(y, n - i, &ybar, &ysd);
        if (ysd == 0.0)
            return 3;
        xcov[i] = 0.0;
        for (j = 0; j < n - i; j++) {
            xcov[i] += (y[j] - ybar) * (x[j + i] - xbar);
        }
        xcov[i] /= n - i;
        xcor[i] = xcov[i] / (xsd * ysd);
    }
    return 0;
}

/* References,

   _Aplied Linear Regression_, Weisberg
   _Elements of Statistical Computing_, Thisted

   Fits y = coef*x + intercept + e

   uses a 2 pass method for means and variances
*/
int linear_regression(int n, double *x, double *y, double *coeff)
{
  double xbar, ybar;            /* sample means */
  double sdx, sdy;              /* sample standard deviations */
  double sxy, rxy;              /* sample covariance and sample correlation */
  double SXX, SYY, SXY;         /* sums of squares */
  double RSS;                   /* residual sum of squares */
  double rms;                   /* residual mean square */
//  double sereg;                       /* standard error of regression */
  double seslope, seintercept;
  double slope, intercept;      /* */
  double SSreg, F /*, R2*/;
  int i;

  if(n < 2)
  {
    return(1);
  } 

  xbar = ybar = 0.0;
  SXX = SYY = SXY = 0.0;

  for(i=0;i<n;i++)
  {
    xbar = xbar + x[i];
    ybar = ybar + y[i];
  }

  xbar = xbar / n;
  ybar = ybar / n;

  for(i=0;i<n;i++)
  {
    SXX = SXX + (x[i] - xbar) * (x[i] - xbar);
    SYY = SYY + (y[i] - ybar) * (y[i] - ybar);
    SXY = SXY + (x[i] - xbar) * (y[i] - ybar);
  }

  sdx = sqrt(SXX / (n - 1));
  sdy = sqrt(SYY / (n - 1));

  if(sdx == 0.0)
  {
    return(2);
  }

  if(sdy == 0.0)
  {
    return(2);
  }

  sxy       = SXY / (n - 1);
  rxy       = sxy / (sdx * sdy);
  slope     = SXY / SXX;
  intercept = ybar - slope * xbar;
  RSS       = SYY - slope * SXY;

  if(isSet(PF_V_FLAG))
  {
    st_report("Number of observations\t\t =    %d", n);
    st_report("Mean of independent variable\t = %12.7f", xbar);
    st_report("Mean of dependent variable\t\t = %12.7f", ybar);
    st_report("Standard dev. of ind. variable\t = %12.7f", sdx);
    st_report("Standard dev. of dep. variable\t = %12.7f", sdy);
    st_report("Correlation coefficient\t\t = %12.7f", rxy);
    st_report("Regression coefficient (SLOPE)\t = %12.7f", slope);
  }

  if(n == 2)
  {
    coeff[1] = (y[1] - y[0])/(x[1] - x[0]);
    coeff[0] = y[0] - coeff[1]*x[0];

    if(isSet(PF_V_FLAG))
    {
      st_report("Regression constant (INTERCEPT)\t = %12.7f\n", intercept);
    }

    return(0);
  } 
    
  rms         = RSS / (n - 2);
//  sereg       = sqrt(rms);
  seintercept = sqrt(rms * (1.0 / n + xbar * xbar / SXX));
  seslope     = sqrt(rms / SXX);
  SSreg       = SYY - RSS;
  F           = SSreg / rms;
//  R2          = SSreg / SYY;

  if(isSet(PF_V_FLAG))
  {
    st_report("Standard error of coefficient\t = %12.7f", seslope);
    st_report("t - value for coefficient\t\t = %12.7f", slope / seslope);
    st_report("Regression constant (INTERCEPT)\t = %12.7f", intercept);
    st_report("Standard error of constant\t\t = %12.7f", seintercept);
    st_report("t - value for constant\t\t = %12.7f", intercept / seintercept);
    st_report("Analysis of variance");

    st_report("Source\t d.f\t Sum of squares\t Mean Square\t F");
    st_report("Regression\t    1\t%12.7f\t%12.7f\t%12.7f", SSreg, SSreg, F);
    st_report("Residual\t%5d\t%12.7f\t%12.7f", n - 2, RSS, rms);
    st_report("Total\t%5d\t%12.7f", n - 1, SYY);
  }

  for(i=0;i<n;i++)
  {
    coeff[0] = intercept;
    coeff[1] = slope;
    coeff[2] = rxy;
  }

  return(0);
}

/*
        a literal translation of the spline routine in
        Forsyth, Malcolm, and Moler
*/
void spline(int n, double *x, double *y, double *b, double *c, double *d)
{
/*
c
c  the coefficients b(i), c(i), and d(i), i=1,2,...,n are computed
c  for a cubic interpolating spline
c
c    s(x) = y(i) + b(i)*(x-x(i)) + c(i)*(x-x(i))**2 + d(i)*(x-x(i))**3
c
c    for  x(i) .le. x .le. x(i+1)
c
c  input..
c
c    n = the number of data points or knots (n.ge.2)
c    x = the abscissas of the knots in strictly increasing order
c    y = the ordinates of the knots
c
c  output..
c
c    b, c, d  = arrays of spline coefficients as defined above.
c
c  using  p  to denote differentiation,
c
c    y(i) = s(x(i))
c    b(i) = sp(x(i))
c    c(i) = spp(x(i))/2
c    d(i) = sppp(x(i))/6  (derivative from the right)
c
c  the accompanying function subprogram  seval  can be used
c  to evaluate the spline.
c
c
*/

    int nm1, ib, i;
    double t;

/*
Gack!
*/
    x--;
    y--;
    b--;
    c--;
    d--;

/*
Fortran 66
*/
    nm1 = n - 1;
    if (n < 2)
        return;
    if (n < 3)
        goto l50;
/*
c
c  set up tridiagonal system
c
c  b = diagonal, d = offdiagonal, c = right hand side.
c
*/
    d[1] = x[2] - x[1];
    c[2] = (y[2] - y[1]) / d[1];
    for (i = 2; i <= nm1; i++) {
        d[i] = x[i + 1] - x[i];
        b[i] = 2.0 * (d[i - 1] + d[i]);
        c[i + 1] = (y[i + 1] - y[i]) / d[i];
        c[i] = c[i + 1] - c[i];
    }
/*
c
c  end conditions.  third derivatives at  x(1)  and  x(n)
c  obtained from divided differences
c
*/
    b[1] = -d[1];
    b[n] = -d[n - 1];
    c[1] = 0.0;
    c[n] = 0.0;
    if (n == 3)
        goto l15;
    c[1] = c[3] / (x[4] - x[2]) - c[2] / (x[3] - x[1]);
    c[n] = c[n - 1] / (x[n] - x[n - 2]) - c[n - 2] / (x[n - 1] - x[n - 3]);
    c[1] = c[1] * d[1] * d[1] / (x[4] - x[1]);
    c[n] = -c[n] * d[n - 1] * d[n - 1] / (x[n] - x[n - 3]);
/*
c
c  forward elimination
c
*/
l15:;
    for (i = 2; i <= n; i++) {
        t = d[i - 1] / b[i - 1];
        b[i] = b[i] - t * d[i - 1];
        c[i] = c[i] - t * c[i - 1];
    }
/*
c
c  back substitution
c
*/
    c[n] = c[n] / b[n];
    for (ib = 1; ib <= nm1; ib++) {
        i = n - ib;
        c[i] = (c[i] - d[i] * c[i + 1]) / b[i];
    }
/*
c
c  c(i) is now the sigma(i) of the text
c
c  compute polynomial coefficients
c
*/
    b[n] = (y[n] - y[nm1]) / d[nm1] + d[nm1] * (c[nm1] + 2.0 * c[n]);
    for (i = 1; i <= nm1; i++) {
        b[i] = (y[i + 1] - y[i]) / d[i] - d[i] * (c[i + 1] + 2.0 * c[i]);
        d[i] = (c[i + 1] - c[i]) / d[i];
        c[i] = 3.0 * c[i];
    }
    c[n] = 3.0 * c[n];
    d[n] = d[n - 1];
    return;

l50:;
    b[1] = (y[2] - y[1]) / (x[2] - x[1]);
    c[1] = 0.0;
    d[1] = 0.0;
    b[2] = b[1];
    c[2] = 0.0;
    d[2] = 0.0;
    return;
}

/***************************************************************************
 * aspline - modified version of David Frey's spline.c                     *
 *                                                                         *    
 * aspline does an Akima spline interpolation.                             *
 ***************************************************************************/

void aspline(int n, double *x, double *y, double *b, double *c, double *d)
{
  int i;
        
  double num, den;
  double m_m1, m_m2, m_p1, m_p2;
  double x_m1, x_m2, x_p1, x_p2;
  double y_m1, y_m2, y_p1, y_p2;

#define dx(i) (x[i+1]-x[i])
#define dy(i) (y[i+1]-y[i])
#define  m(i) (dy(i)/dx(i))

  if (n > 0)                 /* we have data to process */
  {

      /*
       * calculate the coefficients of the spline 
       * (the Akima interpolation itself)                      
       */

      /* b) interpolate the missing points: */

      x_m1 = x[0] + x[1] - x[2]; 
      y_m1 = (x[0]-x_m1) * (m(1) - 2 * m(0)) + y[0];

      m_m1 = (y[0]-y_m1)/(x[0]-x_m1);
       
      x_m2 = 2 * x[0] - x[2];
      y_m2 = (x_m1-x_m2) * (m(0) - 2 * m_m1) + y_m1;
       
      m_m2 = (y_m1-y_m2)/(x_m1-x_m2);

      x_p1 = x[n-1] + x[n-2] - x[n-3];
      y_p1 = (2 * m(n-2) - m(n-3)) * (x_p1 - x[n-1]) + y[n-1];

      m_p1 = (y_p1-y[n-1])/(x_p1-x[n-1]);
      
      x_p2 = 2 * x[n-1] - x[n-3];
      y_p2 = (2 * m_p1 - m(n-2)) * (x_p2 - x_p1) + y_p1;
      
      m_p2 = (y_p2-y_p1)/(x_p2-x_p1);
           
      /* i = 0 */
      num=fabs(m(1) - m(0)) * m_m1 + fabs(m_m1 - m_m2) * m(0);
      den=fabs(m(1) - m(0)) + fabs(m_m1 - m_m2);
        
      if (den != 0.0) b[0]=num / den;
      else            b[0]=0.0;
                
      /* i = 1 */
      num=fabs(m(2) - m(1)) * m(0) + fabs(m(0) - m_m1) * m(1);
      den=fabs(m(2) - m(1)) + fabs(m(0) - m_m1);

      if (den != 0.0) b[1]=num / den;
      else            b[1]=0.0;
                        
      for (i=2; i < n-2; i++)
      {

        num=fabs(m(i+1) - m(i)) * m(i-1) + fabs(m(i-1) - m(i-2)) * m(i);
        den=fabs(m(i+1) - m(i)) + fabs(m(i-1) - m(i-2));

        if (den != 0.0) b[i]=num / den;
        else            b[i]=0.0;
      }

      /* i = n - 2 */
      num=fabs(m_p1 - m(n-2)) * m(n-3) + fabs(m(n-3) - m(n-4)) * m(n-2);
      den=fabs(m_p1 - m(n-2)) + fabs(m(n-3) - m(n-4));

      if (den != 0.0) b[n-2]=num / den;
      else            b[n-2]=0.0;
 
      /* i = n - 1 */
      num=fabs(m_p2 - m_p1) * m(n-2) + fabs(m(n-2) - m(n-3)) * m_p1;
      den=fabs(m_p2 - m_p1) + fabs(m(n-2) - m(n-3));

      if (den != 0.0) b[n-1]=num / den;
      else            b[n-1]=0.0;
 
      for (i=0; i < n-1; i++)
      {
           double dxv = dx(i);
           c[i]=(3 * m(i) - 2 * b[i] - b[i+1]) / dxv;
           d[i]=(b[i] + b[i+1] - 2 * m(i)) / (dxv * dxv);
      }
  }
#undef dx
#undef dy
#undef  m
}

double seval(int n, double u, double *x, double *y, double *b, double *c, double *d)
{
/*
c
c  this subroutine evaluates the cubic spline function
c
c    seval = y(i) + b(i)*(u-x(i)) + c(i)*(u-x(i))**2 + d(i)*(u-x(i))**3
c
c    where  x(i) .lt. u .lt. x(i+1), using horner's rule
c
c  if  u .lt. x(1) then  i = 1  is used.
c  if  u .ge. x(n) then  i = n  is used.
c
c  input..
c
c    n = the number of data points
c    u = the abscissa at which the spline is to be evaluated
c    x,y = the arrays of data abscissas and ordinates
c    b,c,d = arrays of spline coefficients computed by spline
c
c  if  u  is not in the same interval as the previous call, then a
c  binary search is performed to determine the proper interval.
c
*/
    int j, k;
    double dx;
    int i;

/*
c
c  binary search
c
*/
    if (u < x[0]) {
        i = 0;
    } else if (u >= x[n - 1]) {
        i = n - 1;
    } else {
        i = 0;
        j = n;
l20:    ;
        k = (i + j) / 2;
        if (u < x[k])
            j = k;
        if (u >= x[k])
            i = k;
        if (j > i + 1)
            goto l20;
    }
/*
c
c  evaluate spline
c
*/
    dx = u - x[i];
    return (y[i] + dx * (b[i] + dx * (c[i] + dx * d[i])));
}