C DACF 10
C ..................................................................DACF 20
C DACF 30
C SUBROUTINE DACFI DACF 40
C DACF 50
C PURPOSE DACF 60
C TO INTERPOLATE FUNCTION VALUE Y FOR A GIVEN ARGUMENT VALUE DACF 70
C X USING A GIVEN TABLE (ARG,VAL) OF ARGUMENT AND FUNCTION DACF 80
C VALUES. DACF 90
C DACF 100
C USAGE DACF 110
C CALL DACFI (X,ARG,VAL,Y,NDIM,EPS,IER) DACF 120
C DACF 130
C DESCRIPTION OF PARAMETERS DACF 140
C X - DOUBLE PRECISION ARGUMENT VALUE SPECIFIED BY INPUT.DACF 150
C ARG - DOUBLE PRECISION INPUT VECTOR (DIMENSION NDIM) OF DACF 160
C ARGUMENT VALUES OF THE TABLE (POSSIBLY DESTROYED). DACF 170
C VAL - DOUBLE PRECISION INPUT VECTOR (DIMENSION NDIM) OF DACF 180
C FUNCTION VALUES OF THE TABLE (DESTROYED). DACF 190
C Y - RESULTING INTERPOLATED DOUBLE PRECISION FUNCTION DACF 200
C VALUE. DACF 210
C NDIM - AN INPUT VALUE WHICH SPECIFIES THE NUMBER OF DACF 220
C POINTS IN TABLE (ARG,VAL). DACF 230
C EPS - SINGLE PRECISION INPUT CONSTANT WHICH IS USED AS DACF 240
C UPPER BOUND FOR THE ABSOLUTE ERROR. DACF 250
C IER - A RESULTING ERROR PARAMETER. DACF 260
C DACF 270
C REMARKS DACF 280
C (1) TABLE (ARG,VAL) SHOULD REPRESENT A SINGLE-VALUED DACF 290
C FUNCTION AND SHOULD BE STORED IN SUCH A WAY, THAT THE DACF 300
C DISTANCES ABS(ARG(I)-X) INCREASE WITH INCREASING DACF 310
C SUBSCRIPT I. TO GENERATE THIS ORDER IN TABLE (ARG,VAL), DACF 320
C SUBROUTINES DATSG, DATSM OR DATSE COULD BE USED IN A DACF 330
C PREVIOUS STAGE. DACF 340
C (2) NO ACTION BESIDES ERROR MESSAGE IN CASE NDIM LESS DACF 350
C THAN 1. DACF 360
C (3) INTERPOLATION IS TERMINATED EITHER IF THE DIFFERENCE DACF 370
C BETWEEN TWO SUCCESSIVE INTERPOLATED VALUES IS DACF 380
C ABSOLUTELY LESS THAN TOLERANCE EPS, OR IF THE ABSOLUTE DACF 390
C VALUE OF THIS DIFFERENCE STOPS DIMINISHING, OR AFTER DACF 400
C (NDIM-1) STEPS (THE NUMBER OF POSSIBLE STEPS IS DACF 410
C DIMINISHED IF AT ANY STAGE INFINITY ELEMENT APPEARS IN DACF 420
C THE DOWNWARD DIAGONAL OF INVERTED-DIFFERENCES-SCHEME DACF 430
C AND IF IT IS IMPOSSIBLE TO ELIMINATE THIS INFINITY DACF 440
C ELEMENT BY INTERCHANGING OF TABLE POINTS). DACF 450
C FURTHER IT IS TERMINATED IF THE PROCEDURE DISCOVERS TWO DACF 460
C ARGUMENT VALUES IN VECTOR ARG WHICH ARE IDENTICAL. DACF 470
C DEPENDENT ON THESE FOUR CASES, ERROR PARAMETER IER IS DACF 480
C CODED IN THE FOLLOWING FORM DACF 490
C IER=0 - IT WAS POSSIBLE TO REACH THE REQUIRED DACF 500
C ACCURACY (NO ERROR). DACF 510
C IER=1 - IT WAS IMPOSSIBLE TO REACH THE REQUIRED DACF 520
C ACCURACY BECAUSE OF ROUNDING ERRORS. DACF 530
C IER=2 - IT WAS IMPOSSIBLE TO CHECK ACCURACY BECAUSE DACF 540
C NDIM IS LESS THAN 2, OR THE REQUIRED ACCURACY DACF 550
C COULD NOT BE REACHED BY MEANS OF THE GIVEN DACF 560
C TABLE. NDIM SHOULD BE INCREASED. DACF 570
C IER=3 - THE PROCEDURE DISCOVERED TWO ARGUMENT VALUES DACF 580
C IN VECTOR ARG WHICH ARE IDENTICAL. DACF 590
C DACF 600
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED DACF 610
C NONE DACF 620
C DACF 630
C METHOD DACF 640
C INTERPOLATION IS DONE BY CONTINUED FRACTIONS AND INVERTED- DACF 650
C DIFFERENCES-SCHEME. ON RETURN Y CONTAINS AN INTERPOLATED DACF 660
C FUNCTION VALUE AT POINT X, WHICH IS IN THE SENSE OF REMARK DACF 670
C (3) OPTIMAL WITH RESPECT TO GIVEN TABLE. FOR REFERENCE, SEE DACF 680
C F.B.HILDEBRAND, INTRODUCTION TO NUMERICAL ANALYSIS, DACF 690
C MCGRAW-HILL, NEW YORK/TORONTO/LONDON, 1956, PP.395-406. DACF 700
C DACF 710
C ..................................................................DACF 720
C DACF 730
SUBROUTINE DACFI(X,ARG,VAL,Y,NDIM,EPS,IER) DACF 740
C DACF 750
C DACF 760
DIMENSION ARG(1),VAL(1) DACF 770
DOUBLE PRECISION ARG,VAL,X,Y,Z,P1,P2,P3,Q1,Q2,Q3,AUX,H DACF 780
IER=2 DACF 790
IF(NDIM)20,20,1 DACF 800
1 Y=VAL(1) DACF 810
DELT2=0. DACF 820
IF(NDIM-1)20,20,2 DACF 830
C DACF 840
C PREPARATIONS FOR INTERPOLATION LOOP DACF 850
2 P2=1.D0 DACF 860
P3=Y DACF 870
Q2=0.D0 DACF 880
Q3=1.D0 DACF 890
C DACF 900
C DACF 910
C START INTERPOLATION LOOP DACF 920
DO 16 I=2,NDIM DACF 930
II=0 DACF 940
P1=P2 DACF 950
P2=P3 DACF 960
Q1=Q2 DACF 970
Q2=Q3 DACF 980
Z=Y DACF 990
DELT1=DELT2 DACF1000
JEND=I-1 DACF1010
C DACF1020
C COMPUTATION OF INVERTED DIFFERENCES DACF1030
3 AUX=VAL(I) DACF1040
DO 10 J=1,JEND DACF1050
H=VAL(I)-VAL(J) DACF1060
IF(DABS(H)-1.D-13*DABS(VAL(I)))4,4,9 DACF1070
4 IF(ARG(I)-ARG(J))5,17,5 DACF1080
5 IF(J-JEND)8,6,6 DACF1090
C DACF1100
C INTERCHANGE ROW I WITH ROW I+II DACF1110
6 II=II+1 DACF1120
III=I+II DACF1130
IF(III-NDIM)7,7,19 DACF1140
7 VAL(I)=VAL(III) DACF1150
VAL(III)=AUX DACF1160
AUX=ARG(I) DACF1170
ARG(I)=ARG(III) DACF1180
ARG(III)=AUX DACF1190
GOTO 3 DACF1200
C DACF1210
C COMPUTATION OF VAL(I) IN CASE VAL(I)=VAL(J) AND J LESS THAN I-1 DACF1220
8 VAL(I)=1.7D38 DACF1230
GOTO 10 DACF1240
C DACF1250
C COMPUTATION OF VAL(I) IN CASE VAL(I) NOT EQUAL TO VAL(J) DACF1260
9 VAL(I)=(ARG(I)-ARG(J))/H DACF1270
10 CONTINUE DACF1280
C INVERTED DIFFERENCES ARE COMPUTED DACF1290
C DACF1300
C COMPUTATION OF NEW Y DACF1310
P3=VAL(I)*P2+(X-ARG(I-1))*P1 DACF1320
Q3=VAL(I)*Q2+(X-ARG(I-1))*Q1 DACF1330
IF(Q3)11,12,11 DACF1340
11 Y=P3/Q3 DACF1350
GOTO 13 DACF1360
12 Y=1.7D38 DACF1370
13 DELT2=DABS(Z-Y) DACF1380
IF(DELT2-EPS)19,19,14 DACF1390
14 IF(I-10)16,15,15 DACF1400
15 IF(DELT2-DELT1)16,18,18 DACF1410
16 CONTINUE DACF1420
C END OF INTERPOLATION LOOP DACF1430
C DACF1440
C DACF1450
RETURN DACF1460
C DACF1470
C THERE ARE TWO IDENTICAL ARGUMENT VALUES IN VECTOR ARG DACF1480
17 IER=3 DACF1490
RETURN DACF1500
C DACF1510
C TEST VALUE DELT2 STARTS OSCILLATING DACF1520
18 Y=Z DACF1530
IER=1 DACF1540
RETURN DACF1550
C DACF1560
C THERE IS SATISFACTORY ACCURACY WITHIN NDIM-1 STEPS DACF1570
19 IER=0 DACF1580
20 RETURN DACF1590
END DACF1600