function [int, tol1,ix]= gaussq(fun,xlow,xhigh,tol,trace,p1,p2,p3,p4,p5,p6,p7,p8,p9)
%GAUSSQ Numerically evaluates a  integral using a Gauss quadrature.
%       The Quadrature integrates a (2m-1)th order  polynomial exactly 
%       and the  integral is of the form
%              b                              
%             Int (p(x)* Fun(x)) dx 
%              a                 
%
% CALL:
% [int, tol] = gaussq('Fun',xlow,xhigh,[reltol wfun],[trace,gn],p1,p2,....)
%
% [int, tol] = gaussq('Fun',xlow,xhigh,[reltol wfun],[trace,gn],alpha,p1,p2,....)
%
%  [int, tol] = gaussq('Fun',xlow,xhigh,[reltol wfun],[trace,gn],alpha,beta,p1,p2,....)
%
%       int = evaluated integral
%       tol = absolute tolerance abs(int-intold)
%       Fun = string containing the name of the function or a directly given 
%              expression enclosed in parenthesis. 
%xlow,xhigh = integration limits
%    reltol = relative tolerance default=1e-3
%      wfun = weight function
%         1  p(x)=1                       a =-1,   b = 1 Legendre (default)
%         2  p(x)=1/sqrt((x-a)*(b-x)),    a =-1,   b = 1 Chebyshev of the
%                                                             first kind
%         3  p(x)=sqrt((x-a)*(b-x)),      a =-1,   b = 1 Chebyshev of the 
%                                                            second kind
%         4  p(x)=sqrt((x-a)/(b-x)),      a = 0,   b = 1
%         5  p(x)=1/sqrt(b-x),            a = 0,   b = 1
%         6  p(x)=sqrt(b-x),              a = 0,   b = 1
%         7  p(x)=(x-a)^alpha*(b-x)^beta  a =-1,   b = 1 Jacobi 
%                                     alpha, beta >-1 (default alpha=beta=0)
%         8  p(x)=x^alpha*exp(-x)         a = 0,   b = inf generalized Laguerre
%         9  p(x)=exp(-x^2)               a =-inf, b = inf Hermite
%        10  p(x)=1                       a =-1,   b = 1 Legendre (slower than 1)
%
%   trace = for non-zero TRACE traces the function evaluations 
%              with a point plot of the integrand.
%      gn = number of base points and weight points to start the 
%            integration with (default=2)
%p1,p2,...= coefficients to be passed directly to function Fun:   
%                  G = Fun(x,p1,p2,...).
%
% Note that int is the common size of xlow, xhigh and p1,p2,....
% Example: a) integration of x.^2 from 0 to 2 and from 2 to 4 is done by:
%                        gaussq('(x.^2)',[0 2],[2 4])
%          b) integration of x^2*exp(-x) is done by 
%                        gaussq('(1)',0,inf,[1e-3 8],[],2)
%             or
%                        gaussq('(x.^2)',0,inf,[1e-3 8],[],0)
% See also qrule



%This function works just like QUAD or QUAD8 but uses a Gauss-Chebyshev
% quadrature integration scheme. The Quadrature integrates a 
%(2m-1)th order  polynomial exactly and the  integral is of the form
%                   b                                m
%                 Int (p(x)* fun(x)) dx  =  Sum ( wf_j* fun( bp_j ) )
%                  a                                j=1                   
% 
% modified version of quadc and quadg by Per A. Brodtkorb 30.03.99, 17.02.99 :
% -accept multiple integrationlimits, int is the common size of xlow and xhigh
% -optimized by only computing the integrals which did not converge.
%  - enabled the integration of directly given functions enclosed in 
%     parenthesis. Example: integration from 0 to 2 and from 2 to 4 for x is done by:
%                        gaussq('(x.^2)',[0 2],[2 4])


% Reference 
%   wfun 1: from grule.m in NIT toolbox, see ref [2] 
%   wfun 2-6: see ref [4]
%   wfun 7-10:  Adapted from Netlib routine gaussq.f see ref [1,3]
%
% [1]  Golub, G. H. and Welsch, J. H. (1969)
% 'Calculation of Gaussian Quadrature Rules'
%  Mathematics of Computation, vol 23,page 221-230,
%
% [2] Davis and Rabinowitz (1975) 'Methods of Numerical Integration', page 365,
%     Academic Press.
%
% [3]. Stroud and Secrest (1966), 'gaussian quadrature formulas', 
%      prentice-hall, Englewood cliffs, n.j.
% 
% [4] Abromowitz and Stegun (1954) 'Handbook of mathematical functions'


global ALPHA1 BETA1 ALPHA2
wfun=1;
if nargin<4| isempty(tol),
  tol=1e-3;
elseif length(tol)==2,
  wfun=tol(2);
   tol=tol(1);
end



istart=1; alpha1=0;beta1=0;
if (wfun==7)|(wfun==8),
  istart=2;
  if(nargin>=6&~isempty(p1)),
    alpha1=p1;
  end
end    
if wfun==7,
  istart=3;
  if nargin>=7&~isempty(p2),
    beta1=p2;
  end
end


FindWeigths=1;
switch wfun
  case 7,
    if isempty(ALPHA1)|isempty(BETA1),
    elseif ALPHA1==alpha1 & BETA1==beta1,
      FindWeigths=0;
    end
    ALPHA1=alpha1;BETA1=beta1;
      %remember what type of weights are saved as global constants
  case 8,
  if isempty(ALPHA2),
    elseif ALPHA2==alpha1,
      FindWeigths=0;
    end
    ALPHA2=alpha1;
    %remember what type of weights are saved as global constants
  otherwise,FindWeigths=0;
