home *** CD-ROM | disk | FTP | other *** search
/ Frostbyte's 1980s DOS Shareware Collection / floppyshareware.zip / floppyshareware / DOOG / PCSSP1.ZIP / ITRPAPSM.ZIP / ACFI.FOR next >
Text File  |  1986-01-03  |  5KB  |  164 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C     SUBROUTINE ACFI
  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 AND FUNCTION
  9. C        VALUES.
  10. C
  11. C     USAGE
  12. C        CALL ACFI (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 (POSSIBLY DESTROYED).
  18. C        VAL    - THE INPUT VECTOR (DIMENSION NDIM) OF FUNCTION
  19. C             VALUES OF THE TABLE (DESTROYED).
  20. C        Y       - THE RESULTING INTERPOLATED FUNCTION VALUE.
  21. C        NDIM   - AN INPUT VALUE WHICH SPECIFIES THE NUMBER OF
  22. C             POINTS IN TABLE (ARG,VAL).
  23. C        EPS    - AN INPUT CONSTANT WHICH IS USED AS UPPER BOUND
  24. C             FOR THE ABSOLUTE ERROR.
  25. C        IER    - A RESULTING ERROR PARAMETER.
  26. C
  27. C     REMARKS
  28. C        (1) TABLE (ARG,VAL) SHOULD REPRESENT A SINGLE-VALUED
  29. C        FUNCTION AND SHOULD BE STORED IN SUCH A WAY, THAT THE
  30. C        DISTANCES ABS(ARG(I)-X) INCREASE WITH INCREASING
  31. C        SUBSCRIPT I. TO GENERATE THIS ORDER IN TABLE (ARG,VAL),
  32. C        SUBROUTINES ATSG, ATSM OR ATSE COULD BE USED IN A
  33. C        PREVIOUS STAGE.
  34. C        (2) NO ACTION BESIDES ERROR MESSAGE IN CASE NDIM LESS
  35. C        THAN 1.
  36. C        (3) INTERPOLATION IS TERMINATED EITHER IF THE DIFFERENCE
  37. C        BETWEEN TWO SUCCESSIVE INTERPOLATED VALUES IS
  38. C        ABSOLUTELY LESS THAN TOLERANCE EPS, OR IF THE ABSOLUTE
  39. C        VALUE OF THIS DIFFERENCE STOPS DIMINISHING, OR AFTER
  40. C        (NDIM-1) STEPS (THE NUMBER OF POSSIBLE STEPS IS
  41. C        DIMINISHED IF AT ANY STAGE INFINITY ELEMENT APPEARS IN
  42. C        THE DOWNWARD DIAGONAL OF INVERTED-DIFFERENCES-SCHEME
  43. C        AND IF IT IS IMPOSSIBLE TO ELIMINATE THIS INFINITY
  44. C        ELEMENT BY INTERCHANGING OF TABLE POINTS).
  45. C        FURTHER IT IS TERMINATED IF THE PROCEDURE DISCOVERS TWO
  46. C        ARGUMENT VALUES IN VECTOR ARG WHICH ARE IDENTICAL.
  47. C        DEPENDENT ON THESE FOUR CASES, ERROR PARAMETER IER IS
  48. C        CODED IN THE FOLLOWING FORM
  49. C         IER=0 - IT WAS POSSIBLE TO REACH THE REQUIRED
  50. C             ACCURACY (NO ERROR).
  51. C         IER=1 - IT WAS IMPOSSIBLE TO REACH THE REQUIRED
  52. C             ACCURACY BECAUSE OF ROUNDING ERRORS.
  53. C         IER=2 - IT WAS IMPOSSIBLE TO CHECK ACCURACY BECAUSE
  54. C             NDIM IS LESS THAN 2, OR THE REQUIRED ACCURACY
  55. C             COULD NOT BE REACHED BY MEANS OF THE GIVEN
  56. C             TABLE. NDIM SHOULD BE INCREASED.
  57. C         IER=3 - THE PROCEDURE DISCOVERED TWO ARGUMENT VALUES
  58. C             IN VECTOR ARG WHICH ARE IDENTICAL.
  59. C
  60. C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  61. C        NONE
  62. C
  63. C     METHOD
  64. C        INTERPOLATION IS DONE BY CONTINUED FRACTIONS AND INVERTED-
  65. C        DIFFERENCES-SCHEME. ON RETURN Y CONTAINS AN INTERPOLATED
  66. C        FUNCTION VALUE AT POINT X, WHICH IS IN THE SENSE OF REMARK
  67. C        (3) OPTIMAL WITH RESPECT TO GIVEN TABLE. FOR REFERENCE, SEE
  68. C        F.B.HILDEBRAND, INTRODUCTION TO NUMERICAL ANALYSIS,
  69. C        MCGRAW-HILL, NEW YORK/TORONTO/LONDON, 1956, PP.395-406.
  70. C
  71. C     ..................................................................
  72. C
  73.       SUBROUTINE ACFI(X,ARG,VAL,Y,NDIM,EPS,IER)
  74. C
  75. C
  76.       DIMENSION ARG(1),VAL(1)
  77.       IER=2
  78.       IF(NDIM)20,20,1
  79.     1 Y=VAL(1)
  80.       DELT2=0.
  81.       IF(NDIM-1)20,20,2
  82. C
  83. C     PREPARATIONS FOR INTERPOLATION LOOP
  84.     2 P2=1.
  85.       P3=Y
  86.       Q2=0.
  87.       Q3=1.
  88. C
  89. C
  90. C     START INTERPOLATION LOOP
  91.       DO 16 I=2,NDIM
  92.       II=0
  93.       P1=P2
  94.       P2=P3
  95.       Q1=Q2
  96.       Q2=Q3
  97.       Z=Y
  98.       DELT1=DELT2
  99.       JEND=I-1
  100. C
  101. C     COMPUTATION OF INVERTED DIFFERENCES
  102.     3 AUX=VAL(I)
  103.       DO 10 J=1,JEND
  104.       H=VAL(I)-VAL(J)
  105.       IF(ABS(H)-1.E-6*ABS(VAL(I)))4,4,9
  106.     4 IF(ARG(I)-ARG(J))5,17,5
  107.     5 IF(J-JEND)8,6,6
  108. C
  109. C     INTERCHANGE ROW I WITH ROW I+II
  110.     6 II=II+1
  111.       III=I+II
  112.       IF(III-NDIM)7,7,19
  113.     7 VAL(I)=VAL(III)
  114.       VAL(III)=AUX
  115.       AUX=ARG(I)
  116.       ARG(I)=ARG(III)
  117.       ARG(III)=AUX
  118.       GOTO 3
  119. C
  120. C     COMPUTATION OF VAL(I) IN CASE VAL(I)=VAL(J) AND J LESS THAN I-1
  121.     8 VAL(I)=1.E38
  122.       GOTO 10
  123. C
  124. C     COMPUTATION OF VAL(I) IN CASE VAL(I) NOT EQUAL TO VAL(J)
  125.     9 VAL(I)=(ARG(I)-ARG(J))/H
  126.    10 CONTINUE
  127. C     INVERTED DIFFERENCES ARE COMPUTED
  128. C
  129. C     COMPUTATION OF NEW Y
  130.       P3=VAL(I)*P2+(X-ARG(I-1))*P1
  131.       Q3=VAL(I)*Q2+(X-ARG(I-1))*Q1
  132.       IF(Q3)11,12,11
  133.    11 Y=P3/Q3
  134.       GOTO 13
  135.    12 Y=1.E38
  136.    13 DELT2=ABS(Z-Y)
  137.       IF(DELT2-EPS)19,19,14
  138.    14 IF(I-8)16,15,15
  139.    15 IF(DELT2-DELT1)16,18,18
  140.    16 CONTINUE
  141. C     END OF INTERPOLATION LOOP
  142. C
  143. C
  144.       RETURN
  145. C
  146. C     THERE ARE TWO IDENTICAL ARGUMENT VALUES IN VECTOR ARG
  147.    17 IER=3
  148.       RETURN
  149. C
  150. C     TEST VALUE DELT2 STARTS OSCILLATING
  151.    18 Y=Z
  152.       IER=1
  153.       RETURN
  154. C
  155. C     THERE IS SATISFACTORY ACCURACY WITHIN NDIM-1 STEPS
  156.    19 IER=0
  157.    20 RETURN
  158.       END
  159. 
  160. COMPUTED
  161. C
  162. C     COMPUTATION OF NEW Y
  163.       P3=VAL(I)*P2+(X-ARG(I-1))*P1
  164.       Q3