home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / libcruft / villad / intrp.f < prev    next >
Text File  |  1996-09-28  |  2KB  |  94 lines
  1.       SUBROUTINE INTRP ( ND, NT, X, ROOT, DIF1, XINTP )
  2. C
  3.       INTEGER           ND, NT
  4.       DOUBLE PRECISION  ROOT(ND), DIF1(ND), XINTP(ND)
  5. C
  6. C***********************************************************************
  7. C
  8. C     LAGRANGE INTERPOLATION
  9. C
  10. C     VILLADSEN AND MICHELSEN, PAGES 132-133, 420
  11. C
  12. C     INPUT PARAMETERS:
  13. C
  14. C       NT     : THE TOTAL NUMBER OF INTERPOLATION POINTS FOR WHICH THE
  15. C                VALUE OF THE DEPENDENT VARIABLE Y IS KNOWN.  NOTE:
  16. C
  17. C                  NT = N + N0 + N1
  18. C
  19. C       X      : THE ABCISSA X WHERE Y(X) IS DESIRED
  20. C
  21. C       ROOT   : ONE DIMENSIONAL VECTOR CONTAINING ON EXIT THE
  22. C                N + N0 + N1 ZEROS OF THE NODE POLYNOMIAL USED IN THE
  23. C                INTERPOLATION ROUTINE
  24. C
  25. C       DIF1   : ONE DIMENSIONAL VECTOR CONTAINING THE FIRST DERIVATIVE
  26. C                OF THE NODE POLYNOMIAL AT THE ZEROS
  27. C
  28. C     OUTPUT PARAMETERS:
  29. C
  30. C       XINTP  : THE VECTOR OF INTERPOLATION WEIGHTS
  31. C
  32. C                Y(X) IS GIVEN BY:
  33. C
  34. C                            NT
  35. C                  Y(X)  =  SUM  XINTRP(I) * Y(I)
  36. C                           I=1
  37. C
  38. C     COMMON BLOCKS:      NONE
  39. C
  40. C     REQUIRED ROUTINES:  VILERR
  41. C
  42. C***********************************************************************
  43. C
  44.       INTEGER           I,IER
  45.       DOUBLE PRECISION  POL,Y,X
  46.       DOUBLE PRECISION  ZERO,ONE
  47.       LOGICAL           LSTOP
  48. C
  49.       PARAMETER ( ZERO = 0.0D+00, ONE = 1.0D+00 )
  50. C
  51. C -- ERROR CHECKING
  52. C
  53.       IF (ND .LT. NT) THEN
  54.         IER   = 3
  55.         LSTOP = .TRUE.
  56.         CALL VILERR(IER,LSTOP)
  57.       ELSE
  58.       END IF
  59. C
  60.       IF (NT .LT. 1) THEN
  61.         IER   = 7
  62.         LSTOP = .TRUE.
  63.         CALL VILERR(IER,LSTOP)
  64.       ELSE
  65.       END IF
  66. C
  67. C -- EVALUATE LAGRANGIAN INTERPOLATION COEFFICIENTS
  68. C
  69.       POL = ONE
  70. C
  71.       DO 5 I = 1,NT
  72. C
  73.         Y        = X - ROOT(I)
  74.         XINTP(I) = ZERO
  75. C
  76.         IF (Y .EQ. ZERO) THEN
  77.           XINTP(I) = ONE
  78.         ELSE
  79.         END IF
  80. C
  81.         POL = POL*Y
  82. C
  83.     5 CONTINUE
  84. C
  85.       IF (POL .NE. ZERO) THEN
  86.         DO 6 I = 1,NT
  87.           XINTP(I) = POL/DIF1(I)/(X - ROOT(I))
  88.     6   CONTINUE
  89.       ELSE
  90.       END IF
  91. C
  92.       RETURN
  93.       END
  94.