function result=trplquad(intfcn,inmin,inmax,midmin,midmax,outmin,outmax,tol,method) 
%TRPLQUAD Numerically evaluate triple integral. 
%   RESULT = TRPLQUAD('FUN',INMIN,INMAX,MIDMIN<MIDMAX<OUTMIN,OUTMAX) evaluates the
%   double integral FUN(INNER,MIDDLE,OUTER) using the QUAD quadrature function,
%   where INNER is the inner variable ranging from INMIN to INMAX, MIDDLE
%is the middle variable, ranging from MIDMIN to MIDMAX, and
%   OUTER is the outer variable ranging from OUTMIN to OUTMAX. The first
%   argument 'FUN' is a string, or inline object, representing the integrand
%   function.  This function must be a function of three variables of the form
%   FOUT=FUN(INNER,MIDDLE,OUTER).  The function must take a vector INNER and a
%   scalars MIDDLE and OUTER and return a vector FOUT that is the function evaluated
%   at MIDDLE and OUTER and each value of INNER.
%
%   RESULT = TRPLQUAD('FUN',INMIN,INMAX,MIDMIN,MIDMAX,OUTMIN,OUTMAX,TOL) passes
%   TOL to the QUAD function. See the help entry for QUAD for a
%   description of the parameter TOL.
%
%   RESULT = TRPLQUAD('FUN',INMIN,INMAX,MIDMIN,MIDMAX,OUTMIN,OUTMAX,TOL,METHOD)
%   passes TOL to the QUAD or QUAD8 function depending on the value
%   of the string METHOD.  Valid values for METHOD are 'quad' and
%   'quad8', or the name of any user-defined quadrature method with the
%   same calling and return arguments as QUAD or QUAD8.
%
%   %
%   See also QUAD, QUAD8, INLINE. 

%   Adapted from dblquad(Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 1.7 $  $Date: 1997/11/21 23:30:41 $) by Paul Green

% Determine input arguments.

if nargin < 7 
  error('Requires at least seven inputs'); 
elseif ~isstr(intfcn) & ~isa(intfcn,'inline')
  error(sprintf(['First argument must be a string, or inline object,'...
        ' representing\n the integrand function.'])); 
end

if nargin < 9, method = 'quad'; end
if nargin < 8, tol = []; end 

% Pass the necessary information to MIDLOOP. MIDLOOP 
% evaluates the middle loop, calling INNERLOOP. Note that we pass TOL 
% twice; once to QUAD (QUAD8) directly, and once indirectly to 
% MIDLOOP so that MIDLOOP can again pass it to INNERLOOP.  
% directly. 
trace = [];
result= feval(method,'midloop', outmin, outmax, tol, trace, intfcn, ...
   inmin, inmax,midmin,midmax, tol, method); 

