home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l292 / 1.ddi / CURVEFIT.FOR < prev    next >
Encoding:
Text File  |  1990-02-20  |  3.8 KB  |  145 lines
  1.  
  2.       SUBROUTINE CalcSpline (xv, coef, n, x, y)
  3.       INCLUDE 'STDHDR.FOR'
  4.       REAL xv(0:maxv),coef(0:maxr,0:maxc), x, y
  5.       INTEGER n, i
  6.       REAL CoefVector(0:maxc), dif
  7.  
  8.       i = -1
  9.  
  10.       DO WHILE  ((x .LE. xv(i) .AND. x .LT. xv(i + 1)) .AND. i .NE. n)
  11.         i = i + 1
  12.       END DO
  13.       dif = x - xv(i)
  14.       DO j = 0, 3
  15.         CoefVector(j) = coef(i, j)
  16.       END DO
  17.       y = PolyCurveCalc(dif, CoefVector, 3)
  18.       END !SUBROUTINE CalcSpline
  19.  
  20.  
  21.  
  22.       FUNCTION PolyCurveCalc (XIn, CoefVector, order)
  23.       INCLUDE 'STDHDR.FOR'
  24.       REAL xin, CoefVector(0:maxc), YOut
  25.       INTEGER order
  26.  
  27.       YOut = CoefVector(order)
  28.       DO i = order - 1, 0,  -1
  29.         YOut = YOut * XIn + CoefVector(i)
  30.       END DO
  31.       PolyCurveCalc = YOut
  32.       END ! FUNCTION   PolyCurveCalc
  33.  
  34.  
  35.  
  36.  
  37.       SUBROUTINE CubicSplines (x, y, n, s)
  38.       INCLUDE 'STDHDR.FOR'
  39.       REAL x(0:maxr),y(0:maxr),s(0:maxr,0:maxc), nearzero
  40.       INTEGER n,iErr
  41.       PARAMETER ( nearzero = 1E-20)
  42.       REAL eqn[ALLOCATABLE](:,:), h[ALLOCATABLE](:),bvec[ALLOCATABLE](:)
  43.       INTEGER i
  44.       REAL bb, cc, hh,ff, fi, hi, nextf
  45.       ALLOCATE(eqn(0:maxr,0 :maxc), h(0:maxr),bvec(0:maxr),STAT=iErr)
  46.       CALL CheckMem(iErr)
  47.       DO i = 0, n - 1
  48.         h(i) = x(i + 1) - x(i)
  49.       END DO
  50.       hh = h(0)
  51.       ff = y(0)
  52.       nextf = y(1)
  53.       DO i = 1, n - 1
  54.         fi = nextf
  55.         nextf = y(i + 1)
  56.         hi = h(i)
  57.         eqn(i, 2) = 3.0 * ((nextf - fi) / hi - (fi - ff) / hh)
  58.         eqn(i, 0) = 2.0 * (hh + hi)
  59.         eqn(i, 3) = hh
  60.         eqn(i, 1) = hi
  61.         hh = hi
  62.         ff = fi
  63.       END DO
  64.       CALL tridiag(eqn, n - 1, bvec)
  65.       bb = 0.0
  66.       hh = h(0)
  67.       s(0, 2) = bb
  68.       s(0, 0) = y(0)
  69.       cc = (y(1) - s(0, 0)) / hh - bvec(1) * hh / 3.0
  70.       s(0, 1) = cc
  71.       DO i = 1, n - 1
  72.         s(i, 2) = bvec(i)
  73.         s(i, 0) = y(i)
  74.         cc = (s(i, 2) + bb) * hh + cc
  75.         s(i, 1) = cc
  76.         s(i - 1, 3) = (s(i, 2) - bb) / (3.0 * hh)
  77.         bb = s(i, 2)
  78.         hh = h(i)
  79.       END DO
  80.       s(n - 1, 3) = -bb / (3.0 * hh)
  81.       DEALLOCATE(eqn, h,bvec,STAT = iErr)
  82.       CALL CheckDealloc(iErr)
  83.       END !SUBROUTINE Cubic Splines
  84.  
  85.       SUBROUTINE tridiag (mat, n, vec)
  86.       INCLUDE 'STDHDR.FOR'
  87.       REAL  mat(0:maxr, 0:maxc), vec(0:maxr), nearzero
  88.       INTEGER n
  89.       PARAMETER ( nearzero = 1E-20)
  90.       REAL pivot,mult, ci,bi
  91.  
  92.       DO i = 1, n - 1
  93.         pivot = mat(i, 0)
  94.         ci = mat(i, 1)
  95.         bi = mat(i, 2)
  96.         mult = mat(i + 1, 3) / pivot
  97.         IF (ABS(mult) .GT. nearzero) THEN
  98.           mat(i + 1, 3) = mult
  99.           mat(i + 1, 0) = mat(i + 1, 0) - mult * ci
  100.           mat(i + 1, 2) = mat(i + 1, 2) - mult * bi
  101.         ELSE
  102.           mat(i + 1, 3) = 0
  103.         END IF
  104.       END DO
  105.  
  106.       vec(n) = mat(n, 2) / mat(n, 0)
  107.       DO i = n - 1, 1,  -1
  108.         vec(i) = (mat(i, 2) - mat(i, 1) * vec(i + 1)) / mat(i, 0)
  109.       END DO
  110.       END !SUBROUTINE Triag
  111.  
  112.  
  113.  
  114.  
  115.       SUBROUTINE PolyCurveFit (indvar, depvar, numobs, order, coef,
  116.      +       yest, resid, see, coefsig, rsq, r, cferror)
  117.       INCLUDE 'STDHDR.FOR'
  118.  
  119.       REAL indvar(0:maxr), depvar(0:maxr), coef(0:maxc), yest(0:maxr)
  120.       REAL resid(0:maxr), coefsig(0:maxc)
  121.       REAL regmat[ALLOCATABLE](:,:)
  122.       INTEGER numobs, order,iErr
  123.       REAL r, rsq
  124.       LOGICAL cferror
  125.  
  126.       ALLOCATE(regmat(0: maxr, 0:maxc),STAT=iErr)
  127.  
  128.       IF (order + 1 .LT. numobs) THEN
  129.         DO i = 0, numobs - 1
  130.           regmat(i, 0) = indvar(i)
  131.           DO j = 1, order - 1
  132.             regmat(i, j) = indvar(i) * regmat(i, j - 1)
  133.           END DO
  134.         END DO
  135.  
  136.         CALL MultipleReg(regmat, depvar, order, numobs,
  137.      +        coef, yest, resid, see, coefsig, rsq, r, cferror)
  138.       ELSE
  139.         cferror = .TRUE.
  140.       END IF
  141.       DEALLOCATE(regmat,STAT=iErr)
  142.  
  143.       END !SUBROUTINE  PolyCurvefit
  144.  
  145.