
      SUBROUTINE CalcSpline (xv, coef, n, x, y)
      INCLUDE 'STDHDR.FOR'
      REAL xv(0:maxv),coef(0:maxr,0:maxc), x, y
      INTEGER n, i
      REAL CoefVector(0:maxc), dif

      i = -1

      DO WHILE  ((x .LE. xv(i) .AND. x .LT. xv(i + 1)) .AND. i .NE. n)
        i = i + 1
      END DO
      dif = x - xv(i)
      DO j = 0, 3
        CoefVector(j) = coef(i, j)
      END DO
      y = PolyCurveCalc(dif, CoefVector, 3)
      END !SUBROUTINE CalcSpline



      FUNCTION PolyCurveCalc (XIn, CoefVector, order)
      INCLUDE 'STDHDR.FOR'
      REAL xin, CoefVector(0:maxc), YOut
      INTEGER order

      YOut = CoefVector(order)
      DO i = order - 1, 0,  -1
        YOut = YOut * XIn + CoefVector(i)
      END DO
      PolyCurveCalc = YOut
      END ! FUNCTION   PolyCurveCalc




      SUBROUTINE CubicSplines (x, y, n, s)
      INCLUDE 'STDHDR.FOR'
      REAL x(0:maxr),y(0:maxr),s(0:maxr,0:maxc), nearzero
      INTEGER n,iErr
      PARAMETER ( nearzero = 1E-20)
      REAL eqn[ALLOCATABLE](:,:), h[ALLOCATABLE](:),bvec[ALLOCATABLE](:)
      INTEGER i
      REAL bb, cc, hh,ff, fi, hi, nextf
      ALLOCATE(eqn(0:maxr,0 :maxc), h(0:maxr),bvec(0:maxr),STAT=iErr)
      CALL CheckMem(iErr)
      DO i = 0, n - 1
        h(i) = x(i + 1) - x(i)
      END DO
      hh = h(0)
      ff = y(0)
      nextf = y(1)
      DO i = 1, n - 1
        fi = nextf
        nextf = y(i + 1)
        hi = h(i)
        eqn(i, 2) = 3.0 * ((nextf - fi) / hi - (fi - ff) / hh)
        eqn(i, 0) = 2.0 * (hh + hi)
        eqn(i, 3) = hh
        eqn(i, 1) = hi
        hh = hi
        ff = fi
      END DO
      CALL tridiag(eqn, n - 1, bvec)
      bb = 0.0
      hh = h(0)
      s(0, 2) = bb
      s(0, 0) = y(0)
      cc = (y(1) - s(0, 0)) / hh - bvec(1) * hh / 3.0
      s(0, 1) = cc
      DO i = 1, n - 1
        s(i, 2) = bvec(i)
        s(i, 0) = y(i)
        cc = (s(i, 2) + bb) * hh + cc
        s(i, 1) = cc
        s(i - 1, 3) = (s(i, 2) - bb) / (3.0 * hh)
        bb = s(i, 2)
        hh = h(i)
      END DO
      s(n - 1, 3) = -bb / (3.0 * hh)
      DEALLOCATE(eqn, h,bvec,STAT = iErr)
      CALL CheckDealloc(iErr)
      END !SUBROUTINE Cubic Splines

      SUBROUTINE tridiag (mat, n, vec)
      INCLUDE 'STDHDR.FOR'
      REAL  mat(0:maxr, 0:maxc), vec(0:maxr), nearzero
      INTEGER n
      PARAMETER ( nearzero = 1E-20)
      REAL pivot,mult, ci,bi

      DO i = 1, n - 1
        pivot = mat(i, 0)
        ci = mat(i, 1)
        bi = mat(i, 2)
        mult = mat(i + 1, 3) / pivot
        IF (ABS(mult) .GT. nearzero) THEN
          mat(i + 1, 3) = mult
          mat(i + 1, 0) = mat(i + 1, 0) - mult * ci
          mat(i + 1, 2) = mat(i + 1, 2) - mult * bi
        ELSE
          mat(i + 1, 3) = 0
        END IF
      END DO

      vec(n) = mat(n, 2) / mat(n, 0)
      DO i = n - 1, 1,  -1
        vec(i) = (mat(i, 2) - mat(i, 1) * vec(i + 1)) / mat(i, 0)
      END DO
      END !SUBROUTINE Triag




      SUBROUTINE PolyCurveFit (indvar, depvar, numobs, order, coef,
     +       yest, resid, see, coefsig, rsq, r, cferror)
      INCLUDE 'STDHDR.FOR'

      REAL indvar(0:maxr), depvar(0:maxr), coef(0:maxc), yest(0:maxr)
      REAL resid(0:maxr), coefsig(0:maxc)
      REAL regmat[ALLOCATABLE](:,:)
      INTEGER numobs, order,iErr
      REAL r, rsq
      LOGICAL cferror

      ALLOCATE(regmat(0: maxr, 0:maxc),STAT=iErr)

      IF (order + 1 .LT. numobs) THEN
        DO i = 0, numobs - 1
          regmat(i, 0) = indvar(i)
          DO j = 1, order - 1
            regmat(i, j) = indvar(i) * regmat(i, j - 1)
          END DO
        END DO

        CALL MultipleReg(regmat, depvar, order, numobs,
     +        coef, yest, resid, see, coefsig, rsq, r, cferror)
      ELSE
        cferror = .TRUE.
      END IF
      DEALLOCATE(regmat,STAT=iErr)

      END !SUBROUTINE  PolyCurvefit

