% QUADR Numerical evaluation of an integral.
%	Q = QUADS(F,X1,X2) approximates the integral of F(X)
%	from X1 to X2 with relative error  1e-4.  'F' is a string
%	containing the name of the function such as 'sin' or expression
%	depending on x (such as '1/(x+1)^2'). Function or expression F
%	must return a vector of output values of the same size as a
%	vector of input values.
%
%	Q = QUADS(F,X1,X2,TOL) integrates to a relative error of TOL.
%	Q = QUADS(F,X1,X2,TOL,TRACE) also traces the function
%	evaluations with a point plot (with different colors for each
%	recursion level).
%	Q = QUADS(F,X1,X2,TOL,TRACE,P1,P2,...) allows parameters P1, P2
%	to be passed directly to a function F in the form F(X,P1,P2,...).
%	When F stands for expression, its parameter dependence must
%	contain 'p1' ... explicitly, such as '(x-p1)/(x+p2)^p3'.
%	[Q,X,Y] = QUADS(F,...) also returns vector of points X of
%	function evaluation and corresponding function values Y.

%  Kirill K. Pankratov,  kirill@plume.mit.edu
%  01/30/95

require lininsrt

quadr = function (fun, x1, x2, tol, trace, p1, p2, p3, p4, p5)
{
  local (x1, x2, tol, trace, p1, p2, p3, p4, p5)

  % Defaults and parameters ..................................
  tol_dflt = 1e-4; % Default tolerance for integral estimation
  trace_dflt = 0;  % No plotting by default
  n_max = 7;       % Max. number of points inserted at one
                   % iteration between 2 old points
  max_lev = 8;     % Max. number of iterations (recursion level)
  pow = 1/4;       % Power index for estimation of number of
                   % inserted points

  % Coefficients ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  % Function interpolation coefficients ..........
  % (Lagrange polynomials interp. from 4 basis pts to 3 pts.
  ci =  [-1/6, 2/3,  0,  2/3, -1/6];   % [1 2 4 5] to 3

  % Integral approximation coefficients ..........
  cq = [0.31111111111111;
        1.42222222222222;
        0.53333333333334;
        1.42222222222222;
        0.31111111111111];


  % Handle input .............................................
  if (nargs < 5) { trace = trace_dflt; }
  if (nargs < 4) { tol = tol_dflt; }
  if (trace == []) { trace = trace_dflt; }
  if (tol == [] || all([tol[:]; -1]) <= 0) { tol = tol_dflt; }

  % Set up function call .....................................
  if (class (fun) == "function")
  {
    call = "fun ( x";
    for (n in 1:nargs-5)
    {
      call = call + ", p" + num2str(n);
    }
    call = "f = " + call + ");";

  else               % .. Expression

    call = "f = " + fun + ";";
  }

  % Initial layout ..............................
  g = sqrt(3)/3;  % Initial layout section ratio (just some
                  % number close but not equal to 1/2)
  span = x2 - x1;
  x_o = [x1, x1+[1-g, 1+g/117, 1+g]*span/2, x2];

  % Auxillary ...................................
  ci5 = ci[5];
  ci = ci[1:4];
  cq5 = cq[5];
  cq = cq[1:4]';
  tol_f = tol/sqrt(span)*16;
  o4 = ones(4,1);

  % Initialize output and control parameters ....
  f_o = [];
  n_ins = 3;
  n_i = 1;
  last = 1;
  iter = 0;
  is_more = 1;

  % Begin recursively refine approximations ..................
  while (is_more)  % While not enough accuracy ``````````````````0
  {
    iter = iter + 1;

    % Insert new points ..........................
    F = n_ins[o4;];
    </ mask ; x_o /> = lininsrt(x_o, F);
    if (iter == 1) { mask = ones(size(mask)); }

    % Expand "last" vector
    last = last[o4;1:n_i];
    last = [last[:]', last[n_i]];

    % Evaluate function at new points ............
    x = x_o[find (mask)];
    eval(call);

    % Insert new points in the output and control vectors
    n_ins = find(!mask);
    misfit = find(mask);
    if (iter > 1)
    {
      f_o[n_ins] = f_o;
      last[n_ins] = last;
    }
    f_o [misfit] = f;
    last[misfit] = ones(size(misfit));

    % Reshape function values vector ............
    l_o = length(x_o);
    n_i = (l_o - 1)/4;
    F = reshape(f_o[1:l_o-1], 4, n_i);

    % Calculate difference between actual 
    % and interpolated values ...................
    misfit = ci*F+ci5*f_o[5:l_o:4];  % Interpolation
    misfit = abs(misfit-F[3;]);      % Difference

    % Calculate how many points need to be inserted
    scale = mean(abs(f_o));
    n_ins = (misfit/(tol_f*scale)).^pow-1;

    % Smooth (moving average)
    if (iter > 2)
    {
      misfit = cumsum(n_ins);
      n_ins[3:n_i-2] = (misfit[5:n_i] - misfit[1:n_i-4])/4;
      n_ins[[1, 2, n_i-1, n_i]] = n_ins[3, 3, n_i-2, n_i-2];
    }
    
    % Make it integer and upper bounded
    n_ins = min(ceil(n_ins), n_max*ones(size(n_ins)));

    % If previous iteration was accurate
    % enough, do not insert new points .........
    misfit = n_ins > 1;
    n_ins = max(n_ins, last[5:l_o:4]);
    last = misfit;

    is_more = any(n_ins);   % If needs more points inserted

    if (trace)    % Plot new points ..............
    {
      plot([ x, f ]);
    }

    % If still not enough accuracy .............
    if (iter > max_lev)
    {
      disp(" Recursion level limit reached. Singularity likely.");
      is_more = 0;
    }

  }  
  % End while (recursion) ''''''''''''''''''''''''''''''''0

  if (trace) { "hold off" }

  % Now calculate the integral itself .........
  last = 5:l_o:4;
  misfit = cq*F+cq5*f_o[last];
  Q = (x_o[last] - x_o[last - 1])*misfit';

  return << Q=Q; x_o=x_o; f_o=f_o >>;
};