Module daskr

Integrate a system of ordinary differential algebraic equations with or
without events.

Description:

  Solve a system of ordinary differential algebraic equations using the
  FORTRAN library DDASKR.
  
  Uses the backward differentiation formulas of orders one through five
  to solve a system of the form G(T,Y,Y') = 0 where Y can be a vector.

  Values for Y and Y' at the initial time must be given as input.  These
  values should be consistent, that is, if t0, y0, yp0 are the given
  initial values, they should satisfy G(t0,y0,yp0) = 0.  However, if
  consistent values are not known, in many cases you can have DDASKR
  solve for them if the appropriate option in the info array is set.

  Normally, DDASKR solves the system from some intial t0 to a specified
  tout.  This wrapper takes a list of times t as input, and then solves
  the system for the range of times.  In the interval mode of operation,
  only the solution values at time points in t will be returned.
  Intermediate results can also be obtained easily by specifying info(3).

  While integrating the given DAE system, DDASKR also searches for roots
  of the given constraint functions Ri(T,Y,Y') given by RT. If DDASKR
  detects a sign change in any Ri(T,Y,Y'), it will return the intermediate
  value of T and Y for which Ri(T,Y,Y') = 0.

  Caution: If some Ri has a root at or very near the initial time, DDASKR
  may fail to find it, or may find extraneous roots there, because it does
  not yet have a sufficient history of the solution.
  
Inputs:

  res  -- res(t, y, ...) computes the residual function for the
          differential/algebraic system to be solved.
  t    -- a sequence of time points for which to solve for y.  The intial 
          value of the dependent variable(s) y should correspond to the
          first time in this sequence.
  y0   -- this array must contain the initial values of the NEQ solution
          components at the initial time point (note NEQ is defined as the
          length of this array).
  yp0  -- this array must contain the initial values of the NEQ first
          derivatives of the solution components at the initial time point.
  jac  -- this is the name of a routine that you may supply that relates to
          the Jacobian matrix of the nonlinear system that the code must
          solve at each time step.
  rtol -- array of relative error tolerances (must be non-negative).
  atol -- array of absolute error tolerances (must be non-negative).
  rt   -- this is the name of the subroutine for defining the vector
          Ri(T, Y, Y') of constraint functions.
  nrt  -- the number of constraint functions Ri(T, Y, Y').
  jac  -- the name of a routine that you may supply (optionally) that
          relactes to the Jacobian matrix of the nonlinear system that the
          code must solve at each time step. If no Jacobian is passed then
          the system will automatically use a numerical Jacobian.
  args -- parameter for passing extra information to the wrapper.
  tstop -- sometimes it is not possible to integrate past some point tstop
          because the equation changes there or it is not defined past
          tstop.  if tstop is set the integrator will not integrate past
          tstop.  NOTE:  If you select a tstop, the integrator may skip
          some points in your list of time points, t.  This is because it
          may pass tstop before hitting all the points in the list of times
          for an integration that is progressing quickly.
  intermediate_output -- True to have all output returned.
  ineq_constr -- True to have inequality constraint checking on.
          If True, the code will check for non-negativity in Y during the
          integration (not during initial condition calculations).
  calculate_ic -- True to have the DDASKR automatically calculate
         consistent initial conditions.  Given Y_d, the code will
         calculate Y_a and Y'_d (note: you must also pass var_types
         to indicate which variables are algebraic and which are
         differential).  We do not use the functionality of daskr
         that allows calculation of Y given Y'.
  var_types -- this list tells whether each variable in the system is
          differential or algebraic. Vor a variable Yi, var_types[i] = +1 if
          Yi is differential, and var_types[i] = -1 if it is algebraic.
  redir_output -- False to have print statements from DASKR printed to the
          standard output.  If redir_output is True, the output will be
          redirected.
  max_steps -- set this if you want the integrator to be able to take more
          than 500 steps between time points when you're not in
          "intermediate output" mode.  If intermediate_output = True then
          every time step is considered an output point so there are never
          multiple steps between time points.  If intermediate_output =
          False (default) then the integrator will stop after taking 500
          steps unless max_steps is set > 500.  The default value for
          max_steps is 500.
  rpar -- If not None, this sequence is passed to the function being
          evaluated and can be used to pass additional arguments.
  hmax -- If not None, this sets an upper limit to the step sizes that 
          ddaskr will take during the integration. Setting this can help
          integration for difficult systems, but it can also slow things
          down.
     
Outputs:  (yout, tout, ypout, t_root, y_root, i_root)

  yout -- the numerically calculated list of solution points corresponding
          to the times in tout.
  tout -- the list of time points corresponding to yout.
  ypout -- the numerically calculated list of solution points corresponding
  to the times in tout.
  t_root -- the time of a terminating event.
  y_root -- the system state at time t_root.
  i_root -- tells which event(s) fired and how they fired.
            If i_root(i) = 0, then event i did not fire. If nonzero,
            i_root(i) shows the direction of the sign change in Ri.
            i_root(i) = 1  means Ri changed from negative to positive.
            i_root(i) = -1 means Ri changed from positive to negative.