1*DECK AVINT 2 SUBROUTINE AVINT (X, Y, N, XLO, XUP, ANS, IERR) 3C***BEGIN PROLOGUE AVINT 4C***PURPOSE Integrate a function tabulated at arbitrarily spaced 5C abscissas using overlapping parabolas. 6C***LIBRARY SLATEC 7C***CATEGORY H2A1B2 8C***TYPE SINGLE PRECISION (AVINT-S, DAVINT-D) 9C***KEYWORDS INTEGRATION, QUADRATURE, TABULATED DATA 10C***AUTHOR Jones, R. E., (SNLA) 11C***DESCRIPTION 12C 13C Abstract 14C AVINT integrates a function tabulated at arbitrarily spaced 15C abscissas. The limits of integration need not coincide 16C with the tabulated abscissas. 17C 18C A method of overlapping parabolas fitted to the data is used 19C provided that there are at least 3 abscissas between the 20C limits of integration. AVINT also handles two special cases. 21C If the limits of integration are equal, AVINT returns a result 22C of zero regardless of the number of tabulated values. 23C If there are only two function values, AVINT uses the 24C trapezoid rule. 25C 26C Description of Parameters 27C The user must dimension all arrays appearing in the call list 28C X(N), Y(N). 29C 30C Input-- 31C X - real array of abscissas, which must be in increasing 32C order. 33C Y - real array of functional values. i.e., Y(I)=FUNC(X(I)). 34C N - the integer number of function values supplied. 35C N .GE. 2 unless XLO = XUP. 36C XLO - real lower limit of integration. 37C XUP - real upper limit of integration. 38C Must have XLO .LE. XUP. 39C 40C Output-- 41C ANS - computed approximate value of integral 42C IERR - a status code 43C --normal code 44C =1 means the requested integration was performed. 45C --abnormal codes 46C =2 means XUP was less than XLO. 47C =3 means the number of X(I) between XLO and XUP 48C (inclusive) was less than 3 and neither of the two 49C special cases described in the Abstract occurred. 50C No integration was performed. 51C =4 means the restriction X(I+1) .GT. X(I) was violated. 52C =5 means the number N of function values was .LT. 2. 53C ANS is set to zero if IERR=2,3,4,or 5. 54C 55C AVINT is documented completely in SC-M-69-335 56C Original program from "Numerical Integration" by Davis & 57C Rabinowitz. 58C Adaptation and modifications for Sandia Mathematical Program 59C Library by Rondall E. Jones. 60C 61C***REFERENCES R. E. Jones, Approximate integrator of functions 62C tabulated at arbitrarily spaced abscissas, 63C Report SC-M-69-335, Sandia Laboratories, 1969. 64C***ROUTINES CALLED XERMSG 65C***REVISION HISTORY (YYMMDD) 66C 690901 DATE WRITTEN 67C 890831 Modified array declarations. (WRB) 68C 890831 REVISION DATE from Version 3.2 69C 891214 Prologue converted to Version 4.0 format. (BAB) 70C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) 71C 900326 Removed duplicate information from DESCRIPTION section. 72C (WRB) 73C 920501 Reformatted the REFERENCES section. (WRB) 74C***END PROLOGUE AVINT 75C 76 DOUBLE PRECISION R3,RP5,SUM,SYL,SYL2,SYL3,SYU,SYU2,SYU3,X1,X2,X3 77 1,X12,X13,X23,TERM1,TERM2,TERM3,A,B,C,CA,CB,CC 78 DIMENSION X(*),Y(*) 79C***FIRST EXECUTABLE STATEMENT AVINT 80 IERR=1 81 ANS =0.0 82 IF (XLO-XUP) 3,100,200 83 3 IF (N.LT.2) GO TO 215 84 DO 5 I=2,N 85 IF (X(I).LE.X(I-1)) GO TO 210 86 IF (X(I).GT.XUP) GO TO 6 87 5 CONTINUE 88 6 CONTINUE 89 IF (N.GE.3) GO TO 9 90C 91C SPECIAL N=2 CASE 92 SLOPE = (Y(2)-Y(1))/(X(2)-X(1)) 93 FL = Y(1) + SLOPE*(XLO-X(1)) 94 FR = Y(2) + SLOPE*(XUP-X(2)) 95 ANS = 0.5*(FL+FR)*(XUP-XLO) 96 RETURN 97 9 CONTINUE 98 IF (X(N-2).LT.XLO) GO TO 205 99 IF (X(3).GT.XUP) GO TO 205 100 I = 1 101 10 IF (X(I).GE.XLO) GO TO 15 102 I = I+1 103 GO TO 10 104 15 INLFT = I 105 I = N 106 20 IF (X(I).LE.XUP) GO TO 25 107 I = I-1 108 GO TO 20 109 25 INRT = I 110 IF ((INRT-INLFT).LT.2) GO TO 205 111 ISTART = INLFT 112 IF (INLFT.EQ.1) ISTART = 2 113 ISTOP = INRT 114 IF (INRT.EQ.N) ISTOP = N-1 115C 116 R3 = 3.0D0 117 RP5= 0.5D0 118 SUM = 0.0 119 SYL = XLO 120 SYL2= SYL*SYL 121 SYL3= SYL2*SYL 122C 123 DO 50 I=ISTART,ISTOP 124 X1 = X(I-1) 125 X2 = X(I) 126 X3 = X(I+1) 127 X12 = X1-X2 128 X13 = X1-X3 129 X23 = X2-X3 130 TERM1 = DBLE(Y(I-1))/(X12*X13) 131 TERM2 =-DBLE(Y(I)) /(X12*X23) 132 TERM3 = DBLE(Y(I+1))/(X13*X23) 133 A = TERM1+TERM2+TERM3 134 B = -(X2+X3)*TERM1 - (X1+X3)*TERM2 - (X1+X2)*TERM3 135 C = X2*X3*TERM1 + X1*X3*TERM2 + X1*X2*TERM3 136 IF (I-ISTART) 30,30,35 137 30 CA = A 138 CB = B 139 CC = C 140 GO TO 40 141 35 CA = 0.5*(A+CA) 142 CB = 0.5*(B+CB) 143 CC = 0.5*(C+CC) 144 40 SYU = X2 145 SYU2= SYU*SYU 146 SYU3= SYU2*SYU 147 SUM = SUM + CA*(SYU3-SYL3)/R3 + CB*RP5*(SYU2-SYL2) + CC*(SYU-SYL) 148 CA = A 149 CB = B 150 CC = C 151 SYL = SYU 152 SYL2= SYU2 153 SYL3= SYU3 154 50 CONTINUE 155 SYU = XUP 156 ANS = SUM + CA*(SYU**3-SYL3)/R3 + CB*RP5*(SYU**2-SYL2) 157 1 + CC*(SYU-SYL) 158 100 RETURN 159 200 IERR=2 160 CALL XERMSG ('SLATEC', 'AVINT', 161 + 'THE UPPER LIMIT OF INTEGRATION WAS NOT GREATER THAN THE ' // 162 + 'LOWER LIMIT.', 4, 1) 163 RETURN 164 205 IERR=3 165 CALL XERMSG ('SLATEC', 'AVINT', 166 + 'THERE WERE LESS THAN THREE FUNCTION VALUES BETWEEN THE ' // 167 + 'LIMITS OF INTEGRATION.', 4, 1) 168 RETURN 169 210 IERR=4 170 CALL XERMSG ('SLATEC', 'AVINT', 171 + 'THE ABSCISSAS WERE NOT STRICTLY INCREASING. MUST HAVE ' // 172 + 'X(I-1) .LT. X(I) FOR ALL I.', 4, 1) 173 RETURN 174 215 IERR=5 175 CALL XERMSG ('SLATEC', 'AVINT', 176 + 'LESS THAN TWO FUNCTION VALUES WERE SUPPLIED.', 4, 1) 177 RETURN 178 END 179