Mercurial > octave
diff src/quad.cc @ 1:78fd87e624cb
[project @ 1993-08-08 01:13:40 by jwe]
Initial revision
author | jwe |
---|---|
date | Sun, 08 Aug 1993 01:13:40 +0000 |
parents | |
children | d68036bcad4c |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/quad.cc Sun Aug 08 01:13:40 1993 +0000 @@ -0,0 +1,187 @@ +// tc-quad.cc -*- C++ -*- +/* + +Copyright (C) 1993 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. + +*/ + +#ifdef __GNUG__ +#pragma implementation +#endif + +#include "Quad.h" + +#include "tree-const.h" +#include "variables.h" +#include "mappers.h" +#include "gripes.h" +#include "error.h" +#include "utils.h" + +// Global pointer for user defined function required by quadrature functions. +static tree *quad_fcn; + +#ifdef WITH_DLD +tree_constant * +builtin_quad_2 (tree_constant *args, int nargin, int nargout) +{ + return do_quad (args, nargin, nargout); +} +#endif + +double +quad_user_function (double x) +{ + double retval = 0.0; + +// tree_constant name = tree_constant (quad_fcn->name ()); + tree_constant *args = new tree_constant [2]; +// args[0] = name; + args[1] = tree_constant (x); + + if (quad_fcn != NULL_TREE) + { + tree_constant *tmp = quad_fcn->eval (args, 2, 1, 0); + delete [] args; + if (tmp != NULL_TREE_CONST && tmp[0].is_defined ()) + { + retval = tmp[0].to_scalar (); + delete [] tmp; + } + else + { + delete [] tmp; + gripe_user_supplied_eval ("quad"); + jump_to_top_level (); + } + } + + return retval; +} + +tree_constant * +do_quad (tree_constant *args, int nargin, int nargout) +{ +// Assumes that we have been given the correct number of arguments. + + tree_constant *retval = NULL_TREE_CONST; + + quad_fcn = is_valid_function (args[1], "fsolve", 1); + if (quad_fcn == NULL_TREE + || takes_correct_nargs (quad_fcn, 2, "fsolve", 1) != 1) + return retval; + + double a = args[2].to_scalar (); + double b = args[3].to_scalar (); + + int indefinite = 0; + IndefQuad::IntegralType indef_type = IndefQuad::doubly_infinite; + double bound = 0.0; + if ((int) xisinf (a) && (int) xisinf (b)) + { + indefinite = 1; + indef_type = IndefQuad::doubly_infinite; + } + else if ((int) xisinf (a)) + { + indefinite = 1; + bound = b; + indef_type = IndefQuad::neg_inf_to_bound; + } + else if ((int) xisinf (b)) + { + indefinite = 1; + bound = a; + indef_type = IndefQuad::bound_to_inf; + } + + int ier = 0; + int nfun = 0; + double abserr = 0.0; + double val = 0.0; + double abstol = 1e-6; + double reltol = 1e-6; + Vector tol (2); + Vector sing; + int have_sing = 0; + switch (nargin) + { + case 6: + if (indefinite) + { + message ("quad", "sorry, singularities not allowed on infinite intervals"); + return retval; + } + have_sing = 1; + sing = args[5].to_vector (); + case 5: + tol = args[4].to_vector (); + switch (tol.capacity ()) + { + case 2: + reltol = tol.elem (1); + case 1: + abstol = tol.elem (0); + break; + default: + message ("quad", "expecting tol to contain no more than two values"); + return retval; + } + case 4: + if (indefinite) + { + IndefQuad iq (quad_user_function, bound, indef_type, abstol, reltol); + val = iq.integrate (ier, nfun, abserr); + } + else + { + if (have_sing) + { + DefQuad dq (quad_user_function, a, b, sing, abstol, reltol); + val = dq.integrate (ier, nfun, abserr); + } + else + { + DefQuad dq (quad_user_function, a, b, abstol, reltol); + val = dq.integrate (ier, nfun, abserr); + } + } + break; + default: + panic_impossible (); + break; + } + + retval = new tree_constant [5]; + + retval[0] = tree_constant (val); + retval[1] = tree_constant ((double) ier); + retval[2] = tree_constant ((double) nfun); + retval[3] = tree_constant (abserr); + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; page-delimiter: "^/\\*" *** +;;; End: *** +*/ +