## Copyright (C) 2015 Nicholas R. Jankowski <jankowskin@asme.org> **
##
## This file is part of Octave.
##
## Octave is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING.  If not,
## see <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn {} {} humps (@var{X})
## @deftypefnx {} {[@var{X}, @var{Y}] =} humps (@var{X})
## Return the output of the rational function:
##
## @example
## @group
##         1200*@var{X}^4 - 2880*@var{X}^3 + 2036*@var{X}^2 - 348*@var{X} - 88
##  @var{Y} = - ---------------------------------------------
##          200*@var{X}^4 - 480*@var{X}^3 + 406*@var{X}^2 - 138*@var{X} + 17
## @end group
## @end example
##
## @var{X} may be a scalar, vector or array.  humps is vectorized such
## that 
##   @var{Yi} = humps(@var{Xi})
##
## If @var{X} is omitted, the range [0:0.05:1] is used instead.
##
## When called with two output arguments, [@var{X},@var{Y}], @var{X} will 
## return the input value in @var{X} along with the output from humps in @var{Y}.
##
## humps has two local maxima located near @var{x} = 0.300 and 0.893 
## and a local minimum near @var{x} = 0.637, and is a useful function for 
## testing algorithms to find local minima and maxima.
##
## @end deftypefn

function varargout = humps (x)
  
  if (nargin > 1)
    error("humps: Too many input arguments")
    print_usage ();
  
  elseif (nargin == 0)
    x = [0:0.05:1];
  
  endif

 y = - 4*( 300*x.^4 - 720*x.^3 + 509*x.^2 - 87*x - 22) ./ ...
         ((10*x.^2 - 6*x + 1).*(20*x.^2 - 36*x + 17));

    
  switch nargout
    
    case {0,1}
      varargout{1} = y;
    
    case 2
      varargout{1} = x;
      varargout{2} = y;
      
    otherwise
      error("humps: Too many output arguments requested")
      print_usage ();
    
  endswitch
   
endfunction

%%% tests
%%% value checks
%!assert (humps(0), 88/17, 10*eps)
%!assert (humps(1), 16, 10*eps)
%!assert (humps(-1), -6376/1241, 10*eps)
%!assert (humps(),[88/17, 16106/1769, 263/17, 82802/3133, 2432/53, ...
%!   2818/37, 193/2, 10538/137, 1376/29, 36434/1261, 19, 5258/377,  ...
%!   152/13, 24562/2173, 421/34, 250/17, 232/13, 1762/85, 803/37,   ...
%!   58354/2941, 16], 1000*eps)

%%% vector checks
%!assert (humps([0, 1]), [88/17, 16], 10*eps)
%!assert (humps([0, 1]'), [88/17, 16]', 10*eps)
%!assert (humps([0, 1; 1, 0]'), [88/17, 16; 16, 88/17]', 10*eps)

%%% array checks
%!assert (humps(repmat(eye(2), 1, 1, 2)), repmat([16, 88/17; 88/17, 16], ...
%! 1, 1, 2), 10*eps)