## -*- texinfo -*-
## @deftypefn  {} {[@var{q}, @var{num_iter}] =} quadgF (@var{f}, @var{g}, @var{lim_I}, @var{lim_S})
## @deftypefnx  {} {[@var{q}, @var{num_iter}] =} quadgF (@var{f}, @var{g}, @var{lim_I}, @var{lim_S}, @var{tol})
## @deftypefnx  {} {[@var{q}, @var{num_iter}] =} quadgF (@var{f}, @var{g}, @var{lim_I}, @var{lim_S}, @var{tol}, @var{max_iter})
##
## This function evaluates numerically the integral:
## @verbatim    
##        ⌠LIM_S                                    
##        │     F(x) exp[iG(x)] dx               
##        ⌡LIM_I                                    
## @end verbatim
## using an adaptative quadrature described in reference [1]. Where
## @var{f} and @var{g} are function handles and i = sqrt(-1). Both functions 
## @var{f} and @var{g} must be vectorized and return a vector of output values 
## when given a vector of input values.
##
## @var{lim_I} and @var{lim_S} are the lower and upper limits of integration.  
## Both limits must be finite.
##
## The optional argument @var{tol} defines the absolute tolerance with which
## to perform the integration. The default value is 1e-6.
##
## The optional argument @var{max_iter} defines the maximum number of iterations
## the principal loop performs. The default is 5. The user can set a higher 
## number, but a very large number of @var{num_iter} means the algoritm used in 
## this function is not suitable for calculating the integral.
##
## The result of the integration is returned in @var{q}.
##
## The optional output @var{num_iter} indicates the number of iterations the 
## principal loop performs.
##
## Example
## @example
## [q,num_iter] = quadgF(@@(x) sin(2π*x), @@(x) 1000*x, 0, 1);
## @end example
##
## This function, @code{quadgF}, implements the algorithm described in reference [1]. 
## However, instead of using the Fresnel integrals, as in the article,
## the octave's built-in function @strong{erf} is used.
##
## [1] Shampine, L.F. Integrating oscillatory functions in MATLAB, II. Electronic
##     Transactions on Numerical Analysis, volume 39, pp 403-413, 2012.
##


function [q,num_iter] = quadgF(f, g, lim_I, lim_S, tol = 1e-6, max_iter = 5);

  if (nargin < 4)
    print_usage ();
  endif
  if (! isscalar (tol) || tol < 0)
    error ("quadgF: TOL must be a scalar >=0");
  elseif (tol < eps)
    warning ("quadgF: TOL specified is smaller than machine precision, using %g",tol);
    tol = eps;
  endif
  if (lim_S < lim_I)
    ## Reverse integration
    [q, ni] = quadgF (f, g, lim_S, lim_I, tol, max_iter);
    q = -q;
    num_iter = ni;
    return;
  endif

  dom_int = lim_S-lim_I;				# Integration interval width
  Num_intervals = 32;
  q = 0;
  num_iter = 0;
  tol2 = tol^2;
  lint = dom_int/Num_intervals;		# 1st pass sub-interval width

# The values of functions F and G are calculated at 129 equi-spaced points 
# in the integration interval
  x = linspace(lim_I,lim_S,129);
  fm = f(x); 
  gm = g(x); 

# The previous values are arranged in a three-dimensional array, so that each 
# page contains the data of one of the 32 initial intervals
  C = zeros(3,5,Num_intervals);
  for m = 1:Num_intervals
      n = (m-1)*4+1;
      C(:,:,m) = [x(n:n+4); fm(n:n+4); gm(n:n+4)];
  endfor

# The size of F and G in the interval [limI, limS] are calculated using
# the Simpson's rule with 129 equi-spaced points
  fm2 = fm.^2;   gm2 = gm.^2;
  F2 = (lint/3/4)*(fm2(1)+2*sum(fm2(3:2:end-2))+4*sum(fm2(2:2:end))+fm2(end));
  G2 = (lint/12)*(gm(1)+2*sum(gm(3:2:end-2))+4*sum(gm(2:2:end))+gm(end));

