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) .GT. 0.) ETA = ABS(W) S2(I) = S2(I) + W 160 CONTINUE IF ((ETA-EPSLN) .GE. 0.) GOTO 130 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 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