Mercurial > forge
changeset 2716:909915aadbbc octave-forge
Commit derived file gsl_sf.cc to the CVS so that the function help can find it
| author | adb014 |
|---|---|
| date | Tue, 17 Oct 2006 19:08:07 +0000 |
| parents | def1f0152603 |
| children | dc09a0aff019 |
| files | main/gsl/src/gsl_sf.cc |
| diffstat | 1 files changed, 9464 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/gsl/src/gsl_sf.cc Tue Oct 17 19:08:07 2006 +0000 @@ -0,0 +1,9464 @@ +// *** DO NOT EDIT *** this file is generated by buildgsl_sf.sh + +/* + Copyright (C) 2004 Teemu Ikonen <tpikonen@pcu.helsinki.fi> + + This program 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 of the License, or + (at your option) any later version. + + This program 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 this program; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + +#include <octave/oct.h> +#include <gsl/gsl_sf.h> +#include <gsl/gsl_sf_result.h> +#include <gsl/gsl_errno.h> + +void octave_gsl_errorhandler (const char * reason, const char * file, + int line, int gsl_errno) +{ + error("GSL error %d at %s, line %d: %s\n", gsl_errno, file, line, reason); +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(clausen, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} clausen (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} clausen (@dots{})\n\ +\n\ +The Clausen function is defined by the following integral,\n\ +\n\ +Cl_2(x) = - \\int_0^x dt \\log(2 \\sin(t/2))\n\ +\n\ +It is related to the dilogarithm by Cl_2(\\theta) = \\Im Li_2(\\exp(i \\theta)).\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_clausen (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_clausen_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(dawson, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} dawson (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} dawson (@dots{})\n\ +\n\ +The Dawson integral is defined by \\exp(-x^2) \\int_0^x dt \\exp(t^2).\n\ +A table of Dawson's integral can be found in Abramowitz & Stegun, Table 7.5.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_dawson (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_dawson_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(debye_1, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} debye_1 (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} debye_1 (@dots{})\n\ +\n\ +The Debye functions are defined by the integral \n\ +\n\ +D_n(x) = n/x^n \\int_0^x dt (t^n/(e^t - 1)).\n\ +\n\ +For further information see Abramowitz & Stegun, Section 27.1.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_debye_1 (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_debye_1_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(debye_2, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} debye_2 (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} debye_2 (@dots{})\n\ +\n\ +The Debye functions are defined by the integral\n\ +\n\ +D_n(x) = n/x^n \\int_0^x dt (t^n/(e^t - 1)).\n\ +\n\ +For further information see Abramowitz & Stegun, Section 27.1.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_debye_2 (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_debye_2_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(debye_3, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} debye_3 (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} debye_3 (@dots{})\n\ +\n\ +The Debye functions are defined by the integral\n\ +\n\ +D_n(x) = n/x^n \\int_0^x dt (t^n/(e^t - 1)).\n\ +\n\ +For further information see Abramowitz & Stegun, Section 27.1.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_debye_3 (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_debye_3_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(debye_4, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} debye_4 (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} debye_4 (@dots{})\n\ +\n\ +The Debye functions are defined by the integral\n\ +\n\ +D_n(x) = n/x^n \\int_0^x dt (t^n/(e^t - 1)).\n\ +\n\ +For further information see Abramowitz & Stegun, Section 27.1.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_debye_4 (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_debye_4_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(erf_gsl, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} erf_gsl (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} erf_gsl (@dots{})\n\ +\n\ +These routines compute the error function \n\ +erf(x) = (2/\\sqrt(\\pi)) \\int_0^x dt \\exp(-t^2).\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_erf (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_erf_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(erfc_gsl, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} erfc_gsl (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} erfc_gsl (@dots{})\n\ +\n\ +These routines compute the complementary error function \n\ +erfc(x) = 1 - erf(x) = (2/\\sqrt(\\pi)) \\int_x^\\infty \\exp(-t^2).\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_erfc (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_erfc_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(log_erfc, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} log_erfc (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} log_erfc (@dots{})\n\ +\n\ +These routines compute the logarithm of the complementary error\n\ +function \\log(\\erfc(x)).\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_log_erfc (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_log_erfc_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(erf_Z, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} erf_Z (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} erf_Z (@dots{})\n\ +\n\ +These routines compute the Gaussian probability function \n\ +Z(x) = (1/(2\\pi)) \\exp(-x^2/2).\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_erf_Z (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_erf_Z_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(erf_Q, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} erf_Q (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} erf_Q (@dots{})\n\ +\n\ +These routines compute the upper tail of the Gaussian probability\n\ +function Q(x) = (1/(2\\pi)) \\int_x^\\infty dt \\exp(-t^2/2).\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_erf_Q (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_erf_Q_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(hazard, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} hazard (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} hazard (@dots{})\n\ +\n\ +The hazard function for the normal distrbution, also known as the \n\ +inverse Mill's ratio, is defined as \n\ +h(x) = Z(x)/Q(x) = \\sqrt@{2/\\pi \\exp(-x^2 / 2) / \\erfc(x/\\sqrt 2)@}. \n\ +It decreases rapidly as x approaches -\\infty and asymptotes to \n\ +h(x) \\sim x as x approaches +\\infty.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_hazard (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_hazard_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(expm1, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} expm1 (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} expm1 (@dots{})\n\ +\n\ +These routines compute the quantity \\exp(x)-1 using an algorithm that \n\ +is accurate for small x.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_expm1 (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_expm1_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(exprel, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} exprel (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} exprel (@dots{})\n\ +\n\ +These routines compute the quantity (\\exp(x)-1)/x using an algorithm \n\ +that is accurate for small x. For small x the algorithm is based on \n\ +the expansion (\\exp(x)-1)/x = 1 + x/2 + x^2/(2*3) + x^3/(2*3*4) + \\dots.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_exprel (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_exprel_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(exprel_2, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} exprel_2 (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} exprel_2 (@dots{})\n\ +\n\ +These routines compute the quantity 2(\\exp(x)-1-x)/x^2 using an\n\ +algorithm that is accurate for small x. For small x the algorithm is \n\ +based on the expansion \n\ +2(\\exp(x)-1-x)/x^2 = 1 + x/3 + x^2/(3*4) + x^3/(3*4*5) + \\dots.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_exprel_2 (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_exprel_2_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(expint_E1, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} expint_E1 (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} expint_E1 (@dots{})\n\ +\n\ +These routines compute the exponential integral E_1(x),\n\ +\n\ +E_1(x) := Re \\int_1^\\infty dt \\exp(-xt)/t.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_expint_E1 (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_expint_E1_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(expint_E2, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} expint_E2 (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} expint_E2 (@dots{})\n\ +\n\ +These routines compute the second-order exponential integral E_2(x),\n\ +\n\ +E_2(x) := \\Re \\int_1^\\infty dt \\exp(-xt)/t^2.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_expint_E2 (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_expint_E2_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(expint_Ei, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} expint_Ei (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} expint_Ei (@dots{})\n\ +\n\ +These routines compute the exponential integral E_i(x),\n\ +\n\ +Ei(x) := - PV(\\int_@{-x@}^\\infty dt \\exp(-t)/t)\n\ +\n\ +where PV denotes the principal value of the integral.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_expint_Ei (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_expint_Ei_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(Shi, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} Shi (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} Shi (@dots{})\n\ +\n\ +These routines compute the integral Shi(x) = \\int_0^x dt \\sinh(t)/t.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_Shi (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_Shi_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(Chi, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} Chi (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} Chi (@dots{})\n\ +\n\ +These routines compute the integral \n\ +\n\ +Chi(x) := Re[ \\gamma_E + \\log(x) + \\int_0^x dt (\\cosh[t]-1)/t] , \n\ +\n\ +where \\gamma_E is the Euler constant.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_Chi (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_Chi_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(expint_3, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} expint_3 (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} expint_3 (@dots{})\n\ +\n\ +These routines compute the exponential integral \n\ +Ei_3(x) = \\int_0^x dt \\exp(-t^3) for x >= 0.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_expint_3 (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_expint_3_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(Si, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} Si (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} Si (@dots{})\n\ +\n\ +These routines compute the Sine integral Si(x) = \\int_0^x dt \\sin(t)/t.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_Si (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_Si_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(Ci, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} Ci (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} Ci (@dots{})\n\ +\n\ +These routines compute the Cosine integral \n\ +Ci(x) = -\\int_x^\\infty dt \\cos(t)/t for x > 0. \n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_Ci (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_Ci_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(atanint, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} atanint (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} atanint (@dots{})\n\ +\n\ +These routines compute the Arctangent integral \n\ +AtanInt(x) = \\int_0^x dt \\arctan(t)/t. \n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_atanint (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_atanint_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(fermi_dirac_mhalf, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} fermi_dirac_mhalf (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} fermi_dirac_mhalf (@dots{})\n\ +\n\ +These routines compute the complete Fermi-Dirac integral F_@{-1/2@}(x).\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_fermi_dirac_mhalf (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_fermi_dirac_mhalf_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(fermi_dirac_half, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} fermi_dirac_half (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} fermi_dirac_half (@dots{})\n\ +\n\ +These routines compute the complete Fermi-Dirac integral F_@{1/2@}(x).\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_fermi_dirac_half (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_fermi_dirac_half_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(fermi_dirac_3half, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} fermi_dirac_3half (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} fermi_dirac_3half (@dots{})\n\ +\n\ +These routines compute the complete Fermi-Dirac integral F_@{3/2@}(x).\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_fermi_dirac_3half (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_fermi_dirac_3half_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(gamma_gsl, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} gamma_gsl (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} gamma_gsl (@dots{})\n\ +\n\ +These routines compute the Gamma function \\Gamma(x), subject to x not \n\ +being a negative integer. The function is computed using the real \n\ +Lanczos method. The maximum value of x such that \\Gamma(x) is not \n\ +considered an overflow is given by the macro GSL_SF_GAMMA_XMAX and is 171.0.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_gamma (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_gamma_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(lngamma_gsl, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} lngamma_gsl (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} lngamma_gsl (@dots{})\n\ +\n\ +These routines compute the logarithm of the Gamma function, \n\ +\\log(\\Gamma(x)), subject to x not a being negative integer. \n\ +For x<0 the real part of \\log(\\Gamma(x)) is returned, which is \n\ +equivalent to \\log(|\\Gamma(x)|). The function is computed using \n\ +the real Lanczos method.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_lngamma (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_lngamma_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(gammastar, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} gammastar (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} gammastar (@dots{})\n\ +\n\ +These routines compute the regulated Gamma Function \\Gamma^*(x) \n\ +for x > 0. The regulated gamma function is given by,\n\ +\n\ +\\Gamma^*(x) = \\Gamma(x)/(\\sqrt@{2\\pi@} x^@{(x-1/2)@} \\exp(-x))\n\ + = (1 + (1/12x) + ...) for x \\to \\infty\n\ +\n\ +and is a useful suggestion of Temme. \n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_gammastar (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_gammastar_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(gammainv_gsl, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} gammainv_gsl (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} gammainv_gsl (@dots{})\n\ +\n\ +These routines compute the reciprocal of the gamma function, 1/\\Gamma(x) using the real Lanczos method.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_gammainv (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_gammainv_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(lambert_W0, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} lambert_W0 (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} lambert_W0 (@dots{})\n\ +\n\ +These compute the principal branch of the Lambert W function, W_0(x).\n\ +\n\ +Lambert's W functions, W(x), are defined to be solutions of the\n\ +equation W(x) \\exp(W(x)) = x. This function has multiple branches \n\ +for x < 0; however, it has only two real-valued branches. \n\ +We define W_0(x) to be the principal branch, where W > -1 for x < 0, \n\ +and W_@{-1@}(x) to be the other real branch, where W < -1 for x < 0.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_lambert_W0 (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_lambert_W0_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(lambert_Wm1, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} lambert_Wm1 (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} lambert_Wm1 (@dots{})\n\ +\n\ +These compute the secondary real-valued branch of the Lambert \n\ +W function, W_@{-1@}(x).\n\ +\n\ +Lambert's W functions, W(x), are defined to be solutions of the\n\ +equation W(x) \\exp(W(x)) = x. This function has multiple branches \n\ +for x < 0; however, it has only two real-valued branches. \n\ +We define W_0(x) to be the principal branch, where W > -1 for x < 0, \n\ +and W_@{-1@}(x) to be the other real branch, where W < -1 for x < 0.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_lambert_Wm1 (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_lambert_Wm1_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(log_1plusx, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} log_1plusx (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} log_1plusx (@dots{})\n\ +\n\ +These routines compute \\log(1 + x) for x > -1 using an algorithm that\n\ +is accurate for small x.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_log_1plusx (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_log_1plusx_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(log_1plusx_mx, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} log_1plusx_mx (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} log_1plusx_mx (@dots{})\n\ +\n\ +These routines compute \\log(1 + x) - x for x > -1 using an algorithm \n\ +that is accurate for small x.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_log_1plusx_mx (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_log_1plusx_mx_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(psi, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} psi (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} psi (@dots{})\n\ +\n\ +These routines compute the digamma function \\psi(x) for general x, \n\ +x \\ne 0.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_psi (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_psi_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(psi_1piy, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} psi_1piy (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} psi_1piy (@dots{})\n\ +\n\ +These routines compute the real part of the digamma function on \n\ +the line 1+i y, Re[\\psi(1 + i y)].\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_psi_1piy (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_psi_1piy_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(synchrotron_1, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} synchrotron_1 (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} synchrotron_1 (@dots{})\n\ +\n\ +These routines compute the first synchrotron function \n\ +x \\int_x^\\infty dt K_@{5/3@}(t) for x >= 0.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_synchrotron_1 (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_synchrotron_1_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(synchrotron_2, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} synchrotron_2 (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} synchrotron_2 (@dots{})\n\ +\n\ +These routines compute the second synchrotron function \n\ +x K_@{2/3@}(x) for x >= 0.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_synchrotron_2 (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_synchrotron_2_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(transport_2, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} transport_2 (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} transport_2 (@dots{})\n\ +\n\ +These routines compute the transport function J(2,x).\n\ +\n\ +The transport functions J(n,x) are defined by the integral\n\ +representations J(n,x) := \\int_0^x dt t^n e^t /(e^t - 1)^2.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_transport_2 (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_transport_2_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(transport_3, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} transport_3 (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} transport_3 (@dots{})\n\ +\n\ +These routines compute the transport function J(3,x).\n\ +\n\ +The transport functions J(n,x) are defined by the integral\n\ +representations J(n,x) := \\int_0^x dt t^n e^t /(e^t - 1)^2.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_transport_3 (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_transport_3_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(transport_4, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} transport_4 (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} transport_4 (@dots{})\n\ +\n\ +These routines compute the transport function J(4,x).\n\ +\n\ +The transport functions J(n,x) are defined by the integral\n\ +representations J(n,x) := \\int_0^x dt t^n e^t /(e^t - 1)^2.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_transport_4 (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_transport_4_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(transport_5, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} transport_5 (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} transport_5 (@dots{})\n\ +\n\ +These routines compute the transport function J(5,x).\n\ +\n\ +The transport functions J(n,x) are defined by the integral\n\ +representations J(n,x) := \\int_0^x dt t^n e^t /(e^t - 1)^2.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_transport_5 (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_transport_5_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(sinc_gsl, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} sinc_gsl (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} sinc_gsl (@dots{})\n\ +\n\ +These routines compute \\sinc(x) = \\sin(\\pi x) / (\\pi x) for any value of x.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_sinc (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_sinc_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(lnsinh, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} lnsinh (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} lnsinh (@dots{})\n\ +\n\ +These routines compute \\log(\\sinh(x)) for x > 0.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_lnsinh (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_lnsinh_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(lncosh, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} lncosh (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} lncosh (@dots{})\n\ +\n\ +These routines compute \\log(\\cosh(x)) for any x.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_lncosh (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_lncosh_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(zeta, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} zeta (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} zeta (@dots{})\n\ +\n\ +These routines compute the Riemann zeta function \\zeta(s) for \n\ +arbitrary s, s \\ne 1.\n\ +\n\ +The Riemann zeta function is defined by the infinite sum \n\ +\\zeta(s) = \\sum_@{k=1@}^\\infty k^@{-s@}. \n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_zeta (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_zeta_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(eta, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} eta (@var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} eta (@dots{})\n\ +\n\ +These routines compute the eta function \\eta(s) for arbitrary s.\n\ +\n\ +The eta function is defined by \\eta(s) = (1-2^@{1-s@}) \\zeta(s).\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_eta (x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_eta_e (x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +DEFUN_DLD(bessel_Jn, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} bessel_Jn (@var{n}, @var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} bessel_Jn (@dots{})\n\ +\n\ +These routines compute the regular cylindrical Bessel function of\n\ +order n, J_n(x).\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_Jn (static_cast<int>(n.xelem(i)), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_Jn_e (static_cast<int>(n.xelem(i)), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + int nint = static_cast<int>(n.xelem(0)); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_Jn (nint, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_Jn_e (nint, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_Jn (static_cast<int>(n.xelem(i)), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_Jn_e (static_cast<int>(n.xelem(i)), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(bessel_Yn, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} bessel_Yn (@var{n}, @var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} bessel_Yn (@dots{})\n\ +\n\ +These routines compute the irregular cylindrical Bessel function of \n\ +order n, Y_n(x), for x>0.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_Yn (static_cast<int>(n.xelem(i)), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_Yn_e (static_cast<int>(n.xelem(i)), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + int nint = static_cast<int>(n.xelem(0)); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_Yn (nint, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_Yn_e (nint, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_Yn (static_cast<int>(n.xelem(i)), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_Yn_e (static_cast<int>(n.xelem(i)), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(bessel_In, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} bessel_In (@var{n}, @var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} bessel_In (@dots{})\n\ +\n\ +These routines compute the regular modified cylindrical Bessel\n\ +function of order n, I_n(x).\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_In (static_cast<int>(n.xelem(i)), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_In_e (static_cast<int>(n.xelem(i)), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + int nint = static_cast<int>(n.xelem(0)); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_In (nint, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_In_e (nint, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_In (static_cast<int>(n.xelem(i)), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_In_e (static_cast<int>(n.xelem(i)), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(bessel_In_scaled, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} bessel_In_scaled (@var{n}, @var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} bessel_In_scaled (@dots{})\n\ +\n\ +These routines compute the scaled regular modified cylindrical Bessel\n\ +function of order n, \\exp(-|x|) I_n(x)\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_In_scaled (static_cast<int>(n.xelem(i)), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_In_scaled_e (static_cast<int>(n.xelem(i)), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + int nint = static_cast<int>(n.xelem(0)); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_In_scaled (nint, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_In_scaled_e (nint, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_In_scaled (static_cast<int>(n.xelem(i)), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_In_scaled_e (static_cast<int>(n.xelem(i)), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(bessel_Kn, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} bessel_Kn (@var{n}, @var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} bessel_Kn (@dots{})\n\ +\n\ +These routines compute the irregular modified cylindrical Bessel\n\ +function of order n, K_n(x), for x > 0.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_Kn (static_cast<int>(n.xelem(i)), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_Kn_e (static_cast<int>(n.xelem(i)), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + int nint = static_cast<int>(n.xelem(0)); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_Kn (nint, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_Kn_e (nint, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_Kn (static_cast<int>(n.xelem(i)), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_Kn_e (static_cast<int>(n.xelem(i)), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(bessel_Kn_scaled, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} bessel_Kn_scaled (@var{n}, @var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} bessel_Kn_scaled (@dots{})\n\ +\n\ +\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_Kn_scaled (static_cast<int>(n.xelem(i)), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_Kn_scaled_e (static_cast<int>(n.xelem(i)), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + int nint = static_cast<int>(n.xelem(0)); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_Kn_scaled (nint, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_Kn_scaled_e (nint, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_Kn_scaled (static_cast<int>(n.xelem(i)), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_Kn_scaled_e (static_cast<int>(n.xelem(i)), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(bessel_jl, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} bessel_jl (@var{n}, @var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} bessel_jl (@dots{})\n\ +\n\ +These routines compute the regular spherical Bessel function of \n\ +order l, j_l(x), for l >= 0 and x >= 0.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_jl (static_cast<int>(n.xelem(i)), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_jl_e (static_cast<int>(n.xelem(i)), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + int nint = static_cast<int>(n.xelem(0)); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_jl (nint, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_jl_e (nint, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_jl (static_cast<int>(n.xelem(i)), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_jl_e (static_cast<int>(n.xelem(i)), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(bessel_yl, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} bessel_yl (@var{n}, @var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} bessel_yl (@dots{})\n\ +\n\ +These routines compute the irregular spherical Bessel function of\n\ +order l, y_l(x), for l >= 0.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_yl (static_cast<int>(n.xelem(i)), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_yl_e (static_cast<int>(n.xelem(i)), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + int nint = static_cast<int>(n.xelem(0)); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_yl (nint, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_yl_e (nint, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_yl (static_cast<int>(n.xelem(i)), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_yl_e (static_cast<int>(n.xelem(i)), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(bessel_il_scaled, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} bessel_il_scaled (@var{n}, @var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} bessel_il_scaled (@dots{})\n\ +\n\ +These routines compute the scaled regular modified spherical Bessel\n\ +function of order l, \\exp(-|x|) i_l(x)\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_il_scaled (static_cast<int>(n.xelem(i)), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_il_scaled_e (static_cast<int>(n.xelem(i)), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + int nint = static_cast<int>(n.xelem(0)); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_il_scaled (nint, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_il_scaled_e (nint, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_il_scaled (static_cast<int>(n.xelem(i)), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_il_scaled_e (static_cast<int>(n.xelem(i)), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(bessel_kl_scaled, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} bessel_kl_scaled (@var{n}, @var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} bessel_kl_scaled (@dots{})\n\ +\n\ +These routines compute the scaled irregular modified spherical Bessel\n\ +function of order l, \\exp(x) k_l(x), for x>0.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_kl_scaled (static_cast<int>(n.xelem(i)), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_kl_scaled_e (static_cast<int>(n.xelem(i)), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + int nint = static_cast<int>(n.xelem(0)); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_kl_scaled (nint, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_kl_scaled_e (nint, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_kl_scaled (static_cast<int>(n.xelem(i)), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_kl_scaled_e (static_cast<int>(n.xelem(i)), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(exprel_n, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} exprel_n (@var{n}, @var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} exprel_n (@dots{})\n\ +\n\ +These routines compute the N-relative exponential, which is the n-th\n\ +generalization of the functions gsl_sf_exprel and gsl_sf_exprel2. The\n\ +N-relative exponential is given by,\n\ +\n\ +exprel_N(x) = N!/x^N (\\exp(x) - \\sum_@{k=0@}^@{N-1@} x^k/k!)\n\ + = 1 + x/(N+1) + x^2/((N+1)(N+2)) + ...\n\ + = 1F1 (1,1+N,x)\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_exprel_n (static_cast<int>(n.xelem(i)), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_exprel_n_e (static_cast<int>(n.xelem(i)), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + int nint = static_cast<int>(n.xelem(0)); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_exprel_n (nint, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_exprel_n_e (nint, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_exprel_n (static_cast<int>(n.xelem(i)), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_exprel_n_e (static_cast<int>(n.xelem(i)), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(fermi_dirac_int, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} fermi_dirac_int (@var{n}, @var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} fermi_dirac_int (@dots{})\n\ +\n\ +These routines compute the complete Fermi-Dirac integral with an\n\ +integer index of j, F_j(x) = (1/\\Gamma(j+1)) \\int_0^\\infty dt (t^j\n\ +/(\\exp(t-x)+1)).\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_fermi_dirac_int (static_cast<int>(n.xelem(i)), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_fermi_dirac_int_e (static_cast<int>(n.xelem(i)), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + int nint = static_cast<int>(n.xelem(0)); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_fermi_dirac_int (nint, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_fermi_dirac_int_e (nint, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_fermi_dirac_int (static_cast<int>(n.xelem(i)), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_fermi_dirac_int_e (static_cast<int>(n.xelem(i)), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(taylorcoeff, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} taylorcoeff (@var{n}, @var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} taylorcoeff (@dots{})\n\ +\n\ +These routines compute the Taylor coefficient x^n / n! \n\ +for x >= 0, n >= 0.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_taylorcoeff (static_cast<int>(n.xelem(i)), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_taylorcoeff_e (static_cast<int>(n.xelem(i)), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + int nint = static_cast<int>(n.xelem(0)); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_taylorcoeff (nint, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_taylorcoeff_e (nint, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_taylorcoeff (static_cast<int>(n.xelem(i)), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_taylorcoeff_e (static_cast<int>(n.xelem(i)), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(legendre_Pl, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} legendre_Pl (@var{n}, @var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} legendre_Pl (@dots{})\n\ +\n\ +These functions evaluate the Legendre polynomial P_l(x) for a specific\n\ +value of l, x subject to l >= 0, |x| <= 1\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_legendre_Pl (static_cast<int>(n.xelem(i)), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_legendre_Pl_e (static_cast<int>(n.xelem(i)), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + int nint = static_cast<int>(n.xelem(0)); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_legendre_Pl (nint, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_legendre_Pl_e (nint, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_legendre_Pl (static_cast<int>(n.xelem(i)), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_legendre_Pl_e (static_cast<int>(n.xelem(i)), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(legendre_Ql, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} legendre_Ql (@var{n}, @var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} legendre_Ql (@dots{})\n\ +\n\ +These routines compute the Legendre function Q_l(x) for x > -1, x != 1\n\ +and l >= 0.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_legendre_Ql (static_cast<int>(n.xelem(i)), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_legendre_Ql_e (static_cast<int>(n.xelem(i)), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + int nint = static_cast<int>(n.xelem(0)); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_legendre_Ql (nint, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_legendre_Ql_e (nint, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_legendre_Ql (static_cast<int>(n.xelem(i)), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_legendre_Ql_e (static_cast<int>(n.xelem(i)), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(psi_n, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} psi_n (@var{n}, @var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} psi_n (@dots{})\n\ +\n\ +These routines compute the polygamma function \\psi^@{(m)@}(x) \n\ +for m >= 0, x > 0.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_psi_n (static_cast<int>(n.xelem(i)), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_psi_n_e (static_cast<int>(n.xelem(i)), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + int nint = static_cast<int>(n.xelem(0)); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_psi_n (nint, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_psi_n_e (nint, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_psi_n (static_cast<int>(n.xelem(i)), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_psi_n_e (static_cast<int>(n.xelem(i)), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(bessel_Jnu, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{z} =} bessel_Jnu (@var{x}, @var{y})\n\ +@deftypefnx {Loadable Function} {[@var{z}, @var{err}] =} bessel_Jnu (@dots{})\n\ +\n\ +These routines compute the regular cylindrical Bessel function of\n\ +fractional order nu, J_\\nu(x).\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{z}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_Jnu (n.xelem(i), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_Jnu_e (n.xelem(i), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + double ndouble = n.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_Jnu (ndouble, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_Jnu_e (ndouble, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_Jnu (n.xelem(i), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_Jnu_e (n.xelem(i), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(bessel_Ynu, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{z} =} bessel_Ynu (@var{x}, @var{y})\n\ +@deftypefnx {Loadable Function} {[@var{z}, @var{err}] =} bessel_Ynu (@dots{})\n\ +\n\ +These routines compute the irregular cylindrical Bessel function of\n\ +fractional order nu, Y_\\nu(x).\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{z}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_Ynu (n.xelem(i), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_Ynu_e (n.xelem(i), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + double ndouble = n.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_Ynu (ndouble, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_Ynu_e (ndouble, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_Ynu (n.xelem(i), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_Ynu_e (n.xelem(i), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(bessel_Inu, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{z} =} bessel_Inu (@var{x}, @var{y})\n\ +@deftypefnx {Loadable Function} {[@var{z}, @var{err}] =} bessel_Inu (@dots{})\n\ +\n\ +These routines compute the regular modified Bessel function of\n\ +fractional order nu, I_\\nu(x) for x>0, \\nu>0.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{z}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_Inu (n.xelem(i), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_Inu_e (n.xelem(i), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + double ndouble = n.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_Inu (ndouble, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_Inu_e (ndouble, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_Inu (n.xelem(i), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_Inu_e (n.xelem(i), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(bessel_Inu_scaled, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{z} =} bessel_Inu_scaled (@var{x}, @var{y})\n\ +@deftypefnx {Loadable Function} {[@var{z}, @var{err}] =} bessel_Inu_scaled (@dots{})\n\ +\n\ +These routines compute the scaled regular modified Bessel function of\n\ +fractional order nu, \\exp(-|x|)I_\\nu(x) for x>0, \\nu>0.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{z}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_Inu_scaled (n.xelem(i), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_Inu_scaled_e (n.xelem(i), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + double ndouble = n.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_Inu_scaled (ndouble, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_Inu_scaled_e (ndouble, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_Inu_scaled (n.xelem(i), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_Inu_scaled_e (n.xelem(i), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(bessel_Knu, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{z} =} bessel_Knu (@var{x}, @var{y})\n\ +@deftypefnx {Loadable Function} {[@var{z}, @var{err}] =} bessel_Knu (@dots{})\n\ +\n\ +These routines compute the irregular modified Bessel function of\n\ +fractional order nu, K_\\nu(x) for x>0, \\nu>0.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{z}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_Knu (n.xelem(i), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_Knu_e (n.xelem(i), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + double ndouble = n.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_Knu (ndouble, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_Knu_e (ndouble, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_Knu (n.xelem(i), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_Knu_e (n.xelem(i), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(bessel_lnKnu, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{z} =} bessel_lnKnu (@var{x}, @var{y})\n\ +@deftypefnx {Loadable Function} {[@var{z}, @var{err}] =} bessel_lnKnu (@dots{})\n\ +\n\ +These routines compute the logarithm of the irregular modified Bessel\n\ +function of fractional order nu, \\ln(K_\\nu(x)) for x>0, \\nu>0.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{z}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_lnKnu (n.xelem(i), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_lnKnu_e (n.xelem(i), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + double ndouble = n.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_lnKnu (ndouble, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_lnKnu_e (ndouble, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_lnKnu (n.xelem(i), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_lnKnu_e (n.xelem(i), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(bessel_Knu_scaled, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{z} =} bessel_Knu_scaled (@var{x}, @var{y})\n\ +@deftypefnx {Loadable Function} {[@var{z}, @var{err}] =} bessel_Knu_scaled (@dots{})\n\ +\n\ +These routines compute the scaled irregular modified Bessel function\n\ +of fractional order nu, \\exp(+|x|) K_\\nu(x) for x>0, \\nu>0.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{z}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_Knu_scaled (n.xelem(i), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_Knu_scaled_e (n.xelem(i), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + double ndouble = n.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_Knu_scaled (ndouble, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_Knu_scaled_e (ndouble, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_bessel_Knu_scaled (n.xelem(i), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_bessel_Knu_scaled_e (n.xelem(i), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(exp_mult, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{z} =} exp_mult (@var{x}, @var{y})\n\ +@deftypefnx {Loadable Function} {[@var{z}, @var{err}] =} exp_mult (@dots{})\n\ +\n\ +These routines exponentiate x and multiply by the factor y to return\n\ +the product y \\exp(x).\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{z}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_exp_mult (n.xelem(i), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_exp_mult_e (n.xelem(i), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + double ndouble = n.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_exp_mult (ndouble, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_exp_mult_e (ndouble, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_exp_mult (n.xelem(i), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_exp_mult_e (n.xelem(i), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(fermi_dirac_inc_0, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{z} =} fermi_dirac_inc_0 (@var{x}, @var{y})\n\ +@deftypefnx {Loadable Function} {[@var{z}, @var{err}] =} fermi_dirac_inc_0 (@dots{})\n\ +\n\ +These routines compute the incomplete Fermi-Dirac integral with an\n\ +index of zero, F_0(x,b) = \\ln(1 + e^@{b-x@}) - (b-x).\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{z}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_fermi_dirac_inc_0 (n.xelem(i), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_fermi_dirac_inc_0_e (n.xelem(i), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + double ndouble = n.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_fermi_dirac_inc_0 (ndouble, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_fermi_dirac_inc_0_e (ndouble, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_fermi_dirac_inc_0 (n.xelem(i), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_fermi_dirac_inc_0_e (n.xelem(i), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(poch, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{z} =} poch (@var{x}, @var{y})\n\ +@deftypefnx {Loadable Function} {[@var{z}, @var{err}] =} poch (@dots{})\n\ +\n\ +These routines compute the Pochhammer symbol \n\ +\n\ +(a)_x := \\Gamma(a + x)/\\Gamma(a), \n\ +\n\ +subject to a and a+x not being negative integers. The Pochhammer\n\ +symbol is also known as the Apell symbol.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{z}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_poch (n.xelem(i), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_poch_e (n.xelem(i), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + double ndouble = n.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_poch (ndouble, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_poch_e (ndouble, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_poch (n.xelem(i), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_poch_e (n.xelem(i), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(lnpoch, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{z} =} lnpoch (@var{x}, @var{y})\n\ +@deftypefnx {Loadable Function} {[@var{z}, @var{err}] =} lnpoch (@dots{})\n\ +\n\ +These routines compute the logarithm of the Pochhammer symbol,\n\ +\\log((a)_x) = \\log(\\Gamma(a + x)/\\Gamma(a)) for a > 0, a+x > 0.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{z}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_lnpoch (n.xelem(i), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_lnpoch_e (n.xelem(i), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + double ndouble = n.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_lnpoch (ndouble, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_lnpoch_e (ndouble, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_lnpoch (n.xelem(i), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_lnpoch_e (n.xelem(i), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(pochrel, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{z} =} pochrel (@var{x}, @var{y})\n\ +@deftypefnx {Loadable Function} {[@var{z}, @var{err}] =} pochrel (@dots{})\n\ +\n\ +These routines compute the relative Pochhammer symbol ((a,x) - 1)/x\n\ +where (a,x) = (a)_x := \\Gamma(a + x)/\\Gamma(a).\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{z}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_pochrel (n.xelem(i), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_pochrel_e (n.xelem(i), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + double ndouble = n.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_pochrel (ndouble, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_pochrel_e (ndouble, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_pochrel (n.xelem(i), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_pochrel_e (n.xelem(i), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(gamma_inc_Q, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{z} =} gamma_inc_Q (@var{x}, @var{y})\n\ +@deftypefnx {Loadable Function} {[@var{z}, @var{err}] =} gamma_inc_Q (@dots{})\n\ +\n\ +These routines compute the normalized incomplete Gamma Function \n\ +Q(a,x) = 1/\\Gamma(a) \\int_x\\infty dt t^@{a-1@} \\exp(-t) for a > 0, x >= 0.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{z}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_gamma_inc_Q (n.xelem(i), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_gamma_inc_Q_e (n.xelem(i), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + double ndouble = n.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_gamma_inc_Q (ndouble, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_gamma_inc_Q_e (ndouble, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_gamma_inc_Q (n.xelem(i), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_gamma_inc_Q_e (n.xelem(i), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(gamma_inc_P, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{z} =} gamma_inc_P (@var{x}, @var{y})\n\ +@deftypefnx {Loadable Function} {[@var{z}, @var{err}] =} gamma_inc_P (@dots{})\n\ +\n\ +These routines compute the complementary normalized incomplete Gamma\n\ +Function P(a,x) = 1/\\Gamma(a) \\int_0^x dt t^@{a-1@} \\exp(-t) \n\ +for a > 0, x >= 0.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{z}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_gamma_inc_P (n.xelem(i), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_gamma_inc_P_e (n.xelem(i), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + double ndouble = n.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_gamma_inc_P (ndouble, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_gamma_inc_P_e (ndouble, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_gamma_inc_P (n.xelem(i), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_gamma_inc_P_e (n.xelem(i), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(gamma_inc, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{z} =} gamma_inc (@var{x}, @var{y})\n\ +@deftypefnx {Loadable Function} {[@var{z}, @var{err}] =} gamma_inc (@dots{})\n\ +\n\ +These functions compute the incomplete Gamma Function the\n\ +normalization factor included in the previously defined functions:\n\ +\\Gamma(a,x) = \\int_x\\infty dt t^@{a-1@} \\exp(-t) for a real and x >= 0.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{z}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_gamma_inc (n.xelem(i), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_gamma_inc_e (n.xelem(i), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + double ndouble = n.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_gamma_inc (ndouble, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_gamma_inc_e (ndouble, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_gamma_inc (n.xelem(i), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_gamma_inc_e (n.xelem(i), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(beta_gsl, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{z} =} beta_gsl (@var{x}, @var{y})\n\ +@deftypefnx {Loadable Function} {[@var{z}, @var{err}] =} beta_gsl (@dots{})\n\ +\n\ +These routines compute the Beta Function, \n\ +B(a,b) = \\Gamma(a)\\Gamma(b)/\\Gamma(a+b) for a > 0, b > 0.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{z}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_beta (n.xelem(i), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_beta_e (n.xelem(i), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + double ndouble = n.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_beta (ndouble, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_beta_e (ndouble, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_beta (n.xelem(i), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_beta_e (n.xelem(i), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(lnbeta, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{z} =} lnbeta (@var{x}, @var{y})\n\ +@deftypefnx {Loadable Function} {[@var{z}, @var{err}] =} lnbeta (@dots{})\n\ +\n\ +These routines compute the logarithm of the Beta Function,\n\ +\\log(B(a,b)) for a > 0, b > 0.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{z}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_lnbeta (n.xelem(i), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_lnbeta_e (n.xelem(i), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + double ndouble = n.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_lnbeta (ndouble, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_lnbeta_e (ndouble, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_lnbeta (n.xelem(i), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_lnbeta_e (n.xelem(i), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(hyperg_0F1, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{z} =} hyperg_0F1 (@var{x}, @var{y})\n\ +@deftypefnx {Loadable Function} {[@var{z}, @var{err}] =} hyperg_0F1 (@dots{})\n\ +\n\ +These routines compute the hypergeometric function 0F1(c,x).\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{z}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_hyperg_0F1 (n.xelem(i), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_hyperg_0F1_e (n.xelem(i), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + double ndouble = n.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_hyperg_0F1 (ndouble, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_hyperg_0F1_e (ndouble, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_hyperg_0F1 (n.xelem(i), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_hyperg_0F1_e (n.xelem(i), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(conicalP_half, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{z} =} conicalP_half (@var{x}, @var{y})\n\ +@deftypefnx {Loadable Function} {[@var{z}, @var{err}] =} conicalP_half (@dots{})\n\ +\n\ +These routines compute the irregular Spherical Conical Function\n\ +P^@{1/2@}_@{-1/2 + i \\lambda@}(x) for x > -1.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{z}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_conicalP_half (n.xelem(i), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_conicalP_half_e (n.xelem(i), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + double ndouble = n.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_conicalP_half (ndouble, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_conicalP_half_e (ndouble, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_conicalP_half (n.xelem(i), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_conicalP_half_e (n.xelem(i), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(conicalP_mhalf, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{z} =} conicalP_mhalf (@var{x}, @var{y})\n\ +@deftypefnx {Loadable Function} {[@var{z}, @var{err}] =} conicalP_mhalf (@dots{})\n\ +\n\ +These routines compute the regular Spherical Conical Function\n\ +P^@{-1/2@}_@{-1/2 + i \\lambda@}(x) for x > -1.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{z}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_conicalP_mhalf (n.xelem(i), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_conicalP_mhalf_e (n.xelem(i), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + double ndouble = n.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_conicalP_mhalf (ndouble, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_conicalP_mhalf_e (ndouble, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_conicalP_mhalf (n.xelem(i), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_conicalP_mhalf_e (n.xelem(i), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(conicalP_0, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{z} =} conicalP_0 (@var{x}, @var{y})\n\ +@deftypefnx {Loadable Function} {[@var{z}, @var{err}] =} conicalP_0 (@dots{})\n\ +\n\ +These routines compute the conical function P^0_@{-1/2 + i \\lambda@}(x)\n\ +for x > -1.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{z}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_conicalP_0 (n.xelem(i), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_conicalP_0_e (n.xelem(i), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + double ndouble = n.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_conicalP_0 (ndouble, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_conicalP_0_e (ndouble, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_conicalP_0 (n.xelem(i), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_conicalP_0_e (n.xelem(i), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(conicalP_1, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{z} =} conicalP_1 (@var{x}, @var{y})\n\ +@deftypefnx {Loadable Function} {[@var{z}, @var{err}] =} conicalP_1 (@dots{})\n\ +\n\ +These routines compute the conical function P^1_@{-1/2 + i \\lambda@}(x)\n\ +for x > -1.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{z}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_conicalP_1 (n.xelem(i), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_conicalP_1_e (n.xelem(i), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + double ndouble = n.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_conicalP_1 (ndouble, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_conicalP_1_e (ndouble, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_conicalP_1 (n.xelem(i), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_conicalP_1_e (n.xelem(i), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(hzeta, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{z} =} hzeta (@var{x}, @var{y})\n\ +@deftypefnx {Loadable Function} {[@var{z}, @var{err}] =} hzeta (@dots{})\n\ +\n\ +These routines compute the Hurwitz zeta function \\zeta(s,q) \n\ +for s > 1, q > 0.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{z}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray x = args(1).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_hzeta (n.xelem(i), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_hzeta_e (n.xelem(i), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + double ndouble = n.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_hzeta (ndouble, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_hzeta_e (ndouble, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_hzeta (n.xelem(i), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_hzeta_e (n.xelem(i), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("First and second argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(airy_Ai, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} airy_Ai (@var{x}, @var{mode})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} airy_Ai (@dots{})\n\ +\n\ +These routines compute the Airy function Ai(x) with an accuracy\n\ +specified by mode.\n\ +\n\ +The second argument @var{mode} must be an integer corresponding to\n\ +\n\ +@table @asis\n\ +@item 0 = GSL_PREC_DOUBLE\n\ + Double-precision, a relative accuracy of approximately @code{2 * 10^-16}.\n\ +@item 1 = GSL_PREC_SINGLE\n\ + Single-precision, a relative accuracy of approximately @code{10^-7}.\n\ +@item 2 = GSL_PREC_APPROX\n\ + Approximate values, a relative accuracy of approximately @code{5 * 10^-4}.\n\ +@end table\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + if(!args(1).is_scalar_type()) { + error("The mode must be scalar."); + print_usage (); + return octave_value(); + } + int mode = static_cast<int>((args(1).array_value())(0)); + if(mode < 0) + mode = 0; + else if(mode > 2) + mode = 2; + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_airy_Ai (x.xelem(i), mode); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_airy_Ai_e (x.xelem(i), mode, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(airy_Bi, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} airy_Bi (@var{x}, @var{mode})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} airy_Bi (@dots{})\n\ +\n\ +These routines compute the Airy function Bi(x) with an accuracy\n\ +specified by mode.\n\ +\n\ +The second argument @var{mode} must be an integer corresponding to\n\ +\n\ +@table @asis\n\ +@item 0 = GSL_PREC_DOUBLE\n\ + Double-precision, a relative accuracy of approximately @code{2 * 10^-16}.\n\ +@item 1 = GSL_PREC_SINGLE\n\ + Single-precision, a relative accuracy of approximately @code{10^-7}.\n\ +@item 2 = GSL_PREC_APPROX\n\ + Approximate values, a relative accuracy of approximately @code{5 * 10^-4}.\n\ +@end table\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + if(!args(1).is_scalar_type()) { + error("The mode must be scalar."); + print_usage (); + return octave_value(); + } + int mode = static_cast<int>((args(1).array_value())(0)); + if(mode < 0) + mode = 0; + else if(mode > 2) + mode = 2; + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_airy_Bi (x.xelem(i), mode); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_airy_Bi_e (x.xelem(i), mode, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(airy_Ai_scaled, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} airy_Ai_scaled (@var{x}, @var{mode})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} airy_Ai_scaled (@dots{})\n\ +\n\ +These routines compute a scaled version of the Airy function \n\ +S_A(x) Ai(x). For x>0 the scaling factor S_A(x) is \\exp(+(2/3) x^(3/2)), and\n\ +is 1 for x<0.\n\ +\n\ +The second argument @var{mode} must be an integer corresponding to\n\ +\n\ +@table @asis\n\ +@item 0 = GSL_PREC_DOUBLE\n\ + Double-precision, a relative accuracy of approximately @code{2 * 10^-16}.\n\ +@item 1 = GSL_PREC_SINGLE\n\ + Single-precision, a relative accuracy of approximately @code{10^-7}.\n\ +@item 2 = GSL_PREC_APPROX\n\ + Approximate values, a relative accuracy of approximately @code{5 * 10^-4}.\n\ +@end table\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + if(!args(1).is_scalar_type()) { + error("The mode must be scalar."); + print_usage (); + return octave_value(); + } + int mode = static_cast<int>((args(1).array_value())(0)); + if(mode < 0) + mode = 0; + else if(mode > 2) + mode = 2; + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_airy_Ai_scaled (x.xelem(i), mode); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_airy_Ai_scaled_e (x.xelem(i), mode, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(airy_Bi_scaled, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} airy_Bi_scaled (@var{x}, @var{mode})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} airy_Bi_scaled (@dots{})\n\ +\n\ +These routines compute a scaled version of the Airy function \n\ +S_B(x) Bi(x). For x>0 the scaling factor S_B(x) is exp(-(2/3) x^(3/2)), and\n\ +is 1 for x<0.\n\ +\n\ +The second argument @var{mode} must be an integer corresponding to\n\ +\n\ +@table @asis\n\ +@item 0 = GSL_PREC_DOUBLE\n\ + Double-precision, a relative accuracy of approximately @code{2 * 10^-16}.\n\ +@item 1 = GSL_PREC_SINGLE\n\ + Single-precision, a relative accuracy of approximately @code{10^-7}.\n\ +@item 2 = GSL_PREC_APPROX\n\ + Approximate values, a relative accuracy of approximately @code{5 * 10^-4}.\n\ +@end table\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + if(!args(1).is_scalar_type()) { + error("The mode must be scalar."); + print_usage (); + return octave_value(); + } + int mode = static_cast<int>((args(1).array_value())(0)); + if(mode < 0) + mode = 0; + else if(mode > 2) + mode = 2; + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_airy_Bi_scaled (x.xelem(i), mode); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_airy_Bi_scaled_e (x.xelem(i), mode, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(airy_Ai_deriv, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} airy_Ai_deriv (@var{x}, @var{mode})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} airy_Ai_deriv (@dots{})\n\ +\n\ +These routines compute the Airy function derivative Ai'(x) with an\n\ +accuracy specified by mode.\n\ +\n\ +The second argument @var{mode} must be an integer corresponding to\n\ +\n\ +@table @asis\n\ +@item 0 = GSL_PREC_DOUBLE\n\ + Double-precision, a relative accuracy of approximately @code{2 * 10^-16}.\n\ +@item 1 = GSL_PREC_SINGLE\n\ + Single-precision, a relative accuracy of approximately @code{10^-7}.\n\ +@item 2 = GSL_PREC_APPROX\n\ + Approximate values, a relative accuracy of approximately @code{5 * 10^-4}.\n\ +@end table\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + if(!args(1).is_scalar_type()) { + error("The mode must be scalar."); + print_usage (); + return octave_value(); + } + int mode = static_cast<int>((args(1).array_value())(0)); + if(mode < 0) + mode = 0; + else if(mode > 2) + mode = 2; + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_airy_Ai_deriv (x.xelem(i), mode); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_airy_Ai_deriv_e (x.xelem(i), mode, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(airy_Bi_deriv, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} airy_Bi_deriv (@var{x}, @var{mode})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} airy_Bi_deriv (@dots{})\n\ +\n\ +These routines compute the Airy function derivative Bi'(x) with an\n\ +accuracy specified by mode.\n\ +\n\ +The second argument @var{mode} must be an integer corresponding to\n\ +\n\ +@table @asis\n\ +@item 0 = GSL_PREC_DOUBLE\n\ + Double-precision, a relative accuracy of approximately @code{2 * 10^-16}.\n\ +@item 1 = GSL_PREC_SINGLE\n\ + Single-precision, a relative accuracy of approximately @code{10^-7}.\n\ +@item 2 = GSL_PREC_APPROX\n\ + Approximate values, a relative accuracy of approximately @code{5 * 10^-4}.\n\ +@end table\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + if(!args(1).is_scalar_type()) { + error("The mode must be scalar."); + print_usage (); + return octave_value(); + } + int mode = static_cast<int>((args(1).array_value())(0)); + if(mode < 0) + mode = 0; + else if(mode > 2) + mode = 2; + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_airy_Bi_deriv (x.xelem(i), mode); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_airy_Bi_deriv_e (x.xelem(i), mode, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(airy_Ai_deriv_scaled, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} airy_Ai_deriv_scaled (@var{x}, @var{mode})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} airy_Ai_deriv_scaled (@dots{})\n\ +\n\ +These routines compute the derivative of the scaled Airy function\n\ +S_A(x) Ai(x).\n\ +\n\ +The second argument @var{mode} must be an integer corresponding to\n\ +\n\ +@table @asis\n\ +@item 0 = GSL_PREC_DOUBLE\n\ + Double-precision, a relative accuracy of approximately @code{2 * 10^-16}.\n\ +@item 1 = GSL_PREC_SINGLE\n\ + Single-precision, a relative accuracy of approximately @code{10^-7}.\n\ +@item 2 = GSL_PREC_APPROX\n\ + Approximate values, a relative accuracy of approximately @code{5 * 10^-4}.\n\ +@end table\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + if(!args(1).is_scalar_type()) { + error("The mode must be scalar."); + print_usage (); + return octave_value(); + } + int mode = static_cast<int>((args(1).array_value())(0)); + if(mode < 0) + mode = 0; + else if(mode > 2) + mode = 2; + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_airy_Ai_deriv_scaled (x.xelem(i), mode); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_airy_Ai_deriv_scaled_e (x.xelem(i), mode, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(airy_Bi_deriv_scaled, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} airy_Bi_deriv_scaled (@var{x}, @var{mode})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} airy_Bi_deriv_scaled (@dots{})\n\ +\n\ +These routines compute the derivative of the scaled Airy function\n\ +S_B(x) Bi(x).\n\ +\n\ +The second argument @var{mode} must be an integer corresponding to\n\ +\n\ +@table @asis\n\ +@item 0 = GSL_PREC_DOUBLE\n\ + Double-precision, a relative accuracy of approximately @code{2 * 10^-16}.\n\ +@item 1 = GSL_PREC_SINGLE\n\ + Single-precision, a relative accuracy of approximately @code{10^-7}.\n\ +@item 2 = GSL_PREC_APPROX\n\ + Approximate values, a relative accuracy of approximately @code{5 * 10^-4}.\n\ +@end table\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + if(!args(1).is_scalar_type()) { + error("The mode must be scalar."); + print_usage (); + return octave_value(); + } + int mode = static_cast<int>((args(1).array_value())(0)); + if(mode < 0) + mode = 0; + else if(mode > 2) + mode = 2; + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_airy_Bi_deriv_scaled (x.xelem(i), mode); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_airy_Bi_deriv_scaled_e (x.xelem(i), mode, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(ellint_Kcomp, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} ellint_Kcomp (@var{x}, @var{mode})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} ellint_Kcomp (@dots{})\n\ +\n\ +These routines compute the complete elliptic integral K(k) to the\n\ +accuracy specified by the mode variable mode.\n\ +\n\ +The second argument @var{mode} must be an integer corresponding to\n\ +\n\ +@table @asis\n\ +@item 0 = GSL_PREC_DOUBLE\n\ + Double-precision, a relative accuracy of approximately @code{2 * 10^-16}.\n\ +@item 1 = GSL_PREC_SINGLE\n\ + Single-precision, a relative accuracy of approximately @code{10^-7}.\n\ +@item 2 = GSL_PREC_APPROX\n\ + Approximate values, a relative accuracy of approximately @code{5 * 10^-4}.\n\ +@end table\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + if(!args(1).is_scalar_type()) { + error("The mode must be scalar."); + print_usage (); + return octave_value(); + } + int mode = static_cast<int>((args(1).array_value())(0)); + if(mode < 0) + mode = 0; + else if(mode > 2) + mode = 2; + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_ellint_Kcomp (x.xelem(i), mode); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_ellint_Kcomp_e (x.xelem(i), mode, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(ellint_Ecomp, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} ellint_Ecomp (@var{x}, @var{mode})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} ellint_Ecomp (@dots{})\n\ +\n\ +These routines compute the complete elliptic integral E(k) to the\n\ +accuracy specified by the mode variable mode.\n\ +\n\ +The second argument @var{mode} must be an integer corresponding to\n\ +\n\ +@table @asis\n\ +@item 0 = GSL_PREC_DOUBLE\n\ + Double-precision, a relative accuracy of approximately @code{2 * 10^-16}.\n\ +@item 1 = GSL_PREC_SINGLE\n\ + Single-precision, a relative accuracy of approximately @code{10^-7}.\n\ +@item 2 = GSL_PREC_APPROX\n\ + Approximate values, a relative accuracy of approximately @code{5 * 10^-4}.\n\ +@end table\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 2) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + if(!args(1).is_scalar_type()) { + error("The mode must be scalar."); + print_usage (); + return octave_value(); + } + int mode = static_cast<int>((args(1).array_value())(0)); + if(mode < 0) + mode = 0; + else if(mode > 2) + mode = 2; + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_ellint_Ecomp (x.xelem(i), mode); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_ellint_Ecomp_e (x.xelem(i), mode, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(airy_zero_Ai, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} airy_zero_Ai (@var{n})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} airy_zero_Ai (@dots{})\n\ +\n\ +These routines compute the location of the s-th zero of the Airy\n\ +function Ai(x).\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_airy_zero_Ai (static_cast<int>(x.xelem(i))); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_airy_zero_Ai_e (static_cast<int>(x.xelem(i)), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value_list(); + +} + + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(airy_zero_Bi, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} airy_zero_Bi (@var{n})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} airy_zero_Bi (@dots{})\n\ +\n\ +These routines compute the location of the s-th zero of the Airy\n\ +function Bi(x).\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_airy_zero_Bi (static_cast<int>(x.xelem(i))); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_airy_zero_Bi_e (static_cast<int>(x.xelem(i)), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value_list(); + +} + + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(airy_zero_Ai_deriv, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} airy_zero_Ai_deriv (@var{n})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} airy_zero_Ai_deriv (@dots{})\n\ +\n\ +These routines compute the location of the s-th zero of the Airy\n\ +function derivative Ai'(x).\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_airy_zero_Ai_deriv (static_cast<int>(x.xelem(i))); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_airy_zero_Ai_deriv_e (static_cast<int>(x.xelem(i)), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value_list(); + +} + + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(airy_zero_Bi_deriv, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} airy_zero_Bi_deriv (@var{n})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} airy_zero_Bi_deriv (@dots{})\n\ +\n\ +These routines compute the location of the s-th zero of the Airy\n\ +function derivative Bi'(x).\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_airy_zero_Bi_deriv (static_cast<int>(x.xelem(i))); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_airy_zero_Bi_deriv_e (static_cast<int>(x.xelem(i)), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value_list(); + +} + + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(bessel_zero_J0, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} bessel_zero_J0 (@var{n})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} bessel_zero_J0 (@dots{})\n\ +\n\ +These routines compute the location of the s-th positive zero of the\n\ +Bessel function J_0(x).\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_bessel_zero_J0 (static_cast<int>(x.xelem(i))); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_bessel_zero_J0_e (static_cast<int>(x.xelem(i)), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value_list(); + +} + + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(bessel_zero_J1, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} bessel_zero_J1 (@var{n})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} bessel_zero_J1 (@dots{})\n\ +\n\ +These routines compute the location of the s-th positive zero of the\n\ +Bessel function J_1(x).\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_bessel_zero_J1 (static_cast<int>(x.xelem(i))); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_bessel_zero_J1_e (static_cast<int>(x.xelem(i)), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value_list(); + +} + + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(psi_1_int, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} psi_1_int (@var{n})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} psi_1_int (@dots{})\n\ +\n\ +These routines compute the Trigamma function \\psi'(n) for positive\n\ +integer n.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_psi_1_int (static_cast<int>(x.xelem(i))); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_psi_1_int_e (static_cast<int>(x.xelem(i)), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value_list(); + +} + + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(zeta_int, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} zeta_int (@var{n})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} zeta_int (@dots{})\n\ +\n\ +These routines compute the Riemann zeta function \\zeta(n) for \n\ +integer n, n \\ne 1.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_zeta_int (static_cast<int>(x.xelem(i))); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_zeta_int_e (static_cast<int>(x.xelem(i)), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value_list(); + +} + + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(eta_int, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} eta_int (@var{n})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} eta_int (@dots{})\n\ +\n\ +These routines compute the eta function \\eta(n) for integer n.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 1) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type()) { + error("The argument must be real."); + print_usage (); + return octave_value(); + } + + NDArray x = args(0).array_value(); + NDArray y(x.dims()); + int lx = x.length(); +// printf("length: %d\n", lx); +// printf("nargout: %d\n", nargout); + if(nargout < 2) { + for(i = 0; i < lx; i++) { + y.xelem(i) = gsl_sf_eta_int (static_cast<int>(x.xelem(i))); + } + return octave_value(y); + } else { + NDArray err(x.dims()); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < lx; i++) { + gsl_sf_eta_int_e (static_cast<int>(x.xelem(i)), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + + return octave_value_list(); + +} + + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(legendre_Plm, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} legendre_Plm (@var{n}, @var{m}, @var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} legendre_Plm (@dots{})\n\ +\n\ +These routines compute the associated Legendre polynomial P_l^m(x) \n\ +for m >= 0, l >= m, |x| <= 1.\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 3) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type() + || !args(2).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray m = args(1).array_value(); + if(n.length() != m.length()) { + error("Two first arguments must have the same size."); + print_usage (); + return octave_value(); + } + NDArray x = args(2).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_legendre_Plm + (static_cast<int>(n.xelem(i)), + static_cast<int>(m.xelem(i)), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_legendre_Plm_e (static_cast<int>(n.xelem(i)), + static_cast<int>(m.xelem(i)), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + int nint = static_cast<int>(n.xelem(0)); + int mint = static_cast<int>(m.xelem(0)); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_legendre_Plm (nint, mint, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_legendre_Plm_e (nint, mint, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_legendre_Plm (static_cast<int>(n.xelem(i)), + static_cast<int>(m.xelem(i)), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_legendre_Plm_e (static_cast<int>(n.xelem(i)), + static_cast<int>(m.xelem(i)), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("The two first and the third argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +DEFUN_DLD(legendre_sphPlm, args, nargout, "\ + -*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{y} =} legendre_sphPlm (@var{n}, @var{m}, @var{x})\n\ +@deftypefnx {Loadable Function} {[@var{y}, @var{err}] =} legendre_sphPlm (@dots{})\n\ +\n\ +These routines compute the normalized associated Legendre polynomial\n\ +$\\sqrt@{(2l+1)/(4\\pi)@} \\sqrt@{(l-m)!/(l+m)!@} P_l^m(x)$ suitable for use\n\ +in spherical harmonics. The parameters must satisfy m >= 0, l >= m,\n\ +|x| <= 1. Theses routines avoid the overflows that occur for the\n\ +standard normalization of P_l^m(x).\n\ +\n\ +@var{err} contains an estimate of the absolute error in the value @var{y}.\n\ +\n\ +This function is from the GNU Scientific Library,\n\ +see @url{http://www.gnu.org/software/gsl/} for documentation.\n\ +@end deftypefn\n\ +") +{ + int i; + dim_vector dv; + + gsl_set_error_handler (octave_gsl_errorhandler); + + if(args.length() != 3) { + print_usage (); + return octave_value(); + } + if(!args(0).is_real_type() || !args(1).is_real_type() + || !args(2).is_real_type()) { + error("The arguments must be real."); + print_usage (); + return octave_value(); + } + // Nice combinatorial explosion here + NDArray n = args(0).array_value(); + NDArray m = args(1).array_value(); + if(n.length() != m.length()) { + error("Two first arguments must have the same size."); + print_usage (); + return octave_value(); + } + NDArray x = args(2).array_value(); + if(n.length() == x.length()) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_legendre_sphPlm + (static_cast<int>(n.xelem(i)), + static_cast<int>(m.xelem(i)), + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_legendre_sphPlm_e (static_cast<int>(n.xelem(i)), + static_cast<int>(m.xelem(i)), + x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(n.length() == 1) { + dv = x.dims(); + NDArray y(dv); + int len = x.length(); + int nint = static_cast<int>(n.xelem(0)); + int mint = static_cast<int>(m.xelem(0)); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_legendre_sphPlm (nint, mint, + x.xelem(i)); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_legendre_sphPlm_e (nint, mint, x.xelem(i), &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else if(x.length() == 1) { + dv = n.dims(); + NDArray y(dv); + int len = n.length(); + double xdouble = x.xelem(0); + if(nargout < 2) { + for(i = 0; i < len; i++) { + y.xelem(i) = gsl_sf_legendre_sphPlm (static_cast<int>(n.xelem(i)), + static_cast<int>(m.xelem(i)), + xdouble); + } + return octave_value(y); + } else { + NDArray err(dv); + gsl_sf_result result; + octave_value_list retval; + for(i = 0; i < len; i++) { + gsl_sf_legendre_sphPlm_e (static_cast<int>(n.xelem(i)), + static_cast<int>(m.xelem(i)), + xdouble, &result); + y.xelem(i) = result.val; + err.xelem(i) = result.err; + } + retval(1) = octave_value(err); + retval(0) = octave_value(y); + return retval; + } + } else { + error("The two first and the third argument must either have the same size, or one of them must be scalar."); + print_usage (); + } + + return octave_value(); + +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/