*DECK AVINT SUBROUTINE AVINT (X, Y, N, XLO, XUP, ANS, IERR) C***BEGIN PROLOGUE AVINT C***PURPOSE Integrate a function tabulated at arbitrarily spaced C abscissas using overlapping parabolas. C***LIBRARY SLATEC C***CATEGORY H2A1B2 C***TYPE SINGLE PRECISION (AVINT-S, DAVINT-D) C***KEYWORDS INTEGRATION, QUADRATURE, TABULATED DATA C***AUTHOR Jones, R. E., (SNLA) C***DESCRIPTION C C Abstract C AVINT integrates a function tabulated at arbitrarily spaced C abscissas. The limits of integration need not coincide C with the tabulated abscissas. C C A method of overlapping parabolas fitted to the data is used C provided that there are at least 3 abscissas between the C limits of integration. AVINT also handles two special cases. C If the limits of integration are equal, AVINT returns a result C of zero regardless of the number of tabulated values. C If there are only two function values, AVINT uses the C trapezoid rule. C C Description of Parameters C The user must dimension all arrays appearing in the call list C X(N), Y(N). C C Input-- C X - real array of abscissas, which must be in increasing C order. C Y - real array of functional values. i.e., Y(I)=FUNC(X(I)). C N - the integer number of function values supplied. C N .GE. 2 unless XLO = XUP. C XLO - real lower limit of integration. C XUP - real upper limit of integration. C Must have XLO .LE. XUP. C C Output-- C ANS - computed approximate value of integral C IERR - a status code C --normal code C =1 means the requested integration was performed. C --abnormal codes C =2 means XUP was less than XLO. C =3 means the number of X(I) between XLO and XUP C (inclusive) was less than 3 and neither of the two C special cases described in the Abstract occurred. C No integration was performed. C =4 means the restriction X(I+1) .GT. X(I) was violated. C =5 means the number N of function values was .LT. 2. C ANS is set to zero if IERR=2,3,4,or 5. C C AVINT is documented completely in SC-M-69-335 C Original program from "Numerical Integration" by Davis & C Rabinowitz. C Adaptation and modifications for Sandia Mathematical Program C Library by Rondall E. Jones. C C***REFERENCES R. E. Jones, Approximate integrator of functions C tabulated at arbitrarily spaced abscissas, C Report SC-M-69-335, Sandia Laboratories, 1969. C***ROUTINES CALLED XERMSG C***REVISION HISTORY (YYMMDD) C 690901 DATE WRITTEN C 890831 Modified array declarations. (WRB) C 890831 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) C 900326 Removed duplicate information from DESCRIPTION section. C (WRB) C 920501 Reformatted the REFERENCES section. (WRB) C***END PROLOGUE AVINT C DOUBLE PRECISION R3,RP5,SUM,SYL,SYL2,SYL3,SYU,SYU2,SYU3,X1,X2,X3 1,X12,X13,X23,TERM1,TERM2,TERM3,A,B,C,CA,CB,CC DIMENSION X(*),Y(*) C***FIRST EXECUTABLE STATEMENT AVINT IERR=1 ANS =0.0 IF (XLO-XUP) 3,100,200 3 IF (N.LT.2) GO TO 215 DO 5 I=2,N IF (X(I).LE.X(I-1)) GO TO 210 IF (X(I).GT.XUP) GO TO 6 5 CONTINUE 6 CONTINUE IF (N.GE.3) GO TO 9 C C SPECIAL N=2 CASE SLOPE = (Y(2)-Y(1))/(X(2)-X(1)) FL = Y(1) + SLOPE*(XLO-X(1)) FR = Y(2) + SLOPE*(XUP-X(2)) ANS = 0.5*(FL+FR)*(XUP-XLO) RETURN 9 CONTINUE IF (X(N-2).LT.XLO) GO TO 205 IF (X(3).GT.XUP) GO TO 205 I = 1 10 IF (X(I).GE.XLO) GO TO 15 I = I+1 GO TO 10 15 INLFT = I I = N 20 IF (X(I).LE.XUP) GO TO 25 I = I-1 GO TO 20 25 INRT = I IF ((INRT-INLFT).LT.2) GO TO 205 ISTART = INLFT IF (INLFT.EQ.1) ISTART = 2 ISTOP = INRT IF (INRT.EQ.N) ISTOP = N-1 C R3 = 3.0D0 RP5= 0.5D0 SUM = 0.0 SYL = XLO SYL2= SYL*SYL SYL3= SYL2*SYL C DO 50 I=ISTART,ISTOP X1 = X(I-1) X2 = X(I) X3 = X(I+1) X12 = X1-X2 X13 = X1-X3 X23 = X2-X3 TERM1 = DBLE(Y(I-1))/(X12*X13) TERM2 =-DBLE(Y(I)) /(X12*X23) TERM3 = DBLE(Y(I+1))/(X13*X23) A = TERM1+TERM2+TERM3 B = -(X2+X3)*TERM1 - (X1+X3)*TERM2 - (X1+X2)*TERM3 C = X2*X3*TERM1 + X1*X3*TERM2 + X1*X2*TERM3 IF (I-ISTART) 30,30,35 30 CA = A CB = B CC = C GO TO 40 35 CA = 0.5*(A+CA) CB = 0.5*(B+CB) CC = 0.5*(C+CC) 40 SYU = X2 SYU2= SYU*SYU SYU3= SYU2*SYU SUM = SUM + CA*(SYU3-SYL3)/R3 + CB*RP5*(SYU2-SYL2) + CC*(SYU-SYL) CA = A CB = B CC = C SYL = SYU SYL2= SYU2 SYL3= SYU3 50 CONTINUE SYU = XUP ANS = SUM + CA*(SYU**3-SYL3)/R3 + CB*RP5*(SYU**2-SYL2) 1 + CC*(SYU-SYL) 100 RETURN 200 IERR=2 CALL XERMSG ('SLATEC', 'AVINT', + 'THE UPPER LIMIT OF INTEGRATION WAS NOT GREATER THAN THE ' // + 'LOWER LIMIT.', 4, 1) RETURN 205 IERR=3 CALL XERMSG ('SLATEC', 'AVINT', + 'THERE WERE LESS THAN THREE FUNCTION VALUES BETWEEN THE ' // + 'LIMITS OF INTEGRATION.', 4, 1) RETURN 210 IERR=4 CALL XERMSG ('SLATEC', 'AVINT', + 'THE ABSCISSAS WERE NOT STRICTLY INCREASING. MUST HAVE ' // + 'X(I-1) .LT. X(I) FOR ALL I.', 4, 1) RETURN 215 IERR=5 CALL XERMSG ('SLATEC', 'AVINT', + 'LESS THAN TWO FUNCTION VALUES WERE SUPPLIED.', 4, 1) RETURN END