# Functions F and G are approximated by splines
# G ≈ A1 + A2 * t + A3 * t²       with t = x-xm
# F ≈ B1 + B2 * t + B3 * t²       with t = x-xm
  do
    ++num_iter;
    lint2 = lint^2; lint4 = lint/4; lint42 = lint4^2;
    B1 = zeros(1,Num_intervals);
    B2 = B1;  B3 = B1;  A1 = B1;  A2 = B1;  A3 = B1; 
    Hf1 = B1; Hf2 = B1; Hg1 = B1; Hg2 = B1;
    B1(:) = C(2,1,:);                                    # Coefficients B1,B2, 
    B2(:) = (4*C(2,3,:) - 3*C(2,1,:) - C(2,5,:))/lint;   # and B3 are calculated
    B3(:) = (2*C(2,1,:) - 4*C(2,3,:) + 2*C(2,5,:))/lint2;
    Hf1(:) = C(2,2,:); Hf1 = (Hf1 - (B1 + B2 * lint4 + B3 * lint42)).^2;
    Hf2(:) = C(2,4,:); Hf2 = (Hf2 - (B1 + B2 * 3*lint4 + B3 * 9*lint42)).^2;
    Hf = (256/945)*lint*(Hf1+Hf2);
    A1(:) = C(3,1,:);                                    # Coefficients A1,A2, 
    A2(:) = (4*C(3,3,:) - 3*C(3,1,:) - C(3,5,:))/lint;   # and A3 are calculated
    A3(:) = (2*C(3,1,:) - 4*C(3,3,:) + 2*C(3,5,:))/lint2;
    Hg1(:) = C(3,2,:); Hg1 = (Hg1 - (A1 + A2 * lint4 + A3 * lint42)).^2;
    Hg2(:) = C(3,4,:); Hg2 = (Hg2 - (A1 + A2 * 3*lint4 + A3 * 9*lint42)).^2;
    Hg = (256/945)*lint*(Hg1+Hg2);
    condF = (Hf <= tol2*F2*(lint/dom_int));    # Check whether the approach 
    condG = (Hg <= tol2*G2*(lint/dom_int));    # by splines at the different 
    condFG = condF.*condG;                     # intervals is good enough 
    pass = find(condFG);                                       
    nopass = find(~condFG);                                   
    Num_intervals = numel(nopass);                            
    if and(num_iter == max_iter, Num_intervals > 0)
       q = q + int_pass(A1,A2,A3,B1,B2,B3,lint,lint2,tol);
       warning("Number of iterations = max_iter. The result error can be greater than Tolerance");
       Num_intervals = 0;
    else
       q = q + int_pass(A1(pass),A2(pass),A3(pass),B1(pass),B2(pass),B3(pass),lint,lint2,tol);
       C1 = C(:,:,nopass);
       lint = lint/2;
       if Num_intervals > 0                    # Intervals in which the 
          NC = zeros(3,5,2*Num_intervals);     # approach by splines wasn't 
          for m = 1:Num_intervals              # good enough are divided in half
             n = 2*m-1; 
             NC(:,1:2:5,n) = C1(:,1:1:3,m);          NC(:,1:2:5,n+1) = C1(:,3:1:5,m);
             NC(1,2:2:4,n) = C1(1,1:1:2,m)+lint/4;   NC(1,2:2:4,n+1) = C1(1,3:1:4,m)+lint/4;
             NC(2,2:2:4,n) = f(NC(1,2:2:4,n));       NC(2,2:2:4,n+1) = f(NC(1,2:2:4,n+1));
             NC(3,2:2:4,n) = g(NC(1,2:2:4,n));       NC(3,2:2:4,n+1) = g(NC(1,2:2:4,n+1));
          endfor
          C = NC;
          Num_intervals = 2*Num_intervals;
       endif
    endif
  until Num_intervals == 0

endfunction

function out_run = int_pass(a1,a2,a3,b1,b2,b3,lint,lint2,tol);
  norm_st2 = a1.^2*lint+a1.*a2*lint2+lint^3*(2*a1.*a3+a2.^2)/3+a2.*a3*lint2^2/2+(a3.^2*lint^5)/5;
  norm_st = sqrt(norm_st2);
  norm_stm_a1 = sqrt((a2.^2*lint^3)/3+(a2.*a3*lint2^2)/2+(a3.^2*lint^5)/5);
  norm_stm_a1a2 = sqrt((a3.^2*lint^5)/5);
  const_app = norm_stm_a1<=tol*norm_st;
  line_app = (norm_stm_a1a2<=tol*norm_st)-const_app;
  quad_app = ~(const_app + line_app);
  const_app = find(const_app);
  line_app = find(line_app);
  quad_app = find(quad_app);
  pro_int = [0, 0, 0];
  if numel(const_app)>0
     pro_int(1) = int_cte(a1(const_app),b1(const_app),b2(const_app),b3(const_app),lint);
  endif
  if numel(line_app)>0
     pro_int(2) = int_line(a1(line_app),a2(line_app),b1(line_app),b2(line_app),b3(line_app),lint);
  endif
  if numel(quad_app)>0   
     pro_int(3) = int_quad(a1(quad_app),a2(quad_app),a3(quad_app),b1(quad_app),b2(quad_app),b3(quad_app),lint);
  endif   
  out_run = sum(pro_int);
endfunction

function out1 = int_cte(a1,b1,b2,b3,il);               
   int_int = exp(i*a1) .* ((b1.*il) + (b2.*il.^2)/2 + (b3.*il.^3)/3);
   out1 = sum(int_int);
endfunction

function out2 = int_line(a1,a2,b1,b2,b3,il);          
   mu = exp(i*a2.*il);
   I1 = (mu-1)./(i*a2);
   I2 = (il.*mu-I1)./(i*a2);
   I3 = (il.^2 .* mu - 2*I2)./(i*a2);
   int_int = exp(i*a1) .* (b1.*I1 + b2.*I2 + b3.*I3);
   out2 = sum(int_int);
endfunction

function out3 = int_quad(a1,a2,a3,b1,b2,b3,il);        
   c1 = a1-(a2.^2./(4*a3));     c2 = a2./(2*a3);       c3 = a3;
   d1 = b1-b2.*c2+b3.*c2.^2;    d2 = b2-2*c2.*b3;      d3 = b3;
   u1 = il + c2;                u0 = c2;
   v1 = erf(-sqrt(-i*c3).*u1);  v0 = erf(-sqrt(-i*c3).*u0);
   w1 = exp(i*c3.*u1.^2);       w0 = exp(i*c3.*u0.^2);
   intd1 = -0.5 * sqrt(i*pi./c3) .* (v1-v0);
   intd2 = (-i./(2*c3)) .* (w1-w0);
   intd3 = (-i./(4*c3.^(3/2))) .* (sqrt(i*pi)*(v1-v0)+2*sqrt(c3).*(u1.*w1-u0.*w0));
   int_int = exp(i*c1) .* (d1.*intd1 + d2.*intd2 + d3.*intd3);
   out3 = sum(int_int);
endfunction