C
C ****************************************************************
C
FUNCTION SPLIN(T, X, Y, N, NNN)
C
C SPLINE INTERPOLATION OF T IN A GIVEN TABLE OF POINTS(X(I),Y(I)).
C
C Y(I) IS A FUNCTION OF X(I) - X(N)>X(N-1)>....X(1) -OR- X(1)> X(2)
C > ....> X(N) - WHEN X IS NORMALISED X(MIN) = 0, XMAX=1,
C NORM(ABS(X(I+1)-X(I))) < 1.0E-8,(I=1,N-1). THE ARRAY X IS DIVIDED
C INTO OVERLAPPING INTERVALS. EACH INTERVAL CONSISTS OF NNN POINTS,
C THE OVERLAP IS INT(LOG(NNN)+3) POINTS. ( 3<=NNN<=100,N).
C EACH TIME AFTER A CALL WHEN T CHANGES OF INTERVAL OR AFTER A CALL
C WITH ANOTHER SET OF POINTS, A NEW CUBIC SPLINE OF NNN POINTS HAS
C TO BE COMPUTED, THEREFORE THE COMPUTING TIME DEPENDS STRONGLY OF
C THE SUCCESSION OF CALLS.
C FOR EXAMPLE: SUCCESSIVE CALLS, WHERE T IS RANDOMLY CHOSEN WILL COST
C A LOT OF COMPUTING TIME - (T OFTEN CHANGES OF INTERVAL) --->
C TAKE NNN THEN AS LARGE AS POSSIBLE. SUCCESSIVE CALLS WHERE T CHANGES
C SLIGHTLY COST LESS COMPUTING TIME, SO IT WILL BE CLEAR THAT YOUR
C CHOICE OF NNN DEPENDS HEAVILY UPON ITS WAY OF USEAGE.
C
C COMPUTING TIME: FOR THE COMPUTATION OF SPLINE COEFFICIENTS
C THIS ROUTINE NEEDS ABOUT 21*N MULTIPLICATIONS OR DIVISIONS
C FOR THE INTERPOLATION BY MEANS OF THE SPLINE COEFFICIENTS
C THE ROUTINE NEEDS 9 MULTIPLICATIONS OR DIVISIONS.
C
C EXTRAPOLATION: IF T IS NOT IN [X(1),X(N)] T IS LINEAR EXTRAPOLATED
C
C A PART OF THE ALGORITM FOR THE COMPUTATION OF A NATURAL SPLINE
C IS TAKEN FROM T.N.E. GREVILLE, THE LINEAR SYSTEM IS SOLVED BY
C THE METHOD OF SUCCESSIVE OVERRELAXATION. (ERROR = 1.0E-6)
C
C REF: T.N.E. GREVILLE, MATHEMATICS RESEARCH CENTER, U.S. ARMY,
C UNIVERSITY OF WISCONSIN. MATHEMATICAL METHODS.
C
C THIS ROUTINE WAS IMPLEMENTED BY E.B.J. VAN DER ZALM, STERRENWACHT
C UTRECHT (MAY 7TH 1981). ADAPTED FOR VAX\11 BY PAUL KUIN AT OXFORD.
C WRITTEN IN REAL*4 BY P.G. JUDGE, OXFORD
C MODIFIED SO THAT IF THE ARGUMENT IS OUTSIDE THE RANGE OF X-VALUES
C THEN THE FIRST (OR LAST) Y-VALUE IS TAKEN
C
C INPUT: T ARGUMENT (REAL)
C X ARRAY OF ARGUMENTS (REAL)
C Y ARRAY OF FUNCTION VALUES (REAL)
C N LENGTH OF THE ARRAYS X,Y
C IF N = 2 ---> LINEAR INTERPOLATION
C NNN NUMBER OF POINTS OF THE CUBIC SPLINE
C IF NNN > N ---> NNN = N
C IF NNN > 100 ---> NNN = 100
C IF NNN = 2 LINEAR INTERPOLATION ADOPTED.
C
C IF NNN IS ZERO NNN=7 (ALTERED BY P JUDGE)
C (USED TO BE IF NNN OMITTED )
C OUTPUT: SPLIN INTERPOLATED VALUE AT T
C
INCLUDE 'PREC'
DIMENSION X(N), Y(N), S2(100), S3(100), DELY(100), H(100),
* H2(100), DELSQY(100), C(100), B(100), XX(100), YY(100)
LOGICAL SEARCH, FOUND
CHARACTER*4 INIT
DATA INIT/'NOTO'/
SAVE
C
C NN = NUMBER OF POINTS OF THE SPLINE INTERVAL
C
NN = MIN(7,N)
IF (NNN .EQ. 0) GO TO 10
C PGJ IF (%LOC(NNN).EQ.0) GO TO 10
NN = NNN
NN = MIN(N,NNN,100)
C
C IF T IN [X(K)-ERROR,X(K)+ERROR] T = X(K)
C
10 ERROR = 1.0E-12
IF (T.NE.0) ERROR = ABS(1.0E-12*T)
FOUND = SEARCH(X,T,1,N,K,ERROR)
IF (.NOT.FOUND) GOTO 20
SPLIN = Y(K)
RETURN
20 IF (N.NE.2.AND.NN.NE.2) GOTO 30
C
C IF N=2 OR NNN=2 LINEAR INTERPOLATION IS ADOPTED
C
SPLIN = (Y(K+1)-Y(K))*(T-X(K))/(X(K+1)-X(K)) + Y(K)
RETURN
30 IF (N.NE.1) GOTO 40
SPLIN = Y(1)
RETURN
C
C INTERVAL OF THE SPLINE IS COMPUTED
C
40 IOVER = LOG10(FLOAT(NN))*3
INT = (K-IOVER)/(NN-2*IOVER) + 1
I1 = (INT-1)*(NN-2*IOVER)
I2 = I1 + NN -1
IF (I1.GE.1) GOTO 50
I2 = NN
I1 = 1
50 IF (I2.LE.N) GOTO 60
I2 = N
I1 = N - NN + 1
60 ISHIFT = I1 - 1
N1 = NN - 1
C
C INITIALIZATION CHECK
C
IF (INIT.NE.'OKAY') GOTO 80
C
C INTERVAL CHECK
C
DO 70 I=1,NN
IF (XX(I).NE.X(I+ISHIFT)) GOTO 90
IF (YY(I).NE.Y(I+ISHIFT)) GOTO 90
70 CONTINUE
GOTO 200
80 INIT = 'OKAY'
C
C START OF THE COMPUTATION OF SPLINE COEFFICIENTS
C
90 RMIN = 1.0E37
RMAX = -1.0E37
DO 100 I=1,NN
IF (Y(I+ISHIFT).LT.RMIN) RMIN = Y(I+ISHIFT)
IF (Y(I+ISHIFT).GT.RMAX) RMAX = Y(I+ISHIFT)
100 CONTINUE
C
C COMPUTATION OF THE NORM FACTORS XNORM AND YNORM
C
XNORM = 1/(X(NN+ISHIFT)-X(1+ISHIFT))
YNORM = RMAX - RMIN
IF (YNORM.EQ.0) YNORM = RMAX
IF (YNORM.EQ.0) YNORM = 1
YNORM = 1 / YNORM
DO 110 I=1,N1
H(I) = (X(I+1+ISHIFT)-X(I+ISHIFT))*XNORM
DELY(I) = ((Y(I+1+ISHIFT)-Y(I+ISHIFT))/H(I))*YNORM
110 CONTINUE
DO 120 I=2,N1
H2(I) = H(I-1) + H(I)
B(I) = 0.5*H(I-1)/H2(I)
DELSQY(I) = (DELY(I)-DELY(I-1))/H2(I)
S2(I) = 2.*DELSQY(I)
C(I) = 3.*DELSQY(I)
120 CONTINUE
S2(1) = 0
S2(N1+1) = 0
C
C SOLUTION OF THE LINEAR SYSTEM OF THE SPLINE COEFFICIENTS
C BY SUCCESSIVE OVERRELAXATION.
C CONSTANTS: EPS = ERROR CRITERION IN THE ITERATIVE SOLUTION.
C OMEGA = RELAXATION COEFFICIENT.
C
EPSLN = 1.0E-8
OMEGA = 1.0717968
130 ETA = 0.
DO 160 I=2,N1
W = (C(I)-B(I)*S2(I-1)-(0.5-B(I))*S2(I+1)-S2(I))*OMEGA
IF (ABS(W)-ETA) 150, 150, 140
140 ETA = ABS(W)
150 S2(I) = S2(I) + W
160 CONTINUE
IF (ETA-EPSLN) 170, 130, 130
170 DO 180 I=1,N1
S3(I) = (S2(I+1)-S2(I))/H(I)
180 CONTINUE
C
C X AND Y STORED IN XX AND YY FOR INTERVAL CHECK
C
DO 190 I=1,NN
XX(I) = X(I+ISHIFT)
YY(I) = Y(I+ISHIFT)
190 CONTINUE
200 I = K - I1 + 1
HT1 = (T-X(I+ISHIFT))*XNORM
HT2 = (T-X(I+ISHIFT+1) )*XNORM
IF ((T-X(1))*XNORM.GT.0) GOTO 210
C
C EXTRAPOLATION T < X(1)
C
SPLIN = Y(1)
RETURN
210 IF ((T-X(N))*XNORM.GT.0) GOTO 220
C
C INTERPOLATION BY MEANS OF THE SPLINE COEFFICIENTS
C
SS2 = S2(I) + HT1*S3(I)
PROD = HT1*HT2
DELSQS = (S2(I)+S2(I+1)+SS2)/6.
SPLIN = Y(I+ISHIFT) + (HT1*DELY(I)+PROD*DELSQS)/YNORM
RETURN
C
C EXTRAPOLATION T > X(N)
C
220 SPLIN = Y(N)
RETURN
END
C
C **********************************************************************
C
LOGICAL FUNCTION SEARCH(ARR,X,IA,IB,K,ERR)
C
C * SEARCHES THE POINT X WITHIN AN ERROR BOUND -ERR- IN ARRAY ARR *
C * IF THAT POINT IS NOT FOUND, IT GIVES THE *
C * INTERVAL IN THE ARRAY, WHERE THAT POINT FITS. *
C * THE INTERVAL K IF X IN [ARR(K),ARR(K+1)), INTERVAL 1 IF X IN *
C * (-INF,ARR(2)), INTERVAL N-1 IF X IN [ARR(N-1),INF)-(ARR INCREASING)*
C * INTERVAL 1 IF X IN (INF,ARR(2)),INTERVAL N-1 IF X IN [ARR(N-1), *
C * -INF)-(ARR DECREASING) *
C * IN THE WORST CASE A SEARCH COSTS 2*LOG2(N) CYCLES. *
C * THIS ROUTINE IS VERY ECONOMIC, WHEN IN A SEQUENCE OF SEARCHES, *
C * THE DIFFERENCES BETWEEN THE POINTS WHICH ARE TO BE SEARCHED ARE *
C * SMALL. *
C * FOR EXAMPLE: A SEQUENCE OF SEARCHES WHICH GOES ALONG THE ARRAY, *
C * OR A SEQUENCE OF SEARCHES WHICH REMAINS IN A SMALL AREA OF *
C * THE ARRAY. *
C * *
C * INPUT: - ARR ARRAY WITH A NON DECREASING OR A NON *
C * INCREASING FUCNTION. *
C * - IA,IB LOWER AND UPPER LIMIT OF THE ARRAY WITHIN *
C * IS SEARCHED. *
C * - X VALUE WHICH YOU WANT TO SEARCH *
C * - K LAST INDEX WHICH WAS FOUND *
C * SEARCH *
C * - ERR ERROR BOUND *
C * *
C * OUTPUT: - K INDEX OF THE POINT OR INTERVAL OF THE *
C * ARRAY WHICH IS FOUND. *
C * - SEARCH TRUE - X IS FOUND WITHIN THE ERRORBOUND *
C * -ERR- *
C * FALSE - X IS NOT FOUND. *
C * *
C * SEND IN BY E.V.D.ZALM, UTRECHT, STERREWACHT, MARCH 24TH 1981 *
C **********************************************************************
INCLUDE 'PREC'
DIMENSION ARR(IB)
C
ONE=1.0
SEARCH = .TRUE.
IF (ARR(IB)-ARR(IA).LT.0.0) GOTO 100
5 IF (X.LT.ARR(IA).OR.X.GT.ARR(IB)) GOTO 25
IF (K.LT.IA.OR.K.GT.IB-1) GOTO 25
L = K
IP = SIGN(ONE,X-ARR(K))
IF (IP.LT.0.) GOTO 20
10 L = MIN(K+IP,IB)
IF (ARR(L).GE.X) GOTO 30
IP = IP * 2
K = L
GOTO 10
20 K = MAX(L+IP,IA)
IF (ARR(K).LE.X) GOTO 30
IP = IP * 2
L = K
GOTO 20
25 K = IA
L = IB
30 IF ((L-K).LE.1) GOTO 50
I = INT((K+L)/2.)
IF (ARR(I).GE.X) GOTO 40
K = I
GOTO 30
40 L = I
GOTO 30
100 IF (X.LT.ARR(IB).OR.X.GT.ARR(IA)) GOTO 250
IF (K.LT.IA.OR.K.GT.IB-1) GOTO 250
L = K
IP = SIGN(ONE,ARR(K)-X)
IF (IP.LT.0.) GOTO 200
110 L = MIN(K+IP,IB)
IF (ARR(L).LE.X) GOTO 300
IP = IP * 2
K = L
GOTO 110
200 K = MAX(L+IP,IA)
IF (ARR(K).GE.X) GOTO 300
IP = IP * 2
L = K
GOTO 200
250 K = IA
L = IB
300 IF ((L-K).LE.1) GOTO 50
I = INT((K+L)/2.)
IF (ARR(I).LE.X) GOTO 400
K = I
GOTO 300
400 L = I
GOTO 300
50 IF (ABS(X-ARR(K)).LE.ERR) RETURN
IF (ABS(X-ARR(L)).GT.ERR) GOTO 60
K = L
RETURN
60 SEARCH = .FALSE.
RETURN
END