[x,y,n,xlo,xup,ansmlv,ierr]=avint(x,y,n,xlo,xup,ansmlv,ierr);

function [x,y,n,xlo,xup,ansmlv,ierr]=avint(x,y,n,xlo,xup,ansmlv,ierr);
persistent a b c ca cb cc fl fr i inlft inrt istart istop r3 rp5 slope summlv syl syl2 syl3 syu syu2 syu3 term1 term2 term3 x1 x12 x13 x2 x23 x3 ; 

if isempty(fl), fl=0; end;
if isempty(fr), fr=0; end;
if isempty(slope), slope=0; end;
if isempty(i), i=0; end;
if isempty(inlft), inlft=0; end;
if isempty(inrt), inrt=0; end;
if isempty(istart), istart=0; end;
if isempty(istop), istop=0; end;
%***BEGIN PROLOGUE  AVINT
%***PURPOSE  Integrate a function tabulated at arbitrarily spaced
%            abscissas using overlapping parabolas.
%***LIBRARY   SLATEC
%***CATEGORY  H2A1B2
%***TYPE      SINGLE PRECISION (AVINT-S, DAVINT-D)
%***KEYWORDS  INTEGRATION, QUADRATURE, TABULATED DATA
%***AUTHOR  Jones, R. E., (SNLA)
%***DESCRIPTION
%
%     Abstract
%         AVINT integrates a function tabulated at arbitrarily spaced
%         abscissas.  The limits of integration need not coincide
%         with the tabulated abscissas.
%
%         A method of overlapping parabolas fitted to the data is used
%         provided that there are at least 3 abscissas between the
%         limits of integration.  AVINT also handles two special cases.
%         If the limits of integration are equal, AVINT returns a result
%         of zero regardless of the number of tabulated values.
%         If there are only two function values, AVINT uses the
%         trapezoid rule.
%
%     Description of Parameters
%         The user must dimension all arrays appearing in the call list
%              X(N), Y(N).
%
%         Input--
%         X    - real array of abscissas, which must be in increasing
%                order.
%         Y    - real array of functional values. i.e., Y(I)=FUNC(X(I)).
%         N    - the integer number of function values supplied.
%                N .GE. 2 unless XLO = XUP.
%         XLO  - real lower limit of integration.
%         XUP  - real upper limit of integration.
%                Must have XLO .LE. XUP.
%
%         Output--
%         ANS  - computed approximate value of integral
%         IERR - a status code
%              --normal code
%                =1 means the requested integration was performed.
%              --abnormal codes
%                =2 means XUP was less than XLO.
%                =3 means the number of X(I) between XLO and XUP
%                   (inclusive) was less than 3 and neither of the two
%                   special cases described in the Abstract occurred.
%                   No integration was performed.
%                =4 means the restriction X(I+1) .GT. X(I) was violated.
%                =5 means the number N of function values was .LT. 2.
%                ANS is set to zero if IERR=2,3,4,or 5.
%
%     AVINT is documented completely in SC-M-69-335
%     Original program from 'Numerical Integration' by Davis%     Rabinowitz.
%     Adaptation and modifications for Sandia Mathematical Program
%     Library by Rondall E. Jones.
%
%***REFERENCES  R. E. Jones, Approximate integrator of functions
%                 tabulated at arbitrarily spaced abscissas,
%                 Report SC-M-69-335, Sandia Laboratories, 1969.
%***ROUTINES CALLED  XERMSG
%***REVISION HISTORY  (YYMMDD)
%   690901  DATE WRITTEN
%   890831  Modified array declarations.  (WRB)
%   890831  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
%   900326  Removed duplicate information from DESCRIPTION section.
%           (WRB)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  AVINT
%
if isempty(r3), r3=0; end;
if isempty(rp5), rp5=0; end;
if isempty(summlv), summlv=0; end;
if isempty(syl), syl=0; end;
if isempty(syl2), syl2=0; end;
if isempty(syl3), syl3=0; end;
if isempty(syu), syu=0; end;
if isempty(syu2), syu2=0; end;
if isempty(syu3), syu3=0; end;
if isempty(x1), x1=0; end;
if isempty(x2), x2=0; end;
if isempty(x3), x3=0; end;
if isempty(x12), x12=0; end;
if isempty(x13), x13=0; end;
if isempty(x23), x23=0; end;
if isempty(term1), term1=0; end;
if isempty(term2), term2=0; end;
if isempty(term3), term3=0; end;
if isempty(a), a=0; end;
if isempty(b), b=0; end;
if isempty(c), c=0; end;
if isempty(ca), ca=0; end;
if isempty(cb), cb=0; end;
if isempty(cc), cc=0; end;
x_shape=size(x);x=reshape(x,1,[]);
y_shape=size(y);y=reshape(y,1,[]);
%***FIRST EXECUTABLE STATEMENT  AVINT
ierr = 1;
ansmlv = 0.0;
if( xlo<xup )
if( n<2 )
ierr = 5;
xermsg('SLATEC','AVINT','LESS THAN TWO FUNCTION VALUES WERE SUPPLIED.',4,1);
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
else;
for i = 2 : n;
if( x(i)<=x(i-1) )
ierr = 4;
xermsg('SLATEC','AVINT',['THE ABSCISSAS WERE NOT STRICTLY INCREASING.  MUST HAVE ','X(I-1) .LT. X(I) FOR ALL I.'],4,1);
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
elseif( x(i)>xup ) ;
break;
end;
end;
if( n>=3 )
if( x(n-2)<xlo )
ierr = 3;
xermsg('SLATEC','AVINT',['THERE WERE LESS THAN THREE FUNCTION VALUES BETWEEN THE ','LIMITS OF INTEGRATION.'],4,1);
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
if( x(3)>xup )
ierr = 3;
xermsg('SLATEC','AVINT',['THERE WERE LESS THAN THREE FUNCTION VALUES BETWEEN THE ','LIMITS OF INTEGRATION.'],4,1);
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
i = 1;
while(x(i)<xlo);
i = fix(i + 1);
end;
inlft = fix(i);
i = fix(n);
while( x(i)>xup );
i = fix(i - 1);
end;
inrt = fix(i);
if((inrt-inlft)<2 )
ierr = 3;
xermsg('SLATEC','AVINT',['THERE WERE LESS THAN THREE FUNCTION VALUES BETWEEN THE ','LIMITS OF INTEGRATION.'],4,1);
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
istart = fix(inlft);
if( inlft==1 )
istart = 2;
end;
istop = fix(inrt);
if( inrt==n )
istop = fix(n - 1);
end;
%
r3 = 3.0d0;
rp5 = 0.5d0;
summlv = 0.0;
syl = xlo;
syl2 = syl.*syl;
syl3 = syl2.*syl;
%
for i = istart : istop;
x1 = x(i-1);
x2 = x(i);
x3 = x(i+1);
x12 = x1 - x2;
x13 = x1 - x3;
x23 = x2 - x3;
term1 = (y(i-1))./(x12.*x13);
term2 = -(y(i))./(x12.*x23);
term3 = (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 )
ca = a;
cb = b;
cc = c;
else;
ca = 0.5.*(a+ca);
cb = 0.5.*(b+cb);
cc = 0.5.*(c+cc);
end;
syu = x2;
syu2 = syu.*syu;
syu3 = syu2.*syu;
summlv = summlv + ca.*(syu3-syl3)./r3 + cb.*rp5.*(syu2-syl2)+ cc.*(syu-syl);
ca = a;
cb = b;
cc = c;
syl = syu;
syl2 = syu2;
syl3 = syu3;
end; i = fix(istop+1);
syu = xup;
ansmlv = summlv + ca.*(syu.^3-syl3)./r3 + cb.*rp5.*(syu.^2-syl2)+ cc.*(syu-syl);
else;
%
%     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));
ansmlv = 0.5.*(fl+fr).*(xup-xlo);
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
end;
elseif( xlo~=xup ) ;
ierr = 2;
xermsg('SLATEC','AVINT',['THE UPPER LIMIT OF INTEGRATION WAS NOT GREATER THAN THE ','LOWER LIMIT.'],4,1);
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end
%DECK BAKVEC