C DAHI 10
C ..................................................................DAHI 20
C DAHI 30
C SUBROUTINE DAHI DAHI 40
C DAHI 50
C PURPOSE DAHI 60
C TO INTERPOLATE FUNCTION VALUE Y FOR A GIVEN ARGUMENT VALUE DAHI 70
C X USING A GIVEN TABLE (ARG,VAL) OF ARGUMENT, FUNCTION, AND DAHI 80
C DERIVATIVE VALUES. DAHI 90
C DAHI 100
C USAGE DAHI 110
C CALL DAHI (X,ARG,VAL,Y,NDIM,EPS,IER) DAHI 120
C DAHI 130
C DESCRIPTION OF PARAMETERS DAHI 140
C X - DOUBLE PRECISION ARGUMENT VALUE SPECIFIED BY INPUT.DAHI 150
C ARG - DOUBLE PRECISION INPUT VECTOR (DIMENSION NDIM) OF DAHI 160
C ARGUMENT VALUES OF THE TABLE (NOT DESTROYED). DAHI 170
C VAL - DOUBLE PRECISION INPUT VECTOR (DIMENSION 2*NDIM) OFDAHI 180
C FUNCTION AND DERIVATIVE VALUES OF THE TABLE (DES- DAHI 190
C TROYED). FUNCTION AND DERIVATIVE VALUES MUST BE DAHI 200
C STORED IN PAIRS, THAT MEANS BEGINNING WITH FUNCTIONDAHI 210
C VALUE AT POINT ARG(1) EVERY FUNCTION VALUE MUST BE DAHI 220
C FOLLOWED BY THE DERIVATIVE VALUE AT THE SAME POINT.DAHI 230
C Y - RESULTING INTERPOLATED DOUBLE PRECISION FUNCTION DAHI 240
C VALUE. DAHI 250
C NDIM - AN INPUT VALUE WHICH SPECIFIES THE NUMBER OF DAHI 260
C POINTS IN TABLE (ARG,VAL). DAHI 270
C EPS - SINGLE PRECISION INPUT CONSTANT WHICH IS USED AS DAHI 280
C UPPER BOUND FOR THE ABSOLUTE ERROR. DAHI 290
C IER - A RESULTING ERROR PARAMETER. DAHI 300
C DAHI 310
C REMARKS DAHI 320
C (1) TABLE (ARG,VAL) SHOULD REPRESENT A SINGLE-VALUED DAHI 330
C FUNCTION AND SHOULD BE STORED IN SUCH A WAY, THAT THE DAHI 340
C DISTANCES ABS(ARG(I)-X) INCREASE WITH INCREASING DAHI 350
C SUBSCRIPT I. TO GENERATE THIS ORDER IN TABLE (ARG,VAL), DAHI 360
C SUBROUTINES DATSG, DATSM OR DATSE COULD BE USED IN A DAHI 370
C PREVIOUS STAGE. DAHI 380
C (2) NO ACTION BESIDES ERROR MESSAGE IN CASE NDIM LESS DAHI 390
C THAN 1. DAHI 400
C (3) INTERPOLATION IS TERMINATED EITHER IF THE DIFFERENCE DAHI 410
C BETWEEN TWO SUCCESSIVE INTERPOLATED VALUES IS DAHI 420
C ABSOLUTELY LESS THAN TOLERANCE EPS, OR IF THE ABSOLUTE DAHI 430
C VALUE OF THIS DIFFERENCE STOPS DIMINISHING, OR AFTER DAHI 440
C (2*NDIM-2) STEPS. FURTHER IT IS TERMINATED IF THE DAHI 450
C PROCEDURE DISCOVERS TWO ARGUMENT VALUES IN VECTOR ARG DAHI 460
C WHICH ARE IDENTICAL. DEPENDENT ON THESE FOUR CASES, DAHI 470
C ERROR PARAMETER IER IS CODED IN THE FOLLOWING FORM DAHI 480
C IER=0 - IT WAS POSSIBLE TO REACH THE REQUIRED DAHI 490
C ACCURACY (NO ERROR). DAHI 500
C IER=1 - IT WAS IMPOSSIBLE TO REACH THE REQUIRED DAHI 510
C ACCURACY BECAUSE OF ROUNDING ERRORS. DAHI 520
C IER=2 - IT WAS IMPOSSIBLE TO CHECK ACCURACY BECAUSE DAHI 530
C NDIM IS LESS THAN 2, OR THE REQUIRED ACCURACY DAHI 540
C COULD NOT BE REACHED BY MEANS OF THE GIVEN DAHI 550
C TABLE. NDIM SHOULD BE INCREASED. DAHI 560
C IER=3 - THE PROCEDURE DISCOVERED TWO ARGUMENT VALUES DAHI 570
C IN VECTOR ARG WHICH ARE IDENTICAL. DAHI 580
C DAHI 590
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED DAHI 600
C NONE DAHI 610
C DAHI 620
C METHOD DAHI 630
C INTERPOLATION IS DONE BY MEANS OF AITKENS SCHEME OF DAHI 640
C HERMITE INTERPOLATION. ON RETURN Y CONTAINS AN INTERPOLATED DAHI 650
C FUNCTION VALUE AT POINT X, WHICH IS IN THE SENSE OF REMARK DAHI 660
C (3) OPTIMAL WITH RESPECT TO GIVEN TABLE. FOR REFERENCE, SEE DAHI 670
C F.B.HILDEBRAND, INTRODUCTION TO NUMERICAL ANALYSIS, DAHI 680
C MCGRAW-HILL, NEW YORK/TORONTO/LONDON, 1956, PP.314-317, AND DAHI 690
C GERSHINSKY/LEVINE, AITKEN-HERMITE INTERPOLATION, DAHI 700
C JACM, VOL.11, ISS.3 (1964), PP.352-356. DAHI 710
C DAHI 720
C ..................................................................DAHI 730
C DAHI 740
SUBROUTINE DAHI(X,ARG,VAL,Y,NDIM,EPS,IER) DAHI 750
C DAHI 760
C DAHI 770
DIMENSION ARG(1),VAL(1) DAHI 780
DOUBLE PRECISION ARG,VAL,X,Y,H,H1,H2 DAHI 790
IER=2 DAHI 800
H2=X-ARG(1) DAHI 810
IF(NDIM-1)2,1,3 DAHI 820
1 Y=VAL(1)+VAL(2)*H2 DAHI 830
2 RETURN DAHI 840
C DAHI 850
C VECTOR ARG HAS MORE THAN 1 ELEMENT. DAHI 860
C THE FIRST STEP PREPARES VECTOR VAL SUCH THAT AITKEN SCHEME CAN BE DAHI 870
C USED. DAHI 880
3 I=1 DAHI 890
DO 5 J=2,NDIM DAHI 900
H1=H2 DAHI 910
H2=X-ARG(J) DAHI 920
Y=VAL(I) DAHI 930
VAL(I)=Y+VAL(I+1)*H1 DAHI 940
H=H1-H2 DAHI 950
IF(H)4,13,4 DAHI 960
4 VAL(I+1)=Y+(VAL(I+2)-Y)*H1/H DAHI 970
5 I=I+2 DAHI 980
VAL(I)=VAL(I)+VAL(I+1)*H2 DAHI 990
C END OF FIRST STEP DAHI1000
C DAHI1010
C PREPARE AITKEN SCHEME DAHI1020
DELT2=0. DAHI1030
IEND=I-1 DAHI1040
C DAHI1050
C START AITKEN-LOOP DAHI1060
DO 9 I=1,IEND DAHI1070
DELT1=DELT2 DAHI1080
Y=VAL(1) DAHI1090
M=(I+3)/2 DAHI1100
H1=ARG(M) DAHI1110
DO 6 J=1,I DAHI1120
K=I+1-J DAHI1130
L=(K+1)/2 DAHI1140
H=ARG(L)-H1 DAHI1150
IF(H)6,14,6 DAHI1160
6 VAL(K)=(VAL(K)*(X-H1)-VAL(K+1)*(X-ARG(L)))/H DAHI1170
DELT2=DABS(Y-VAL(1)) DAHI1180
IF(DELT2-EPS)11,11,7 DAHI1190
7 IF(I-8)9,8,8 DAHI1200
8 IF(DELT2-DELT1)9,12,12 DAHI1210
9 CONTINUE DAHI1220
C END OF AITKEN-LOOP DAHI1230
C DAHI1240
10 Y=VAL(1) DAHI1250
RETURN DAHI1260
C DAHI1270
C THERE IS SUFFICIENT ACCURACY WITHIN 2*NDIM-2 ITERATION STEPS DAHI1280
11 IER=0 DAHI1290
GOTO 10 DAHI1300
C DAHI1310
C TEST VALUE DELT2 STARTS OSCILLATING DAHI1320
12 IER=1 DAHI1330
RETURN DAHI1340
C DAHI1350
C THERE ARE TWO IDENTICAL ARGUMENT VALUES IN VECTOR ARG DAHI1360
13 Y=VAL(1) DAHI1370
14 IER=3 DAHI1380
RETURN DAHI1390
END DAHI1400