* * $Id: splin3.F,v 1.1.1.1 1996/04/01 15:03:14 mclareni Exp $ * * $Log: splin3.F,v $ * Revision 1.1.1.1 1996/04/01 15:03:14 mclareni * Mathlib gen * * #include "sys/CERNLIB_machine.h" #include "_gen/pilot.h" SUBROUTINE SPLIN3(X,Y,DERIV,N,NC,Z,FVALUE,FDERIV,M,MC,IOP) C C CERN LIBRARY PROGRAM NO E-209. C C REVISED VERSION JULY 1973. C C PURPOSE = TO COMPUTE A NATURAL SPLINE APPROXIMATION OF THIRD ORDER C FOR A FUNCTION Y(X) GIVEN IN THE N POINTS (X(I),Y(I)) , C I=1(1)N. C C PARAMETERS (IN LIST). C C X = AN ARRAY STORING THE INPUT ARGUMENTS.DIMENSION X(N). C Y = AN ARRAY STORING THE INPUT FUNCTION VALUES.THE ELEMENT C Y(I) REPRESENT THE FUNCTION VALUE Y(X) FOR X=X(I). C DERIV = AN ARRAY USED FOR STORING THE COMPUTED DERIVATIVES OF C THE FUNCTION Y(X).IN DERIV(I,1) AND DERIV(I,2) ARE STOR- C ED THE FIRST-AND SECOND ORDER DERIVATIVES OF Y(X) FOR C X=X(I) RESPECTIVELY. C N = NUMBER OF INPUT FUNCTION VALUES. C NC = ARRAY DERIV IS DIMENSIONED DERIV(NC,2) IN CALLING C PROGRAM. C Z = AN ARRAY STORING THE ARGUMENTS FOR THE INTERPOLATED C VALUES TO BE COMPUTED. C FVALUE = AN ARRAY STORING THE COMPUTED INTERPOLATED VALUES. C FVALUE(J) REPRESENT THE FUNCTION VALUE FVALUE(Z) FOR C Z=Z(J). C FDERIV = AN ARRAY USED FOR STORING THE DERIVATIVES OF THE COM- C PUTED INTERPOLATED VALUES.EXPLANATION AS FOR DERIV. C M = NUMBER OF INTERPOLATED VALUES TO BE COMPUTED. C MC = ARRAY FDERIV IS DIMENSIONED FDERIV(MC,2) IN CALLING C PROGRAM. C IOP = OPTION PARAMETER.FOR IOP.LE.0 THE DERIVATIVES FOR EACH C SUB-INTERVAL IN THE SPLINE APPROXIMATION ARE COMPUTED. C IOP=-1, THE SECOND ORDER END-POINT C DERIVATIVES ARE COMPUTED BY C LINEAR EXTRAPOLATION. C IOP=0 , THE SECOND ORDER END-POINT C DERIVATIVES ASSUMED TO BE GI- C VEN (SEE COMMON /SPAPPR/). C IOP=1 , COMPUTE SPLINE APPROXIMATIONS C FOR THE ARGUMENTS GIVEN IN C THE ARRAY Z,THE DERIVATIVES C BEEING ASSUMED TO HAVE BEEN C CALCULATED IN A PREVIOUS CALL C ON THE ROUTINE. C C PARAMETERS (IN COMMON BLOCK / SPAPPR /). C C SECD1 = VALUE OF THE SECOND DERIVATIVE D2Y(X)/DX2 FOR THE INPUT C ARGUMENT X=X(1). C SECDN = VALUE OF THE SECOND DERIVATIVE D2Y(X)/DX2 FOR THE INPUT C ARGUMENT X=X(N). C NB. VALUES HAVE TO BE ASSIGNED TO SECD1 AND SECDN IN THE C CALLING PROGRAM.IF A NATURAL SPLINE FIT IS WANTED PUT C SECD1=SECDN=0. C VOFINT = COMPUTED APPROXIMATION FOR THE INTEGRAL OF Y(X) TAKEN C FROM X(1) TO X(N). C IERR = ERROR PARAMETER.IERR=0,NO ERRORS OCCURED. C IERR=1,THE NUMBER OF POINTS TOO SMALL C I.E.N LESS THAN 4. C IERR=2,THE ARGUMENTS X(I) NOT IN INCREA- C SING ORDER. C IERR=3,ARGUMENT TO BE USED IN INTERPOLA- C TION ABOVE RANGE. C IERR=4,ARGUMENT TO BE USED IN INTERPOLA- C TION BELOW RANGE. C NXY = N (SEE ABOVE),HAS TO BE STORED FOR ENTRIES CORRESPONDING C TO IOP=1. C DIMENSION X(NC) , Y(NC) , DERIV(NC,2) , Z(MC) , FVALUE(MC) , 1 FDERIV(MC,2) C COMMON / SPAPPR / SECD1 , SECDN , VOFINT , IERR , NXY DATA ZERO,HALF,ONE,THREE/0.,.5,1.,3./ DATA THIRD , SIXTH / .333333333333333 , .166666666666667 / C C*UL 1000 IF (IOP.GT.0) GO TO 1110 IF (IOP.GT.0) GO TO 1110 C IERR=0 C C CHECK IF ENOUGH DATA-POINTS ARE AVAILABLEI.E. IF N LESS THAN 4 NO C THIRD ORDER SPLINE APPROXIMATION IS POSSIBLE. C IF (N.GE.4) GO TO 1010 C IERR=1 GO TO 2000 C C START CALCULATION OF COEFFICIENTS TO BE USED IN THE SYSTEM OF EQU- C ATIONS FOR THE SECOND ORDER DERIVATIVES OF Y(X). C 1010 IF (IOP.NE.-1) GO TO 1015 SECD1=ZERO SECDN = ZERO BET1=ONE/(ONE+HALF*(X(2)-X(1))/(X(3)-X(2))) ALF1=BET1*(ONE- ((X(2)-X(1))/(X(3)-X(2)))**2) BETN=ONE/(ONE+HALF*(X(N)-X(N-1))/(X(N-1)-X(N-2))) ALFN=BETN*(ONE- ((X(N)-X(N-1))/(X(N-1)-X(N-2)))**2) C 1015 DERIV(1,2)=SECD1 DERIV(N,2)=SECDN DERIV(1,1)=ZERO DXPLUS=X(2)-X(1) C C CHECK IF ARGUMENTS ARE IN INCREASING ORDER.IF NOT PRINT ERROR C MESSAGE AND STOP. C IF ( DXPLUS.GT.ZERO) GO TO 1020 IN=1 IERR=2 GO TO 2000 C 1020 DYPLUS=(Y(2)-Y(1))/DXPLUS IU=N-1 DO 1040 I=2,IU DXMIN =DXPLUS DYMIN =DYPLUS DXPLUS=X(I+1)-X(I) C C CHECK IF ARGUMENTS ARE IN INCREASING ORDER.IF NOT PRINT ERROR C MESSAGE AND STOP. C IF (DXPLUS.GT.ZERO) GO TO 1030 C IN=I IERR=2 GO TO 2000 C 1030 DXINV =ONE/(DXPLUS+DXMIN) DYPLUS=(Y(I+1)-Y(I))/DXPLUS DIVDIF=DXINV*(DYPLUS-DYMIN) ALF =HALF*DXINV*DXMIN BET =HALF-ALF C IF (I.EQ.2) DIVDIF=DIVDIF-THIRD*ALF*DERIV(1,2) IF (I.EQ.IU) DIVDIF=DIVDIF-THIRD*BET*DERIV(N,2) IF (I.EQ.2) ALF=ZERO C IF (IOP.NE.-1) GO TO 1035 IF (I.NE.2) GO TO 1032 BET=BET*ALF1 DIVDIF=DIVDIF*BET1 GO TO 1035 1032 IF (I.NE.IU) GO TO 1035 ALF=ALF*ALFN DIVDIF=DIVDIF*BETN C 1035 DXINV =ONE/(ONE+ALF*DERIV(I-1,1)) DERIV(I,1)=-DXINV*BET DERIV(I,2)= DXINV*(THREE*DIVDIF-ALF*DERIV(I-1,2)) 1040 CONTINUE C C COMPUTE THE SECOND DERIVATIVES BY BACKWARDS RECURRENCE RELATION. C THE SECOND ORDER DERIVATIVES FOR X=X(N-1) ALREADY COMPUTED. C C*UL 1050 DO 1060 I=2,IU DO 1060 I=2,IU J=N-I DERIV(J,2)=DERIV(J,1)*DERIV(J+1,2)+DERIV(J,2) 1060 CONTINUE C IF (IOP.NE.-1) GO TO 1070 DERIV(1,2)=((X(3)-X(1))/(X(3)-X(2)))*DERIV(2,2)-((X(2)-X(1))/(X(3) 1-X(2)))*DERIV(3,2) DERIV(N,2)=-((X(N)-X(N-1))/(X(N-1)-X(N-2)))*DERIV(N-2,2)+((X(N)-X( 1N-2))/(X(N-1)-X(N-2)))*DERIV(N-1,2) C C CALCULATION OF THE SECOND ORDER DERIVATIVES FINISHED.START CAL- C CULATION OF THE FIRST ORDER DERIVATIVES AND OF THE INTEGRAL. C 1070 VOFINT=ZERO DO 1080 I=1,IU DXPLUS=X(I+1)-X(I) DYPLUS=Y(I+1)-Y(I) DIVDIF=DYPLUS/DXPLUS DERIV(I,1)=DIVDIF-DXPLUS*(THIRD*DERIV(I,2)+SIXTH*DERIV(I+1,2)) DXPLUS=HALF*DXPLUS VOFINT=VOFINT+DXPLUS*(Y(I+1)+Y(I)-THIRD*(DERIV(I+1,2)+DERIV(I,2))* 1DXPLUS**2) 1080 CONTINUE C C COMPUTE THE LAST FIRST ORDER DERIVATIVE. C DXPLUS=X(N)-X(N-1) DYPLUS=Y(N)-Y(N-1) DIVDIF=DYPLUS/DXPLUS DERIV(N,1)=DIVDIF+DXPLUS*(SIXTH*DERIV(N-1,2)+THIRD*DERIV(N,2)) C C CALCULATION OF FIRST ORDER DERIVATIVES AND INTEGRAL FINISHED. C C SET VALUE OF N IN COMMON BLOCK / SPAPPR /. C NXY=N C C COMPUTE INTERPOLATED VALUES IF ANY. C 1110 IF (M.LT.1) RETURN C XL=X(1) XU=X(2) IP=3 IL=0 C C--- Transform this DO-loop into a GOTO to avoid illegal jumps into it C C DO 1160 J=1,M J=0 1121 J=J+1 IF(J.GT.M) GOTO 1160 ARG=Z(J) IF (ARG.GT.XU) GO TO 1170 IF (ARG.LT.XL) GO TO 1190 C C ARGUMENT IN CORRECT RANGE.CHECK IF POLYNOMIAL COEFFICIENTS HAVE C TO BE CALCULATED. C C*UL 1130 IF (IL.GT.0) GO TO 1150 IF (IL.GT.0) GO TO 1150 C C COMPUTE POLYNOMIAL COEFFICIENTS. C 1140 II=IP-2 A0=Y(II) A1=DERIV(II,1) A4=DERIV(II,2) A6=(DERIV(II+1,2)-A4)/(XU-XL) A2=HALF*A4 A3=SIXTH*A6 A5=HALF*A6 IL=1 C C CALCULATION OF POLYNOMIAL COEFFICIENTS FINISHED.COMPUTE VALUES. C 1150 ARG=ARG-XL FVALUE(J)=((A3*ARG+A2)*ARG+A1)*ARG+A0 FDERIV(J,1)=(A5*ARG+A4)*ARG+A1 FDERIV(J,2)=A6*ARG+A4 C 1155 CONTINUE C GOTO 1121 1160 CONTINUE C C CALCULATION OF INTERPOLATED VALUES FINISHED. C RETURN C C ARGUMENT ABOVE PRESENT RANGE.SHIFT RANGE UPWARDS. C 1170 IF(IP.GT.NXY) GO TO 1185 IPP=IP DO 1180 I=IPP,NXY IF (ARG.GT.X(I)) GO TO 1180 XL=X(I-1) XU=X(I) IP=I+1 IL=0 GO TO 1140 C 1180 CONTINUE C C ARGUMENT OUT OF RANGE,I.E. ARG GREATER THAN X(N). C 1185 IERR=3 IP=NXY+1 GO TO 2010 C C ARGUMENT BELOW PRESENT RANGE.SHIFT DOWNWARDS. C 1190 IPP=IP DO 1200 I=1,IPP II=IP-I-2 IF (II.EQ.0) GO TO 1210 IF (ARG.LT.X(II)) GO TO 1200 XL=X(II) XU=X(II+1) IP=II+2 IL=0 GO TO 1140 C 1200 CONTINUE C C ARGUMENT OUT OF RANGE,I.E. ARG LESS THAN X(1). C 1210 IERR=4 IP=3 GO TO 2010 C C PRINT ERROR MESSAGES. C 2000 IF (IERR.EQ.1) WRITE(6,3000) IERR IF (IERR.EQ.2) WRITE(6,3000) IERR , X(IN) , X(IN+1) RETURN C 2010 WRITE(6,3000) IERR , ARG C FVALUE(J)=ZERO FDERIV(J,1)=ZERO FDERIV(J,2)=ZERO C II=IP-2 XL=X(II) XU=X(II+1) IL=0 GO TO 1155 C 3000 FORMAT(//5X,'*** SUBROUTINE SPLIN3 ERROR NO ',I2,' ***', 1 2(4X,E21.14)) C END