end


gn=2;
if nargin <5|isempty(trace) ,
  trace=0; 
elseif length(trace)==2,
  gn=trace(2);
  trace=trace(1);
end



if prod(size(xlow))==1,% make sure the integration limits have correct size
  xlow=xlow(ones(size(xhigh)));;
elseif prod(size(xhigh))==1,
  xhigh=xhigh(ones(size(xlow)));;
elseif any( size(xhigh)~=size(xlow) )
  error('The input must have equal size!')
end
[N M]=size(xlow);%remember the size of input


if any(fun=='('), %  & any(fun=='x'),
  exec_string=['y=',fun ';']; %the call function is already setup
else
  %setup string to call the function
  exec_string=['y=',fun,'(x'];
  num_parameters=nargin-5;
  for i=istart:num_parameters,
    xvar=['p' int2str(i)]; % variable # i
    if eval(['~ischar(' ,xvar,')  & length(',xvar,'(:)) ~=1' ]) ,
        if N*M==1,     
          eval(['[N M]=size(', xvar, ');']); 
        elseif  eval(['N*M~=prod(size(', xvar, '));']);
          error('The input must have equal size!')
        end     
      eval([xvar, '=' ,xvar ,'(:);']); %make sure it is a column 
      exec_string=[exec_string,',' xvar '(k,ones(1,gn) )']; %enable integration with multiple arguments
    else
      exec_string=[exec_string,',' xvar];
    end
  end
  exec_string=[exec_string,');'];
end



nk=N*M;% # of integrals we have to compute
k=(1:nk)';
int=zeros(nk,1);
tol1=int;


%gntxt=int2str(gn);% # of weights
wfuntxt= int2str(wfun);% weightfunction used
wtxt=[int2str(gn),'_',wfuntxt];% weightfunction used
eval(['global cb',wtxt,' cw',wtxt,';']);
if isempty(eval(['cb',wtxt]))|FindWeigths ,  
  % calculate the weights if necessary
  eval(['[cb',wtxt,',cw',wtxt,']=qrule(gn,wfun,alpha1,beta1);']);
end


%setup mapping parameters and execution strings
xlow=xlow(:);
jacob=(xhigh(:)-xlow(:))/2;


switch wfun,% this is clumsy and can written more tidy
  case {1 ,10},
    calcx_string=['(ones(nk,1),:)+1).*jacob(k,ones(1, gn ))+xlow(k,ones(1,gn ));'];
    int_string=['(ones(nk,1),:).*y,2).*jacob(k);'];
  case 2,  
    calcx_string=['(ones(nk,1),:)+1).*jacob(k,ones(1, gn ))+xlow(k,ones(1,gn ));'];
    int_string=['(ones(nk,1),:).*y,2);'];
  case 3, 
    calcx_string=['(ones(nk,1),:)+1).*jacob(k,ones(1, gn ))+xlow(k,ones(1,gn ));'];
    int_string=['(ones(nk,1),:).*y,2).*jacob(k).^2;'];
  case 4,  
    calcx_string=['(ones(nk,1),:)).*jacob(k,ones(1, gn ))*2+xlow(k,ones(1,gn ));'];
    int_string=['(ones(nk,1),:).*y,2).*jacob(k)*2;'];
  case 5,  
    calcx_string=['(ones(nk,1),:)).*jacob(k,ones(1, gn ))*2+xlow(k,ones(1,gn ));'];
    int_string=['(ones(nk,1),:).*y,2).*sqrt(jacob(k)*2);'];
  case 6,  
    calcx_string=['(ones(nk,1),:)).*jacob(k,ones(1, gn ))*2+xlow(k,ones(1,gn ));'];
    int_string=['(ones(nk,1),:).*y,2).*sqrt(jacob(k)*2).^3;'];
  case 7, 
    calcx_string=['(ones(nk,1),:)+1).*jacob(k,ones(1, gn ))+xlow(k,ones(1,gn ));'];
    int_string=['(ones(nk,1),:).*y,2).*jacob(k).^(alpha1+beta1+1);'];
  case {8,9}
     calcx_string=['(ones(nk,1),:));'];
     int_string=['(ones(nk,1),:).*y,2);'];
otherwise error('unknown option')
end



eval(['x=(cb',wtxt, calcx_string]); % calculate the x values            
eval(exec_string); % calculate function values  y=fun(x)                        
eval(['int(k)=sum(cw',wtxt, int_string]);       % do the integration
int_old=int;


if trace==1,
  x_trace=x(:);
  y_trace=y(:);
end


% Break out of the iteration loop for three reasons:
%  1) the last update is very small (compared to int  and  compared to tol)
%  2) There are more than 11 iterations. This should NEVER happen. 


converge='n';
for i=1:10,
  gn=gn*2;% double the # of weights
 % gntxt=int2str(gn);% # of weights
  wtxt=[int2str(gn),'_',wfuntxt];
  eval(['global cb',wtxt,' cw',wtxt]);
  if isempty(eval(['cb',wtxt]))|FindWeigths ,  
    % calculate the weights if necessary
    eval(['[cb',wtxt,',cw',wtxt,']=qrule(gn,wfun,alpha1,beta1);']);
  end
 
  eval(['x=(cb',wtxt, calcx_string]); % calculate the x values  
  eval(exec_string);                 % calculate function values  y=fun(x)