C #SPLINE V1C 16-MAY-73. C LAST UPDATE 30-MAY-73. C FUNCTION SPLINE ( FN, GN, XN, N, X, I) DIMENSION FN(N), GN(N), XN(N) C C *PURPOSE. C TO PERFORM HERMITE (CUBIC) INTERPOLATION OF A FUNCTION C GIVEN BY ITS VALUES AND DERIVATIVES AT A NUMBER OF POINTS. C THIS IS A STANDARD WAY OF REPRESENTING SPLINE FUNCTIONS. C C *PARAMETERS: C C FN - REAL ARRAY OF FUNCTION VALUES FN(I) AT XN(I) FOR I=1 TO N C GN - REAL ARRAY OF DERIVATIVE VALUES GN(I) AT XN(I) FOR I=1 TO N C XN - REAL ARRAY OF X POSITIONS - KNOTS IN CASE OF SPLINE. C N - INTEGER NUMBER OF FUNCTION POINTS STORED - NUMBER OF KNOTS C X - REAL X VALUE AT WHICH INTERPOLATION REQUIRED. C I - INTEGER INDEX TO KNOT INTERVAL USED, WILL BE RETURNED IN RANGE C 1 TO N. MUST BE SET TO ZERO ON FIRST CALL AND WHENEVER C FN OR GN CHANGED. C C INPUT PARAMETERS: C FN, GN, XN, N, X C OUTPUT PARAMETERS: C I, (FUNCTION VALUE) C C *METHOD. C THE ROUTINE SEARCHES FROM THE "I" VALUE LAST USED - C UNLESS I=0 - TO FIND THE CORRECT INTERVAL. IT THEN PERFORMS C A STANDARD HERMITE INTERPOLATION. MUCH OF THE SEARCHING CAN C BE ELIMINATED IN THIS WAY WHEN SUCCESSIVE X VALUES ARE CLOSE, C AS CAN THE CALCULATION OF THE CUBIC COEFFICIENTS WHEN THE C SAME INTERVAL IS USED AGAIN. C THE VALUE RETURNED WHEN THE X CO-ORDINATE IS OUT OF RANGE C IS THAT OBTAINED BY LINEAR INTERPOLATION FROM THE CLOSEST C VALUE AND DERIVATIVE AT AN ENDPOINT. C THE SEARCH ASSUMES THE CO-ORDINATE VALUES XN(I) ARE C IN INCREASING ORDER! C C *ACCURACY. C C *RESTRICTIONS. C THE ARRAY XN(I) MUST BE ORDERED. C C *ERROR CONDITIONS. C THE BEHAVIOUR IS NOT DEFINED IF THE ARRAY XN(I) IS NOT C ORDERED PROPERLY, AND NO INDICATION IS GIVEN OF THIS. THE C OVERFLOW CONDITION - I.E. X > XN(N) - IS SHOWN BY I=N, BUT C THE UNDERFLOW CONDITION IS NOT SHOWN BY I=0 AS THIS IS THE C INDICATOR FOR PRIMING THE INTERNAL VARIABLES. C C *NON STANDARD ROUTINES CALLED. C C *TYPICAL TIMES. C C *ORIGIN. M.R.MANNING. C C *COMMENTS. C C #END. C DATA ILAST/ -1/ C C PHASE 1 - LOCATE INTERVAL TO INTERPOLATE IN. C IF (I .GT. 0) GO TO 5 I= 1 ILAST= -1 5 DX= X - XN(I) IF(DX) 10, 50, 20 C TOO HIGH IN XN TABLE - STEP BACK IF POSSIBLE. 10 I= I - 1 IF( I .LE. 0) GO TO 50 DX= X - XN(I) IF(DX) 10, 50, 50 C C TOO LOW IN XN TABLE - MOVE UP TO NEXT ENTRY 15 I= I+ 1 DX= DEL C TOO LOW IN XN TABLE - CHECK NEXT ENTRY UP. 20 IF( I .GE. N) GO TO 50 DEL= X - XN(I+1) IF(DEL) 50, 15, 15 C C FOUND INTERVAL XN(I) TO XN(I+1), I=0 IS UNDERFLOW AND C I=N IS OVERFLOW. C PHASE 2 - CALCULATE COEFFICIENTS OF LOCAL CUBIC IF NECC. C 50 IF( I .EQ. ILAST) GO TO 100 C ILAST= I IF( I .GT. 0) GO TO 60 C0= FN(1) C1= GN(1) GO TO 65 C 60 C0=FN(I) C1= GN(I) IF( I .LT. N) GO TO 70 65 C2= 0. C3= 0. GO TO 100 C 70 SX= XN(I+1) - XN(I) C2= ( FN(I+1) - C0)/SX C3= GN(I+1) + C1 - C2 - C2 C2= (C2 - C1 - C3)/SX C3= C3/SX/SX C C NOW HAVE CUBIC COEFFICIENTS C PHASE 3 - EVALUATE CUBIC AND RETURN. C 100 SPLINE= C0 + DX*(C1 + DX*(C2 + DX*C3)) C C DO NOT RETURN I=0 UNDERFLOW INDICATION AS THIS CONFUSING C WITH INITIALIZATION. I= MAX0(I, 1) RETURN END