home *** CD-ROM | disk | FTP | other *** search
/ Frostbyte's 1980s DOS Shareware Collection / floppyshareware.zip / floppyshareware / DOOG / PCSSP1.ZIP / ITRPAPSM.ZIP / AHI.FOR < prev    next >
Text File  |  1985-11-29  |  5KB  |  139 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE AHI
  5. C
  6. C        PURPOSE
  7. C           TO INTERPOLATE FUNCTION VALUE Y FOR A GIVEN ARGUMENT VALUE
  8. C           X USING A GIVEN TABLE (ARG,VAL) OF ARGUMENT, FUNCTION, AND
  9. C           DERIVATIVE VALUES.
  10. C
  11. C        USAGE
  12. C           CALL AHI (X,ARG,VAL,Y,NDIM,EPS,IER)
  13. C
  14. C        DESCRIPTION OF PARAMETERS
  15. C           X      - THE ARGUMENT VALUE SPECIFIED BY INPUT.
  16. C           ARG    - THE INPUT VECTOR (DIMENSION NDIM) OF ARGUMENT
  17. C                    VALUES OF THE TABLE (NOT DESTROYED).
  18. C           VAL    - THE INPUT VECTOR (DIMENSION 2*NDIM) OF FUNCTION
  19. C                    AND DERIVATIVE VALUES OF THE TABLE (DESTROYED).
  20. C                    FUNCTION AND DERIVATIVE VALUES MUST BE STORED IN
  21. C                    PAIRS, THAT MEANS BEGINNING WITH FUNCTION VALUE AT
  22. C                    POINT ARG(1) EVERY FUNCTION VALUE MUST BE FOLLOWED
  23. C                    BY THE VALUE OF DERIVATIVE AT THE SAME POINT.
  24. C           Y      - THE RESULTING INTERPOLATED FUNCTION VALUE.
  25. C           NDIM   - AN INPUT VALUE WHICH SPECIFIES THE NUMBER OF
  26. C                    POINTS IN TABLE (ARG,VAL).
  27. C           EPS    - AN INPUT CONSTANT WHICH IS USED AS UPPER BOUND
  28. C                    FOR THE ABSOLUTE ERROR.
  29. C           IER    - A RESULTING ERROR PARAMETER.
  30. C
  31. C        REMARKS
  32. C           (1) TABLE (ARG,VAL) SHOULD REPRESENT A SINGLE-VALUED
  33. C               FUNCTION AND SHOULD BE STORED IN SUCH A WAY, THAT THE
  34. C               DISTANCES ABS(ARG(I)-X) INCREASE WITH INCREASING
  35. C               SUBSCRIPT I. TO GENERATE THIS ORDER IN TABLE (ARG,VAL),
  36. C               SUBROUTINES ATSG, ATSM OR ATSE COULD BE USED IN A
  37. C               PREVIOUS STAGE.
  38. C           (2) NO ACTION BESIDES ERROR MESSAGE IN CASE NDIM LESS
  39. C               THAN 1.
  40. C           (3) INTERPOLATION IS TERMINATED EITHER IF THE DIFFERENCE
  41. C               BETWEEN TWO SUCCESSIVE INTERPOLATED VALUES IS
  42. C               ABSOLUTELY LESS THAN TOLERANCE EPS, OR IF THE ABSOLUTE
  43. C               VALUE OF THIS DIFFERENCE STOPS DIMINISHING, OR AFTER
  44. C               (2*NDIM-2) STEPS. FURTHER IT IS TERMINATED IF THE
  45. C               PROCEDURE DISCOVERS TWO ARGUMENT VALUES IN VECTOR ARG
  46. C               WHICH ARE IDENTICAL. DEPENDENT ON THESE FOUR CASES,
  47. C               ERROR PARAMETER IER IS CODED IN THE FOLLOWING FORM
  48. C                IER=0 - IT WAS POSSIBLE TO REACH THE REQUIRED
  49. C                        ACCURACY (NO ERROR).
  50. C                IER=1 - IT WAS IMPOSSIBLE TO REACH THE REQUIRED
  51. C                        ACCURACY BECAUSE OF ROUNDING ERRORS.
  52. C                IER=2 - IT WAS IMPOSSIBLE TO CHECK ACCURACY BECAUSE
  53. C                        NDIM IS LESS THAN 2, OR THE REQUIRED ACCURACY
  54. C                        COULD NOT BE REACHED BY MEANS OF THE GIVEN
  55. C                        TABLE. NDIM SHOULD BE INCREASED.
  56. C                IER=3 - THE PROCEDURE DISCOVERED TWO ARGUMENT VALUES
  57. C                        IN VECTOR ARG WHICH ARE IDENTICAL.
  58. C
  59. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  60. C           NONE
  61. C
  62. C        METHOD
  63. C           INTERPOLATION IS DONE BY MEANS OF AITKENS SCHEME OF
  64. C           HERMITE INTERPOLATION. ON RETURN Y CONTAINS AN INTERPOLATED
  65. C           FUNCTION VALUE AT POINT X, WHICH IS IN THE SENSE OF REMARK
  66. C           (3) OPTIMAL WITH RESPECT TO GIVEN TABLE. FOR REFERENCE, SEE
  67. C           F.B.HILDEBRAND, INTRODUCTION TO NUMERICAL ANALYSIS,
  68. C           MCGRAW-HILL, NEW YORK/TORONTO/LONDON, 1956, PP.314-317, AND
  69. C           GERSHINSKY/LEVINE, AITKEN-HERMITE INTERPOLATION,
  70. C           JACM, VOL.11, ISS.3 (1964), PP.352-356.
  71. C
  72. C     ..................................................................
  73. C
  74.       SUBROUTINE AHI(X,ARG,VAL,Y,NDIM,EPS,IER)
  75. C
  76. C
  77.       DIMENSION ARG(1),VAL(1)
  78.       IER=2
  79.       H2=X-ARG(1)
  80.       IF(NDIM-1)2,1,3
  81.     1 Y=VAL(1)+VAL(2)*H2
  82.     2 RETURN
  83. C
  84. C     VECTOR ARG HAS MORE THAN 1 ELEMENT.
  85. C     THE FIRST STEP PREPARES VECTOR VAL SUCH THAT AITKEN SCHEME CAN BE
  86. C     USED.
  87.     3 I=1
  88.       DO 5 J=2,NDIM
  89.       H1=H2
  90.       H2=X-ARG(J)
  91.       Y=VAL(I)
  92.       VAL(I)=Y+VAL(I+1)*H1
  93.       H=H1-H2
  94.       IF(H)4,13,4
  95.     4 VAL(I+1)=Y+(VAL(I+2)-Y)*H1/H
  96.     5 I=I+2
  97.       VAL(I)=VAL(I)+VAL(I+1)*H2
  98. C     END OF FIRST STEP
  99. C
  100. C     PREPARE AITKEN SCHEME
  101.       DELT2=0.
  102.       IEND=I-1
  103. C
  104. C     START AITKEN-LOOP
  105.       DO 9 I=1,IEND
  106.       DELT1=DELT2
  107.       Y=VAL(1)
  108.       M=(I+3)/2
  109.       H1=ARG(M)
  110.       DO 6 J=1,I
  111.       K=I+1-J
  112.       L=(K+1)/2
  113.       H=ARG(L)-H1
  114.       IF(H)6,14,6
  115.     6 VAL(K)=(VAL(K)*(X-H1)-VAL(K+1)*(X-ARG(L)))/H
  116.       DELT2=ABS(Y-VAL(1))
  117.       IF(DELT2-EPS)11,11,7
  118.     7 IF(I-5)9,8,8
  119.     8 IF(DELT2-DELT1)9,12,12
  120.     9 CONTINUE
  121. C     END OF AITKEN-LOOP
  122. C
  123.    10 Y=VAL(1)
  124.       RETURN
  125. C
  126. C     THERE IS SUFFICIENT ACCURACY WITHIN 2*NDIM-2 ITERATION STEPS
  127.    11 IER=0
  128.       GOTO 10
  129. C
  130. C     TEST VALUE DELT2 STARTS OSCILLATING
  131.    12 IER=1
  132.       RETURN
  133. C
  134. C     THERE ARE TWO IDENTICAL ARGUMENT VALUES IN VECTOR ARG
  135.    13 Y=VAL(1)
  136.    14 IER=3
  137.       RETURN
  138.       END
  139.