home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 7 / 07.iso / c / c129 / 1.ddi / CURVEFIT.C < prev    next >
Encoding:
C/C++ Source or Header  |  1989-08-17  |  3.4 KB  |  164 lines
  1.  
  2. # define nearzero 1.0E-20
  3. # include <math.h>
  4. # include <stdlib.h>
  5.  
  6. void PolyCurveFit(float *indvar,float *depvar,int numobs,int order,
  7.                   float *coef,float *yest,float *resid,float *see,
  8.                   float *coefsig, float *rsq,float *r,char *cferror)
  9. {
  10.    int i;
  11.    int j;
  12.    int numcol;
  13.    float *regmat;
  14.    float seexx;
  15.    float rsqxx;
  16.    float rxx;
  17.    float cferrorxx;
  18.  
  19.  
  20.  
  21.    regmat = (float *) calloc(numobs*(order),4);
  22.  
  23.    if ( order + 1 < numobs ) {
  24.       for ( i = 0; i <= numobs-1; ++i ) {
  25.          regmat[i*order] = indvar[i];
  26.          for ( j = 1; j <= order-1; ++j ) {
  27.             regmat[i*order+j] = indvar[i] * regmat[i*order+j - 1];
  28.          }
  29.       }
  30.       MultipleReg(regmat,depvar,order,numobs,
  31.       coef,yest,resid,see,coefsig,rsq,r,cferror);
  32.    }
  33.    else {
  34.       (*cferror) = 1;
  35.    }
  36.    free(regmat);
  37. }
  38.  
  39. void tridiag(float *mat,int n,float *vec)
  40. {
  41.    int i;
  42.    float pivot;
  43.    float mult;
  44.    float ci;
  45.    float bi;
  46.    char c;
  47.  
  48.    for ( i = 1; i <= n - 1; ++i ) {
  49.       pivot = mat[i*4+0];
  50.       ci = mat[i*4+1];
  51.       bi = mat[i*4+2];
  52.       mult = mat[(i+1)*4+3] / pivot;
  53.     
  54.       if ( fabs(mult) > nearzero ) {
  55.          mat[(i+1)*4+3] = mult;
  56.          mat[(i+1)*4+0] = mat[(i+1)*4+0] - mult * ci;
  57.          mat[(i+1)*4+2] = mat[(i+1)*4+2] - mult * bi;
  58.  
  59.        }
  60.       else {
  61.          mat[(i+1)*4+3] = 0.0;
  62.       }
  63.    }
  64.     vec[n] = mat[n*4+2] / mat[n*4+0];
  65.    for ( i = n - 1; i >= 1; --i ) {
  66.       vec[i] = (mat[i*4+2] - mat[i*4+1] * vec[i+1]) / mat[i*4+0];
  67.  } }
  68.  
  69. void CubicSplines(float *x,float *y,int n,float *s)
  70. {
  71.    int i;
  72.    float *eqn;
  73.    float bb;
  74.    float cc;
  75.    float hh;
  76.    float ff;
  77.    float fi;
  78.    float hi;
  79.    float nextf;
  80.    float *h;
  81.    float *bvec;
  82.    char c;
  83.  
  84.    eqn = (float *) calloc((n+1)*4,4);
  85.  
  86.    h = (float *) calloc(n+1,4);
  87.    bvec = (float *) calloc(n+1,4);
  88.   
  89.    for ( i = 0; i <= n - 1; ++i ) {
  90.       h[i] = x[i+1] - x[i];
  91.    }
  92.    hh = h[0];
  93.    ff = y[0];
  94.    nextf = y[1];
  95.    
  96.    for ( i = 1; i <= n - 1; ++i ) {
  97.       fi = nextf;
  98.       nextf = y[i+1];
  99.       hi = h[i];
  100.       eqn[i*4+2] = 3.0 * ((nextf - fi) / hi - (fi - ff) / hh);
  101.       eqn[i*4+0] = 2.0 * (hh + hi);
  102.       eqn[i*4+3] = hh;
  103.       eqn[i*4+1] = hi;
  104.       hh = hi;
  105.       ff = fi;
  106.    }
  107.     tridiag(eqn,n - 1,bvec);
  108.     hh = h[0];
  109.    s[2] = bb;
  110.    s[0] = y[0];
  111.    cc = (y[1] - s[0]) / hh - bvec[1] * hh / 3.0;
  112.    s[1] = cc;
  113.    for ( i = 1; i <= n - 1; ++i ) {
  114.       s[i*4+2] = bvec[i];
  115.       s[i*4+0] = y[i];
  116.       cc = (s[i*4+2] + bb) * hh + cc;
  117.   
  118.        
  119.       s[i*4+1] = cc;
  120.       s[(i-1)*4+3] = (s[i*4+2] - bb) / (3.0 * hh);
  121.       bb = s[i*4+2];
  122.       hh = h[i];
  123.     }
  124.    s[(n-1)*4+3] = -bb / (3.0 * hh);
  125.    free(eqn); free(h); free(bvec);
  126. }
  127.  
  128. float PolyCalc(float XIn,float *CoefVector,int order)
  129.  
  130. {
  131.    float _fret;
  132.    int i;
  133.    float YOut;
  134.  
  135.    YOut = CoefVector[order];
  136.    for ( i = order - 1; i >= 0; --i ) {
  137.       YOut = YOut * XIn + CoefVector[i];
  138.    }
  139.    _fret = YOut;
  140.    return(_fret);
  141. }
  142.  
  143. void CalcSpline(float *xv,float *coef,int n,float x,float *y)
  144. {
  145.    int i;
  146.    int j;
  147.    int k;
  148.    float *CoefVector;
  149.    float dif;
  150.  
  151.    CoefVector = (float *) calloc(4,4);
  152.  
  153.    i = -1;
  154.    do {
  155.       i = i + 1;
  156.    } while ( ! (((x >= xv[i]) && (x < xv[i + 1])) | (i == n)) );
  157.    dif = x - xv[i];
  158.    for ( j = 0; j <= 3; ++j ) {
  159.       CoefVector[j] = coef[i*4+j];
  160.    }
  161.    (*y) = PolyCalc(dif,CoefVector,3);
  162.    free(CoefVector);
  163. }
  164.