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