########################################################################
##
## Copyright (C) 2002-2022 The Octave Project Developers
##
## See the file COPYRIGHT.md in the top-level directory of this
## distribution or <https://octave.org/copyright/>.
##
## 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
## <https://www.gnu.org/licenses/>.
##
########################################################################

CLASS = "LSODE"

INCLUDE = "ODE.h"

OPTION
  NAME = "absolute tolerance"
  DOC_ITEM
Absolute tolerance.  May be either vector or scalar.  If a vector, it
must match the dimension of the state vector.

  END_DOC_ITEM
  TYPE = "Array<double>"
  SET_ARG_TYPE = "const $TYPE&"
  INIT_BODY
    $OPTVAR.resize (dim_vector (1, 1));
    $OPTVAR(0) = ::sqrt (std::numeric_limits<double>::epsilon ());
  END_INIT_BODY
  SET_CODE
    void set_$OPT (double val)
      {
        $OPTVAR.resize (dim_vector (1, 1));
        $OPTVAR(0) = (val > 0.0) ? val : ::sqrt (std::numeric_limits<double>::epsilon ());
        m_reset = true;
      }

    void set_$OPT (const $TYPE& val)
      { $OPTVAR = val; m_reset = true; }
  END_SET_CODE
END_OPTION

OPTION
  NAME = "relative tolerance"
  DOC_ITEM
Relative tolerance parameter.  Unlike the absolute tolerance, this
parameter may only be a scalar.

The local error test applied at each integration step is

@example
@group
  abs (local error in x(i)) <= ...
      rtol * abs (y(i)) + atol(i)
@end group
@end example

  END_DOC_ITEM
  TYPE = "double"
  INIT_VALUE = "::sqrt (std::numeric_limits<double>::epsilon ())"
  SET_EXPR = "(val > 0.0) ? val : ::sqrt (std::numeric_limits<double>::epsilon ())"
END_OPTION

OPTION
  NAME = "integration method"
  DOC_ITEM
A string specifying the method of integration to use to solve the ODE
system.  Valid values are

@table @asis
@item  @qcode{"adams"}
@itemx @qcode{"non-stiff"}
No Jacobian used (even if it is available).

@item  @qcode{"bdf"}
@itemx @qcode{"stiff"}
Use stiff backward differentiation formula (BDF) method.  If a
function to compute the Jacobian is not supplied, @code{lsode} will
compute a finite difference approximation of the Jacobian matrix.
@end table

  END_DOC_ITEM
  TYPE = "std::string"
  SET_ARG_TYPE = "const $TYPE&"
  INIT_VALUE = ""stiff""
  SET_BODY
    if (val == "stiff" || val == "bdf")
      $OPTVAR = "stiff";
    else if (val == "non-stiff" || val == "adams")
      $OPTVAR = "non-stiff";
    else
      (*current_liboctave_error_handler)
        ("lsode_options: method must be \"stiff\", \"bdf\", \"non-stiff\", or \"adams\"");
  END_SET_BODY
END_OPTION

OPTION
  NAME = "initial step size"
  DOC_ITEM
The step size to be attempted on the first step (default is determined
automatically).

  END_DOC_ITEM
  TYPE = "double"
  INIT_VALUE = "-1.0"
  SET_EXPR = "(val >= 0.0) ? val : -1.0"
END_OPTION

OPTION
  NAME = "maximum order"
  DOC_ITEM
Restrict the maximum order of the solution method.  If using the Adams
method, this option must be between 1 and 12.  Otherwise, it must be
between 1 and 5, inclusive.

  END_DOC_ITEM
  TYPE = "octave_idx_type"
  INIT_VALUE = "-1"
  SET_EXPR = "val"
END_OPTION

OPTION
  NAME = "maximum step size"
  DOC_ITEM
Setting the maximum stepsize will avoid passing over very large
regions  (default is not specified).

  END_DOC_ITEM
  TYPE = "double"
  INIT_VALUE = "-1.0"
  SET_EXPR = "(val >= 0.0) ? val : -1.0"
END_OPTION

OPTION
  NAME = "minimum step size"
  DOC_ITEM
The minimum absolute step size allowed (default is 0).

  END_DOC_ITEM
  TYPE = "double"
  INIT_VALUE = "0.0"
  SET_EXPR = "(val >= 0.0) ? val : 0.0"
END_OPTION

OPTION
  NAME = "step limit"
  DOC_ITEM
Maximum number of steps allowed (default is 100000).
  END_DOC_ITEM
  TYPE = "octave_idx_type"
  INIT_VALUE = "100000"
  SET_EXPR = "val"
END_OPTION

OPTION
  NAME = "jacobian type"
  DOC_ITEM
A string specifying the type of Jacobian used with the stiff backward
differentiation formula (BDF) integration method. Valid values are

@table @asis
@item @qcode{"full"}
The default. All partial derivatives are approximated or used from the
user-supplied Jacobian function.

@item @qcode{"banded"}
Only the diagonal and the number of lower and upper subdiagonals specified by
the options @qcode{"lower jacobian subdiagonals"} and @qcode{"upper jacobian
subdiagonals"}, respectively, are approximated or used from the user-supplied
Jacobian function. A user-supplied Jacobian function may set all other
partial derivatives to arbitrary values.

@item @qcode{"diagonal"}
If a Jacobian function is supplied by the user, this setting has no effect.
A Jacobian approximated by @code{lsode} is restricted to the diagonal, where
each partial derivative is computed by applying a finite change to all
elements of the state together; if the real Jacobian is indeed always diagonal,
this has the same effect as applying the finite change only to the respective
element of the state, but is more efficient.
@end table

  END_DOC_ITEM
  TYPE = "std::string"
  SET_ARG_TYPE = "const $TYPE&"
  INIT_VALUE = ""full""
  SET_BODY
    if (val == "full" || val == "banded" || val == "diagonal")
      $OPTVAR = val;
    else
      (*current_liboctave_error_handler)
        ("lsode_options: jacobian type must be \"full\", \"banded\", or \"diagonal\"");
  END_SET_BODY
END_OPTION

OPTION
  NAME = "lower jacobian subdiagonals"
  DOC_ITEM
Number of lower subdiagonals used if option @qcode{"jacobian type"} is set to
@qcode{"banded"}. The default is zero.

  END_DOC_ITEM
  TYPE = "octave_idx_type"
  INIT_VALUE = "0"
  SET_EXPR = "val"
END_OPTION

OPTION
  NAME = "upper jacobian subdiagonals"
  DOC_ITEM
Number of upper subdiagonals used if option @qcode{"jacobian type"} is set to
@qcode{"banded"}. The default is zero.

  END_DOC_ITEM
  TYPE = "octave_idx_type"
  INIT_VALUE = "0"
  SET_EXPR = "val"
END_OPTION