CTITLESASPLIN -- CUBIC SPLINE INTERPOLATOR 00010000 C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CA 00020000 CA AUTHOR B. VERWEST 00030000 CA DESIGNER B. VERWEST 00040000 CA LANGUAGE VS FORTRAN 00050000 C WRITTEN 08/21/87 00060000 C REVISED 02/05/90 JJC - RENAMED SPLINE TO SASPLIN. 00070000 C RENAMED SPCOFF TO SASPCOF. 00080000 C REVISED 03/24/92 ESN - CONVERTED TO THE IBM. 00090000 CA 00100000 CA 00110000 CA CALLING SEQUENCE: 00120000 CA 00130000 CA CALL SASPLIN (XN,FN,NPT,XE,DX,LPT,VAL,WORK) 00140000 CA 00150000 CA IN R4 XN AN INDEPENDENT VARIABLE LENGTH NPT 00160000 CA IN R4 FN A DEPENDENT VARIABLE LENGTH NPT 00170000 CA IN I4 NPT NO. OF GIVEN DATA POINTS 00180000 CA IN R4 XE START VALUE OF INDEPENDENT VARIABLE 00190000 CA TO INTERPOLATE 00200000 CA IN R4 DX INCREMENT OF INDEPENDENT VARIABLE 00210000 CA TO INTERPOLATE 00220000 CA IN I4 LPT NO. OF POINTS TO GENERATE ARRAY OF 00230000 CA DEPENDENT VARIABLE 00240000 CA OUT R4 VAL GENERATED DEPENDENT VARIABLE LENGTH LPT 00250000 CA OUT R4 WORK WORK ARRAY LENGTH NPT*2 00260000 CA THE FIRST NPT VALUES IN WORK ARE 00270000 CA THE SECOND DERIVATIVE 00280000 CA 00290000 CA 00300000 CA PURPOSE: CUBIC SASPLIN INTERPOLATION 00310000 CA 00320000 CA EXTRAPOLATION PAST ENDS OF XN IS TURNED OFF AND FN=0.0 00330000 CA (B. VERWEST - AUG. 21, 1987) 00340000 C 00350000 SUBROUTINE SASPLIN (XN,FN,NPT,XE,DX,LPT,VAL,WORK) 00360000 C 00370001 IMPLICIT INTEGER(A-Z) 00380000 C 00390000 REAL A 00400000 REAL B 00410000 REAL DX 00420000 REAL FN 00430000 REAL HL 00440000 REAL R6 00450000 REAL VAL 00460000 REAL WORK 00470000 REAL X 00480000 REAL XE 00490000 REAL XN 00500000 C 00510000 DIMENSION XN(NPT),FN(NPT),VAL(LPT),WORK(2*NPT) 00520000 C 00530000 C SET UP COEFFICIENTS 00540000 C 00550000 R6=1./6. 00560000 CALL SASPCOF(NPT,XN,FN,WORK(1),WORK(NPT+1)) 00570000 C 00580000 C 00590000 NDX = 1 00600000 X = XE 00610000 C 00620000 C INITIALIZE IF TO EXTRAPOLATE ONLY AFTER XN(NPT) 00630000 C 00640000 IF (X .LE. XN(NPT)) GO TO 100 00650000 GO TO 180 00660000 C 00670000 C EXTRAPOLATE UP TO XN(1) 00680000 C 00690000 100 IF (X .GE. XN(1)) GO TO 140 00700000 120 CONTINUE 00710000 VAL(NDX) = 0.0 00720000 NDX = NDX + 1 00730000 X = X + DX 00740000 IF (NDX .GT. LPT) RETURN 00750000 IF (X .LT. XN(1)) GO TO 120 00760000 C 00770000 C INTERPOLATE BETWEEN XN(1) AND XN(NPT) 00780000 C 00790000 140 IF (X. GT. XN(NPT)) GO TO 180 00800000 N = 1 00810000 160 CONTINUE 00820000 IF (N .GT. NPT-1) GO TO 180 00830000 C 00840000 C CHECK IF X IN SEGMENT N - IF YES THEN INTERPOLATE ELSE BUMP N 00850000 C AND TRY AGAIN 00860000 C 00870000 IF (X .GE. XN(N) .AND. X .LE. XN(N+1)) THEN 00880000 A=XN(N+1)-X 00890000 B=X-XN(N) 00900000 HL=XN(N+1)-XN(N) 00910000 VAL(NDX) = A*WORK(N)*(A**2/HL-HL)*R6 00920000 & +B*WORK(N+1)*(B**2/HL-HL)*R6 00930000 & +(A*FN(N)+B*FN(N+1))/HL 00940000 X = X+DX 00950000 NDX = NDX+1 00960000 IF (NDX .GT. LPT) RETURN 00970000 ELSE 00980000 N = N+1 00990000 ENDIF 01000000 GO TO 160 01010000 C 01020000 C EXTRAPOLATE AFTER XN(NPT) 01030000 C 01040000 180 IF (NDX .GT. LPT) RETURN 01050000 DO 200 N=NDX,LPT 01060000 VAL(N) = 0.0 01070000 200 X = X + DX 01080000 RETURN 01090000 END 01100000