;+
; NAME:
;   DDEABM
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   craigm@lheamail.gsfc.nasa.gov
;   UPDATED VERSIONs can be found on my WEB PAGE: 
;      http://cow.physics.wisc.edu/~craigm/idl/idl.html
;
; PURPOSE:
;   Integrate Ordinary Differential Equation with Adams-Bashforth-Moulton
;
; MAJOR TOPICS:
;   Numerical Analysis.
;
; CALLING SEQUENCE:
;   DDEABM, FUNCT, T0, F0, TOUT, [ PRIVATE, FUNCTARGS=, STATE= , ]
;     [ /INIT, /INTERMEDIATE, TSTOP=, EPSREL=, EPSABS=, ]
;     [ TGRID=, YGRID=, YPGRID=, NOUTGRID=, NGRID=, NFEV=, ]
;     [ TIMPULSE=, YIMPULSE=, ]
;     [ MAX_STEPSIZE=, /CONTROL, ]
;     [ STATUS=, ERRMSG= ]
;
; DESCRIPTION:
;
;  DDEABM performs integration of a system of one or more ordinary
;  differential equations using a Predictor-Corrector technique.  An
;  adaptive Adams-Bashforth-Moulton method of variable order between
;  one and twelve, adaptive stepsize, and error control, is used to
;  integrate equations of the form:
;
;     DF_DT = FUNCT(T, F)
;
;  T is the independent variable, F is the (possibly vector) function
;  value at T, and DF_DT is the derivative of F with respect to T,
;  evaluated at T.  FUNCT is a user function which returns the
;  derivative of one or more equations.  
;
;  DDEABM is based on the public domain procedure DDEABM.F written by
;  L. F. Shampine and M. K. Gordon, and available in the DEPAC package
;  of solvers within SLATEC library.
;
;  DDEABM is used primarily to solve non-stiff and mildly stiff
;  ordinary differential equations, where evaluation of the user
;  function is expensive, or high precision is demanded (close to the
;  machine precision).  A Runge-Kutta technique may be more
;  appropriate if lower precision is acceptable.  For stiff problems
;  neither Adams-Bashforth-Moulton nor Runge-Kutta techniques are
;  appropriate.  DDEABM has code which detects the stiffness of the
;  problem and reports it.
;
;  The user can operate DDEABM in three different modes:
;
;    * the user requests one output at a specific point (the default),
;      and maintains the integrator state using the STATE keyword;
;
;    * the user requests a grid of points (by passing an array to
;      TOUT), and optionally maintains the integrator state for
;      subsequent integrations using the STATE keyword;
;
;    * the user requests output at the natural adaptive stepsizes
;      using the /INTERMEDIATE keyword.
;
;  When the user requests explicit output points using TOUT, it is
;  likely that the output will be interpolated between two natural
;  stepsize points.
;
;  If the user requests a grid of outputs, and the integration fails
;  before reaching the requested endpoint, then control returns
;  immediately to the user with the appropriate STATUS code.
;
;  The initial conditions are given by F0, when T = T0.  The number of
;  equations is determined by the number of elements in F0.
;  Integration stops when the independent variable T reaches the final
;  value of TOUT.  If the user wants to continue the integration, it
;  can be restarted with a new call to DDEABM, and a new value of
;  TOUT.
;
; USER FUNCTION
;
;  The user function FUNCT must be of the following form:
;
;  FUNCTION FUNCT, T, F, PRIVATE, ...  [ CONTROL=CONTROL ] ...
;     ; ... compute derivatives ...
;     RETURN, DF_DT
;  END
;
;  The function must accept at least two, but optionally three,
;  parameters.  The first, 'T', is the scalar independent variable
;  where the derivatives are to be evaluated.  The second, 'F', is the
;  vector of function values.  The function must return an array of
;  same size as 'F'.  The third positional parameter, 'PRIVATE', is a
;  purely optional PRIVATE parameter.  FUNCT may accept more
;  positional parameters, but DDEABM will not use them.  Any number of
;  keyword parameters can be passed using the FUNCTARGS keyword.
;
;  The user function FUNCT must not modify either the independent
;  variable 'T' or the function values 'F'.
;
; PASSING 'CONTROL' MESSAGES TO THE USER FUNCTION
;
;  DDEABM may pass control information to the user function, other
;  than requests for function evaluation.  DDEABM will only do this if
;  the /CONTROL keyword is set when DDEABM was invoked.
;
;  When control information has been enabled, the user function *must*
;  accept a keyword named CONTROL.  A message may be passed in this
;  keyword.  If the user function is invoked with the CONTROL keyword
;  defined, the user function should not evaluate the function, but
;  rather it must respond to the message and return the appropriate
;  value.
;
;  The CONTROL message will be a structure of the form,
;     { MESSAGE: 'name', ... }
;  where the MESSAGE field indicates a control message.  Additional
;  fields may carry more information, depending on the message.
;
;  The following control messages are implemented:
;    * { MESSAGE: 'INITIALIZE' } - the integrator has been initialized
;      and the user function may also initialize any state variables,
;      load large data tables, etc.
;      RETURN: 0 for success, negative number to indicate failure.
;
;  If the user function does not recognize a message, a value of 0
;  should be returned.
;
;
; RESTARTING THE INTEGRATOR
;
;  When restarting the integrator, the STATE keyword can be used to
;  save some computation time.  Upon return from DDEABM the STATE
;  keyword will contain a structure which describes the state of the
;  integrator at the output point.  The user need merely pass this
;  variable back into the next call to continue integration.  The
;  value of T, and the function value F, must not be adjusted between
;  calls to DDEABM.
;
;  If T or F0 are to be adjusted, then the integrator must be
;  restarted afresh *without* the previous state.  This can be
;  achieved either by reseting STATE to undefined, or by setting the
;  /INIT keyword to DDEABM.
;
; ERROR CONTROL
;
;  Local error control is governed by two keywords, EPSABS and EPSREL.
;  The former governs the absolute error, while the latter governs the
;  relative or fractional error.  The error test at each iteration is:
;
;     ABS(ERROR)   LE   EPSREL*ABS(F) + EPSABS
;  
;  A scalar value indicates the same constraint should be applied to
;  every equation; a vector array indicates constraints should be
;  applied individually to each component.
;
;  Setting EPSABS=0.D0 results in a pure relative error test on that
;  component. Setting EPSREL=0. results in a pure absolute error test
;  on that component.  A mixed test with non-zero EPSREL and EPSABS
;  corresponds roughly to a relative error test when the solution
;  component is much bigger than EPSABS and to an absolute error test
;  when the solution component is smaller than the threshold EPSABS.
;  
;  Proper selection of the absolute error control parameters EPSABS
;  requires you to have some idea of the scale of the solution
;  components.  To acquire this information may mean that you will
;  have to solve the problem more than once. In the absence of scale
;  information, you should ask for some relative accuracy in all the
;  components (by setting EPSREL values non-zero) and perhaps impose
;  extremely small absolute error tolerances to protect against the
;  danger of a solution component becoming zero.
;  
;  The code will not attempt to compute a solution at an accuracy
;  unreasonable for the machine being used.  It will advise you if you
;  ask for too much accuracy and inform you as to the maximum accuracy
;  it believes possible.
;
; HARD LIMIT ON T
;
;  If for some reason there is a hard limit on the independent
;  variable T, then the user should specify TSTOP.  For efficiency
;  DDEABM may try to integrate *beyond* the requested point of
;  termination, TOUT, and then interpolate backwards to obtain the
;  function value at the requested point.  If this is not possible
;  because the function because the equation changes, or if there is a
;  discontinuity, then users should specify a value for TSTOP; the
;  integrator will not go past this point.
;
; DISCONTINUITIES
;
;  If the ODE or solution has discontinuities at known points, these
;  should be passed to DDEABM in order to aid the solution.  The
;  TIMPULSE and YIMPULSE keyword variables allow the user to identify
;  the positions of the discontinuities and their amplitudes.  As T
;  crosses TIMPULSE(i) the solution will change from Y to
;  Y+YIMPULSE(*,i) in a stepwise fashion.  
;
;  Discontinuities in the function to be integrated can also be
;  entered in this way.  Although DDEABM can adapt the integration
;  step size to accomodate changes in the user function, it may be
;  better to identify such discontinuities.  In that case TIMPULSE(i)
;  should still identify the position of discontinuity, and
;  YIMPULSE(*,i) should be 0.
;
;  Currently this functionality is implemented with restarts of the
;  integrator at the crossing points of the discontinuities.
;
;  This technique handles only discontinuities at explicitly known
;  values of T.  If the discontinuity condition is defined in terms of
;  Y  (or Y and T), then the condition is implicit.  DDEABM does not
;  handle that type of condition.
;
;  You may list the TIMPULSE points in the TOUT output grid.  If an
;  impulse point appears once in TOUT, the corresponding function
;  values reported in YGRID and YPGRID will be from *before* crossing
;  the discontinuity.  If the same TIMPULSE point appears *twice* in
;  TOUT, then the first and second values will correspond to before
;  and after crossing the discontinuity, respectively.
;
;
; INPUTS:
;
;  FUNCT - by default, a scalar string containing the name of an IDL
;          function to be integrated.  See above for the formal
;          definition of FUNCT.  (No default).
;
;  T0 - scalar number, upon input the starting value of the
;       independent variable T.  Upon output, the final value of T.
;
;  Y - vector.  Upon input the starting values of the function for T =
;      T0.  Upon output, the final (vector) value of the function.
;
;  TOUT - must be at least a scalar, but optionally a vector,
;         specifies the desired points of output.
;
;         If TOUT is a scalar and INTERMEDIATE is not set, then DDEABM
;         integrates up to TOUT.  A scalar value of the function at
;         TOUT is returned in F.
;
;         If TOUT is a scalar and /INTERMEDIATE is set, then DDEABM
;         computes a grid of function values at the optimal step
;         points of the solution.  The grid of values is returned in
;         TGRID, YGRID, and YPGRID.  The final function value,
;         evaluated at the last point in TOUT, is returned in F.
;
;         If TOUT is a vector, then DDEABM computes a grid of function
;         values at the requested points.  The grid of values is
;         returned in TGRID, YGRID and YPGRID.  The final function
;         value, evaluated at the last point in TOUT, is returned in
;         F.  If integrating forwards (TOUT greater than T0), TOUT
;         must be strictly increasing.  Generally speaking, TOUT
;         output points should not repeat, except for discontinuities
;         as noted above.
;
;         It is possible for TOUT to be less than T0, i.e., to
;         integrate backwards, in which case TOUT must be strictly
;         decreasing instead.
; 
;  PRIVATE - any optional variable to be passed on to the function to
;            be integrated.  For functions, PRIVATE is passed as the
;            second positional parameter.  DDEABM does not examine or
;            alter PRIVATE.
;
; KEYWORD PARAMETERS:
;
;   CONTROL - if set, then control messages will be set to the user
;             function as described above.  If not set, then no
;             control messages will be passed.
;
;   EPSABS - a scalar number, the absolute error tolerance requested
;            in the integral computation.  If less than or equal to
;            zero, then the value is ignored.
;            Default: 0
;
;   EPSREL - a scalar number, the relative (i.e., fractional) error
;            tolerance in the integral computation.  If less than or
;            equal to zero, then the value is ignored.
;            Default: 1e-4 for float, 1d-6 for double
;
;   ERRMSG - upon return, a descriptive error message.
;
;   FUNCTARGS - A structure which contains the parameters to be passed
;               to the user-supplied function specified by FUNCT via
;               the _EXTRA mechanism.  This is the way you can pass
;               additional data to your user-supplied function without
;               using common blocks.  By default, no extra parameters
;               are passed to the user-supplied function.
;
;   INIT - if set, indicates that the integrator is to be restarted
;          afresh.
;
;   INTERMEDIATE - if set, indicates that the integrator is to compute
;                  at the natural step points.
;
;   MAX_STEPSIZE - a positive scalar value, the maximum integration
;                  step size to take per iteration.  The lesser of the
;                  "natural" step size and MAX_STEPSIZE is used.  If
;                  MAX_STEPSIZE is not specified, there is no maximum.
;
;   NFEV - upon output, the scalar number of function evaluations.
;
;   NGRID - if /INTERMEDIATE is set, the requested number of points to
;           compute before returning.  DDEABM uses this value to
;           allocate storage for TGRID, YGRID, and YPGRID.  Note that
;           DDEABM may not actually calculate this many points.  The
;           user must use NOUTGRID upon return to determine how many
;           points are valid.
;           Default: 1
;
;   NOUTGRID - upon output, the number of grid points computed.  This
;              may be smaller than the requested number of grid points
;              (either via NGRID or TOUT) if an error occurs.
;
;   STATE - upon input and return, the integrator state.  Users should
;           not modify this structure.  If the integrator is to be
;           restarted afresh, then the /INIT keyword should be set, in
;           order to clear out the old state information.
;
;   STATUS - upon output, the integer status of the integration.
;
;           1 - A step was successfully taken in the
;               intermediate-output mode.  The code has not yet
;               reached TOUT.
; 
;           2 - The integration to TOUT was successfully completed
;               (T=TOUT) by stepping exactly to TOUT.
; 
;           3 - The integration to TOUT was successfully completed
;               (T=TOUT) by stepping past TOUT.  Y(*) is obtained by
;               interpolation.
; 
;                       *** Task Interrupted ***
;                 Reported by negative values of STATUS
; 
;           -1 - A large amount of work has been expended.  (500 steps
;                attempted)
; 
;           -2 - The error tolerances are too stringent.
; 
;           -3 - The local error test cannot be satisfied because you
;                specified a zero component in EPSABS and the
;                corresponding computed solution component is zero.
;                Thus, a pure relative error test is impossible for
;                this component.
; 
;           -4 - The problem appears to be stiff.
; 
;                       *** Task Terminated ***
; 
;           -33 - The code has encountered trouble from which it
;                 cannot recover.  A error message is printed
;                 explaining the trouble and control is returned to
;                 the calling program. For example, this occurs when
;                 invalid input is detected.
;
;           -16 - The user function returned a non-finite 
;
;   TGRID - upon output, the grid of values of T for which the
;           integration is provided.  Upon return, only values
;           TGRID(0:NOUTGRID-1) are valid.  The remaining values are
;           undefined.
;
;   TIMPULSE - array of values of T where discontinuities occur.  The
;              array should be in ascending order.  TIMPULSE must
;              match YIMPULSE.
;
;   TSTOP - a scalar, specifies a hard limit for T, beyond which the
;           integration must not proceed.
;           Default:  none, i.e., integration is allowed to
;                    "overshoot"
;
;   YGRID - upon output, the grid of function values for which the
;           integration is provided.  Upon return, only values
;           YGRID(*,0:NOUTGRID-1) are valid.  The remaining values are
;           undefined.
;
;   YIMPULSE - array of discontinuity offset values, of the form
;              DBLARR(NEQ,NIMPULSE), where NEQ is the size of Y and
;              NIMPULSE is the size of TIMPULSE.  YIMPULSE(*,I) is the
;              amount to *add* to Y when T crosses TIMPULSE(I) going
;              in the positive direction.
;
;   YPGRID - upon ouptut, the grid of function derivative values at
;            the points where the integration is provided.  Upon
;            return, only values YPGRID(*,0:NOUTGRID-1) are valid.
;            The remaining values are undefined.
;
; EXAMPLES:
;
;  This is a simple example showing how to computes orbits of a
;  satellite around the earth using Newton's law of gravity.  The
;  earth is assumed to be a central point mass, modeled by the
;  NEWTON_G function which follows.  We assume that the satellite is
;  orbiting at a radius of 7000 km.  The state vector F has six
;  elements consisting of the position and velocity of the satellite.
;
;           POSITION    VELOCITY
;     F = [ X, Y, Z,    VX, VY, VZ]
;
;  The function NEWTON_G below computes the derivative of F, that is,
;
;               VELOCITY       ACCELERATION
;     dF_dt = [ VX, VY, VZ,    AX, AY, AZ]
;
;  Where the acceleration vector [AX,AY,AZ] is computed using Newton's
;  laws.
;
;  GM = 3.986005d14 ; [MKS] - gravitational constant for earth
;  
;  a0 = 7000d3       ; [m] - initial radius
;  v0 = sqrt(GM/a0)  ; [m/s] - initial circular velocity
;  ;     POSITION  VELOCITY
;  f0 = [a0,0,0,   0,-v0,0] ; initial state vector
;
;  t0 = 100d ; [s] Initial time value, meaningless in this case
;
;  ; Initial output time grid (10000 seconds)
;  tout = dindgen(10000) + t0
;  
;  ; Integrate equations of motion, starting at T0, and proceeding to
;  ; the maximum time of TOUT.  Here the variable GM is passed using
;  ; the PRIVATE mechanism.
;  f = f0 & t = t0
;  ddeabm, 'newton_g', t, f, tout, GM, $
;    tgrid=tgrid, ygrid=ygrid, ypgrid=ypgrid, noutgrid=noutgrid, $
;    status=status, errmsg=errmsg
;
;  Now YGRID(0:2,*)  contains the 3D position of the satellite
;      YGRID(3:5,*)  contains the 3D velocity of the satellite
;      YPGRID(3:5,*) contains the 3D acceleration of the satellite
;
;  An alternate way to call DDEABM is to use its natural gridpoints
;  rather than requesting explicit gridpoints.  In that case, we need
;  to specify the maximum time value we are expecing with TOUT, and
;  a maximum number of output grid values using NGRID.
;
;  f = f0 & t = t0
;  tout = 10000d ;; Maximum requested time
;  ddeabm, 'newton_g', t, f, tout, GM, $
;    ngrid=3000, noutgrid=noutgrid, $
;    tgrid=tgrid, ygrid=ygrid, ypgrid=ypgrid, noutgrid=noutgrid, $
;    status=status, errmsg=errmsg
;
;  NOUTGRID contains the actual number of grid values returned by 
;  DDEABM.  If NOUTGRID is less than NGRID, then the remaining values
;  are to be ignored.
;
;  TGRID = TGRID(0:NOUTGRID-1)
;  YGRID = YGRID(*,0:NOUTGRID-1)
;  YPGRID = YPGRID(*,0:NOUTGRID-1)
;
;  The user can then plot these values are use them as desired.  The
;  result should be a circular orbit at radius 7000000d meters, with
;  constant speed given by V0.
;
;
;  ; WORK FUNCTION ----
;
;  ; The acceleration of Newtonian gravity by a central body
;  ; of mass M.
;  ;   T - time (not used)
;  ;   f - state vector
;  ;     f(0:2) = position vector
;  ;     f(3:5) = velocity vector
;  ;   GM - central body newtonian constant
; FUNCTION NEWTON_G, t, f, GM
;    r = f(0:2) ; Position vector
;    v = f(3:5) ; Velocity vector
;    rsq = total(r^2,1) ;; central body distance, squared
;    rr = sqrt(rsq)     ;; central body distance
;   
;    ;; Newtonian gravitational acceleration, three components
;    a = - GM/rsq * r/rr
;
;    ;; Assemble final differential vector
;    df_dt = [v, a]
;    return, df_dt
; end
;
;
; REFERENCES:
;
;   SAND79-2374 , DEPAC - Design of a User Oriented Package of ODE
;                         Solvers.
;
;   "Solving Ordinary Differential Equations with ODE, STEP, and INTRP",
;       by L. F. Shampine and M. K. Gordon, SLA-73-1060.
;
;    SLATEC Common Mathematical Library, Version 4.1, July 1993
;    a comprehensive software library containing over
;    1400 general purpose mathematical and statistical routines
;    written in Fortran 77.  (http://www.netlib.org/slatec/)
;
;
; MODIFICATION HISTORY:
;    Fix bug in TSTOP keyword, 09 May 2002, CM
;    Fix behavior of KSTEPS, which caused premature termination, 26
;      May 2002, CM
;    Fix two errors in the DDEABM_DINTP interpolation step, 04 Jul 2002, CM
;    Handle case of IMPULSES more correctly, 25 Sep 2002, CM
;    Handle case when INIT is not set (default to 1); detect
;      non-finite user function values and error out with STATUS code
;      -16; promote integer values to LONG integers; some internal 
;      function renaming, 28 Jun 2003, CM
;    Fixed bug in handling of DOIMPULSE and INTERMEDIATE, 08 Mar 2004, CM
;    Corrected interface error in usage of NGRID.  Now NGRID is
;      actually the number of INTERMEDIATE points to compute (and is
;      input only).  NOUTGRID is a new keyword, which provides the
;      number of output grid points upon return.  08 Mar 2004, CM
;    Early termination is possible for INTERMEDIATE case.  Handle it
;      properly , 08 Mar 2004, CM
;    Fix a bug in the handling of INIT (strangely the internal
;      code keeps two different INIT variables!); this really only 
;      had an effect when continuing a previous integration; handle
;      impulses properly when integrating in the negative direction; 
;      document the TIMPULSE/YIMPULSE keyword parameters; some other
;      small code cleanups;  16 Jul 2008, CM
;    Handle the case when TOUT EQ TIMPULSE, 05 Sep 2008, CM
;    Further work on TOUT EQ TIMPULSE, also allowing reporting of
;      function values on either side of a discontinuity, 07 Sep 2008, CM
;    Add the MAX_STEPSIZE keyword, 01 Oct 2008, CM
;    Make sure new impulse checks work when integrating in reverse
;      direction, 09 Oct 2008, CM
;    New interface requirement: user function must be able to handle 
;      control messages from DDEABM; first one is INITIALIZE,
;      20 Oct 2008, CM
;    Change the control message interface so that it is
;      backward-compatible; the user must now set the /CONTROL keyword
;      to enable control messages; they are passed to the user
;      function via the CONTROL keyword, 08 Nov 2008, CM
;   Update the documentation; the largest change is the inclusion of
;      a new example, 16 Jan 2010, CM
;   Update documentation with correct spelling of Bashforth,
;      2012-12-18, CM
;
;  $Id: ddeabm.pro,v 1.33 2012/12/19 01:16:44 cmarkwar Exp $
;-
; Portions Copyright (C) 2002, 2003, 2004, 2008, 2010, 2012, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-

; *DECK DDEABM
; C***BEGIN PROLOGUE  DDEABM
; C***PURPOSE  Solve an initial value problem in ordinary differential
; C            equations using an Adams-Bashforth method.
; C***LIBRARY   SLATEC (DEPAC)
; C***CATEGORY  I1A1B
; C***TYPE      DOUBLE PRECISION (DEABM-S, DDEABM-D)
; C***KEYWORDS  ADAMS-BASHFORTH METHOD, DEPAC, INITIAL VALUE PROBLEMS,
; C             ODE, ORDINARY DIFFERENTIAL EQUATIONS, PREDICTOR-CORRECTOR
; C***AUTHOR  Shampine, L. F., (SNLA)
; C           Watts, H. A., (SNLA)
; C***DESCRIPTION
; C
; C   This is the Adams code in the package of differential equation
; C   solvers DEPAC, consisting of the codes DDERKF, DDEABM, and DDEBDF.
; C   Design of the package was by L. F. Shampine and H. A. Watts.
; C   It is documented in
; C        SAND79-2374 , DEPAC - Design of a User Oriented Package of ODE
; C                              Solvers.
; C   DDEABM is a driver for a modification of the code ODE written by
; C             L. F. Shampine and M. K. Gordon
; C             Sandia Laboratories
; C             Albuquerque, New Mexico 87185
;
;   $Id: ddeabm.pro,v 1.33 2012/12/19 01:16:44 cmarkwar Exp $
;
; C
; C **********************************************************************
; C * ABSTRACT *
; C ************
; C
; C   Subroutine DDEABM uses the Adams-Bashforth-Moulton
; C   Predictor-Corrector formulas of orders one through twelve to
; C   integrate a system of NEQ first order ordinary differential
; C   equations of the form
; C                         DU/DX = DF(X,U)
; C   when the vector Y(*) of initial values for U(*) at X=T is given.
; C   The subroutine integrates from T to TOUT. It is easy to continue the
; C   integration to get results at additional TOUT.  This is the interval
; C   mode of operation.  It is also easy for the routine to return with
; C   the solution at each intermediate step on the way to TOUT.  This is
; C   the intermediate-output mode of operation.
; C
; C   DDEABM uses subprograms DDEABM_DDES, DDEABM_DSTEPS, DDEABM_DINTP, DHSTRT, DHVNRM,
; C   D1MACH, and the error handling routine XERMSG.  The only machine
; C   dependent parameters to be assigned appear in D1MACH.
; C
; C **********************************************************************
; C * Description of The Arguments To DDEABM (An Overview) *
; C **********************************************************************
; C
; C   The Parameters are
; C
; C      DF -- This is the name of a subroutine which you provide to
; C             define the differential equations.
; C
; C      NEQ -- This is the number of (first order) differential
; C             equations to be integrated.
; C
; C      T -- This is a DOUBLE PRECISION value of the independent
; C           variable.
; C
; C      Y(*) -- This DOUBLE PRECISION array contains the solution
; C              components at T.
; C
; C      TOUT -- This is a DOUBLE PRECISION point at which a solution is
; C              desired.
; C
; C      INFO(*) -- The basic task of the code is to integrate the
; C             differential equations from T to TOUT and return an
; C             answer at TOUT.  INFO(*) is an INTEGER array which is used
; C             to communicate exactly how you want this task to be
; C             carried out.
; C
; C      RTOL, ATOL -- These DOUBLE PRECISION quantities represent
; C                    relative and absolute error tolerances which you
; C                    provide to indicate how accurately you wish the
; C                    solution to be computed.  You may choose them to be
; C                    both scalars or else both vectors.
; C
; C      IDID -- This scalar quantity is an indicator reporting what
; C             the code did.  You must monitor this INTEGER variable to
; C             decide what action to take next.
; C
; C      RWORK(*), LRW -- RWORK(*) is a DOUBLE PRECISION work array of
; C             length LRW which provides the code with needed storage
; C             space.
; C
; C      IWORK(*), LIW -- IWORK(*) is an INTEGER work array of length LIW
; C             which provides the code with needed storage space and an
; C             across call flag.
; C
; C      RPAR, IPAR -- These are DOUBLE PRECISION and INTEGER parameter
; C             arrays which you can use for communication between your
; C             calling program and the DF subroutine.
; C
; C  Quantities which are used as input items are
; C             NEQ, T, Y(*), TOUT, INFO(*),
; C             RTOL, ATOL, RWORK(1), LRW and LIW.
; C
; C  Quantities which may be altered by the code are
; C             T, Y(*), INFO(1), RTOL, ATOL,
; C             IDID, RWORK(*) and IWORK(*).
; C
; C **********************************************************************
; C * INPUT -- What To Do On The First Call To DDEABM *
; C **********************************************************************
; C
; C   The first call of the code is defined to be the start of each new
; C   problem.  Read through the descriptions of all the following items,
; C   provide sufficient storage space for designated arrays, set
; C   appropriate variables for the initialization of the problem, and
; C   give information about how you want the problem to be solved.
; C
; C
; C      DF -- Provide a subroutine of the form
; C                               DF(X,U,UPRIME,PAR,IPAR)
; C             to define the system of first order differential equations
; C             which is to be solved.  For the given values of X and the
; C             vector  U(*)=(U(1),U(2),...,U(NEQ)) , the subroutine must
; C             evaluate the NEQ components of the system of differential
; C             equations  DU/DX=DF(X,U)  and store the derivatives in the
; C             array UPRIME(*), that is,  UPRIME(I) = * DU(I)/DX *  for
; C             equations I=1,...,NEQ.
; C
; C             Subroutine DF must NOT alter X or U(*).  You must declare
; C             the name df in an external statement in your program that
; C             calls DDEABM.  You must dimension U and UPRIME in DF.
; C
; C             RPAR and IPAR are DOUBLE PRECISION and INTEGER parameter
; C             arrays which you can use for communication between your
; C             calling program and subroutine DF. They are not used or
; C             altered by DDEABM.  If you do not need RPAR or IPAR,
; C             ignore these parameters by treating them as dummy
; C             arguments. If you do choose to use them, dimension them in
; C             your calling program and in DF as arrays of appropriate
; C             length.
; C
; C      NEQ -- Set it to the number of differential equations.
; C             (NEQ .GE. 1)
; C
; C      T -- Set it to the initial point of the integration.
; C             You must use a program variable for T because the code
; C             changes its value.
; C
; C      Y(*) -- Set this vector to the initial values of the NEQ solution
; C             components at the initial point.  You must dimension Y at
; C             least NEQ in your calling program.
; C
; C      TOUT -- Set it to the first point at which a solution
; C             is desired.  You can take TOUT = T, in which case the code
; C             will evaluate the derivative of the solution at T and
; C             return. Integration either forward in T  (TOUT .GT. T)  or
; C             backward in T  (TOUT .LT. T)  is permitted.
; C
; C             The code advances the solution from T to TOUT using
; C             step sizes which are automatically selected so as to
; C             achieve the desired accuracy.  If you wish, the code will
; C             return with the solution and its derivative following
; C             each intermediate step (intermediate-output mode) so that
; C             you can monitor them, but you still must provide TOUT in
; C             accord with the basic aim of the code.
; C
; C             The first step taken by the code is a critical one
; C             because it must reflect how fast the solution changes near
; C             the initial point.  The code automatically selects an
; C             initial step size which is practically always suitable for
; C             the problem. By using the fact that the code will not step
; C             past TOUT in the first step, you could, if necessary,
; C             restrict the length of the initial step size.
; C
; C             For some problems it may not be permissible to integrate
; C             past a point TSTOP because a discontinuity occurs there
; C             or the solution or its derivative is not defined beyond
; C             TSTOP.  When you have declared a TSTOP point (see INFO(4)
; C             and RWORK(1)), you have told the code not to integrate
; C             past TSTOP.  In this case any TOUT beyond TSTOP is invalid
; C             input.
; C
; C      INFO(*) -- Use the INFO array to give the code more details about
; C             how you want your problem solved.  This array should be
; C             dimensioned of length 15 to accommodate other members of
; C             DEPAC or possible future extensions, though DDEABM uses
; C             only the first four entries.  You must respond to all of
; C             the following items which are arranged as questions.  The
; C             simplest use of the code corresponds to answering all
; C             questions as YES ,i.e. setting ALL entries of INFO to 0.
; C
; C        INFO(1) -- This parameter enables the code to initialize
; C               itself.  You must set it to indicate the start of every
; C               new problem.
; C
; C            **** Is this the first call for this problem ...
; C                  YES -- set INFO(1) = 0
; C                   NO -- not applicable here.
; C                         See below for continuation calls.  ****
; C
; C        INFO(2) -- How much accuracy you want of your solution
; C               is specified by the error tolerances RTOL and ATOL.
; C               The simplest use is to take them both to be scalars.
; C               To obtain more flexibility, they can both be vectors.
; C               The code must be told your choice.
; C
; C            **** Are both error tolerances RTOL, ATOL scalars ...
; C                  YES -- set INFO(2) = 0
; C                         and input scalars for both RTOL and ATOL
; C                   NO -- set INFO(2) = 1
; C                         and input arrays for both RTOL and ATOL ****
; C
; C        INFO(3) -- The code integrates from T in the direction
; C               of TOUT by steps.  If you wish, it will return the
; C               computed solution and derivative at the next
; C               intermediate step (the intermediate-output mode) or
; C               TOUT, whichever comes first.  This is a good way to
; C               proceed if you want to see the behavior of the solution.
; C               If you must have solutions at a great many specific
; C               TOUT points, this code will compute them efficiently.
; C
; C            **** Do you want the solution only at
; C                 TOUT (and not at the next intermediate step) ...
; C                  YES -- set INFO(3) = 0
; C                   NO -- set INFO(3) = 1 ****
; C
; C        INFO(4) -- To handle solutions at a great many specific
; C               values TOUT efficiently, this code may integrate past
; C               TOUT and interpolate to obtain the result at TOUT.
; C               Sometimes it is not possible to integrate beyond some
; C               point TSTOP because the equation changes there or it is
; C               not defined past TSTOP.  Then you must tell the code
; C               not to go past.
; C
; C            **** Can the integration be carried out without any
; C                 Restrictions on the independent variable T ...
; C                  YES -- set INFO(4)=0
; C                   NO -- set INFO(4)=1
; C                         and define the stopping point TSTOP by
; C                         setting RWORK(1)=TSTOP ****
; C
; C      RTOL, ATOL -- You must assign relative (RTOL) and absolute (ATOL)
; C             error tolerances to tell the code how accurately you want
; C             the solution to be computed.  They must be defined as
; C             program variables because the code may change them.  You
; C             have two choices --
; C                  Both RTOL and ATOL are scalars. (INFO(2)=0)
; C                  Both RTOL and ATOL are vectors. (INFO(2)=1)
; C             In either case all components must be non-negative.
; C
; C             The tolerances are used by the code in a local error test
; C             at each step which requires roughly that
; C                     ABS(LOCAL ERROR) .LE. RTOL*ABS(Y)+ATOL
; C             for each vector component.
; C             (More specifically, a Euclidean norm is used to measure
; C             the size of vectors, and the error test uses the magnitude
; C             of the solution at the beginning of the step.)
; C
; C             The true (global) error is the difference between the true
; C             solution of the initial value problem and the computed
; C             approximation.  Practically all present day codes,
; C             including this one, control the local error at each step
; C             and do not even attempt to control the global error
; C             directly.  Roughly speaking, they produce a solution Y(T)
; C             which satisfies the differential equations with a
; C             residual R(T),    DY(T)/DT = DF(T,Y(T)) + R(T)   ,
; C             and, almost always, R(T) is bounded by the error
; C             tolerances.  Usually, but not always, the true accuracy of
; C             the computed Y is comparable to the error tolerances. This
; C             code will usually, but not always, deliver a more accurate
; C             solution if you reduce the tolerances and integrate again.
; C             By comparing two such solutions you can get a fairly
; C             reliable idea of the true error in the solution at the
; C             bigger tolerances.
; C
; C             Setting ATOL=0.D0 results in a pure relative error test on
; C             that component. Setting RTOL=0. results in a pure absolute
; C             error test on that component.  A mixed test with non-zero
; C             RTOL and ATOL corresponds roughly to a relative error
; C             test when the solution component is much bigger than ATOL
; C             and to an absolute error test when the solution component
; C             is smaller than the threshold ATOL.
; C
; C             Proper selection of the absolute error control parameters
; C             ATOL  requires you to have some idea of the scale of the
; C             solution components.  To acquire this information may mean
; C             that you will have to solve the problem more than once. In
; C             the absence of scale information, you should ask for some
; C             relative accuracy in all the components (by setting  RTOL
; C             values non-zero) and perhaps impose extremely small
; C             absolute error tolerances to protect against the danger of
; C             a solution component becoming zero.
; C
; C             The code will not attempt to compute a solution at an
; C             accuracy unreasonable for the machine being used.  It will
; C             advise you if you ask for too much accuracy and inform
; C             you as to the maximum accuracy it believes possible.
; C
; C      RWORK(*) -- Dimension this DOUBLE PRECISION work array of length
; C             LRW in your calling program.
; C
; C      RWORK(1) -- If you have set INFO(4)=0, you can ignore this
; C             optional input parameter.  Otherwise you must define a
; C             stopping point TSTOP by setting   RWORK(1) = TSTOP.
; C             (for some problems it may not be permissible to integrate
; C             past a point TSTOP because a discontinuity occurs there
; C             or the solution or its derivative is not defined beyond
; C             TSTOP.)
; C
; C      LRW -- Set it to the declared length of the RWORK array.
; C             You must have  LRW .GE. 130+21*NEQ
; C
; C      IWORK(*) -- Dimension this INTEGER work array of length LIW in
; C             your calling program.
; C
; C      LIW -- Set it to the declared length of the IWORK array.
; C             You must have  LIW .GE. 51
; C
; C      RPAR, IPAR -- These are parameter arrays, of DOUBLE PRECISION and
; C             INTEGER type, respectively.  You can use them for
; C             communication between your program that calls DDEABM and
; C             the  DF subroutine.  They are not used or altered by
; C             DDEABM.  If you do not need RPAR or IPAR, ignore these
; C             parameters by treating them as dummy arguments.  If you do
; C             choose to use them, dimension them in your calling program
; C             and in DF as arrays of appropriate length.
; C
; C **********************************************************************
; C * OUTPUT -- After Any Return From DDEABM *
; C **********************************************************************
; C
; C   The principal aim of the code is to return a computed solution at
; C   TOUT, although it is also possible to obtain intermediate results
; C   along the way.  To find out whether the code achieved its goal
; C   or if the integration process was interrupted before the task was
; C   completed, you must check the IDID parameter.
; C
; C
; C      T -- The solution was successfully advanced to the
; C             output value of T.
; C
; C      Y(*) -- Contains the computed solution approximation at T.
; C             You may also be interested in the approximate derivative
; C             of the solution at T.  It is contained in
; C             RWORK(21),...,RWORK(20+NEQ).
; C
; C      IDID -- Reports what the code did
; C
; C                         *** Task Completed ***
; C                   Reported by positive values of IDID
; C
; C             IDID = 1 -- A step was successfully taken in the
; C                       intermediate-output mode.  The code has not
; C                       yet reached TOUT.
; C
; C             IDID = 2 -- The integration to TOUT was successfully
; C                       completed (T=TOUT) by stepping exactly to TOUT.
; C
; C             IDID = 3 -- The integration to TOUT was successfully
; C                       completed (T=TOUT) by stepping past TOUT.
; C                       Y(*) is obtained by interpolation.
; C
; C                         *** Task Interrupted ***
; C                   Reported by negative values of IDID
; C
; C             IDID = -1 -- A large amount of work has been expended.
; C                       (500 steps attempted)
; C
; C             IDID = -2 -- The error tolerances are too stringent.
; C
; C             IDID = -3 -- The local error test cannot be satisfied
; C                       because you specified a zero component in ATOL
; C                       and the corresponding computed solution
; C                       component is zero.  Thus, a pure relative error
; C                       test is impossible for this component.
; C
; C             IDID = -4 -- The problem appears to be stiff.
; C
; C             IDID = -5,-6,-7,..,-32  -- Not applicable for this code
; C                       but used by other members of DEPAC or possible
; C                       future extensions.
; C
; C                         *** Task Terminated ***
; C                   Reported by the value of IDID=-33
; C
; C             IDID = -33 -- The code has encountered trouble from which
; C                       it cannot recover.  A message is printed
; C                       explaining the trouble and control is returned
; C                       to the calling program. For example, this occurs
; C                       when invalid input is detected.
; C
; C      RTOL, ATOL -- These quantities remain unchanged except when
; C             IDID = -2.  In this case, the error tolerances have been
; C             increased by the code to values which are estimated to be
; C             appropriate for continuing the integration.  However, the
; C             reported solution at T was obtained using the input values
; C             of RTOL and ATOL.
; C
; C      RWORK, IWORK -- Contain information which is usually of no
; C             interest to the user but necessary for subsequent calls.
; C             However, you may find use for
; C
; C             RWORK(11)--which contains the step size H to be
; C                        attempted on the next step.
; C
; C             RWORK(12)--if the tolerances have been increased by the
; C                        code (IDID = -2) , they were multiplied by the
; C                        value in RWORK(12).
; C
; C             RWORK(13)--Which contains the current value of the
; C                        independent variable, i.e. the farthest point
; C                        integration has reached. This will be different
; C                        from T only when interpolation has been
; C                        performed (IDID=3).
; C
; C             RWORK(20+I)--Which contains the approximate derivative
; C                        of the solution component Y(I).  In DDEABM, it
; C                        is obtained by calling subroutine DF to
; C                        evaluate the differential equation using T and
; C                        Y(*) when IDID=1 or 2, and by interpolation
; C                        when IDID=3.
; C
; C **********************************************************************
; C * INPUT -- What To Do To Continue The Integration *
; C *             (calls after the first)             *
; C **********************************************************************
; C
; C        This code is organized so that subsequent calls to continue the
; C        integration involve little (if any) additional effort on your
; C        part. You must monitor the IDID parameter in order to determine
; C        what to do next.
; C
; C        Recalling that the principal task of the code is to integrate
; C        from T to TOUT (the interval mode), usually all you will need
; C        to do is specify a new TOUT upon reaching the current TOUT.
; C
; C        Do not alter any quantity not specifically permitted below,
; C        in particular do not alter NEQ, T, Y(*), RWORK(*), IWORK(*) or
; C        the differential equation in subroutine DF. Any such alteration
; C        constitutes a new problem and must be treated as such, i.e.
; C        you must start afresh.
; C
; C        You cannot change from vector to scalar error control or vice
; C        versa (INFO(2)) but you can change the size of the entries of
; C        RTOL, ATOL.  Increasing a tolerance makes the equation easier
; C        to integrate.  Decreasing a tolerance will make the equation
; C        harder to integrate and should generally be avoided.
; C
; C        You can switch from the intermediate-output mode to the
; C        interval mode (INFO(3)) or vice versa at any time.
; C
; C        If it has been necessary to prevent the integration from going
; C        past a point TSTOP (INFO(4), RWORK(1)), keep in mind that the
; C        code will not integrate to any TOUT beyond the currently
; C        specified TSTOP.  Once TSTOP has been reached you must change
; C        the value of TSTOP or set INFO(4)=0.  You may change INFO(4)
; C        or TSTOP at any time but you must supply the value of TSTOP in
; C        RWORK(1) whenever you set INFO(4)=1.
; C
; C        The parameter INFO(1) is used by the code to indicate the
; C        beginning of a new problem and to indicate whether integration
; C        is to be continued.  You must input the value  INFO(1) = 0
; C        when starting a new problem.  You must input the value
; C        INFO(1) = 1  if you wish to continue after an interrupted task.
; C        Do not set  INFO(1) = 0  on a continuation call unless you
; C        want the code to restart at the current T.
; C
; C                         *** Following A Completed Task ***
; C         If
; C             IDID = 1, call the code again to continue the integration
; C                     another step in the direction of TOUT.
; C
; C             IDID = 2 or 3, define a new TOUT and call the code again.
; C                     TOUT must be different from T. You cannot change
; C                     the direction of integration without restarting.
; C
; C                         *** Following An Interrupted Task ***
; C                     To show the code that you realize the task was
; C                     interrupted and that you want to continue, you
; C                     must take appropriate action and reset INFO(1) = 1
; C         If
; C             IDID = -1, the code has attempted 500 steps.
; C                     If you want to continue, set INFO(1) = 1 and
; C                     call the code again. An additional 500 steps
; C                     will be allowed.
; C
; C             IDID = -2, the error tolerances RTOL, ATOL have been
; C                     increased to values the code estimates appropriate
; C                     for continuing.  You may want to change them
; C                     yourself.  If you are sure you want to continue
; C                     with relaxed error tolerances, set INFO(1)=1 and
; C                     call the code again.
; C
; C             IDID = -3, a solution component is zero and you set the
; C                     corresponding component of ATOL to zero.  If you
; C                     are sure you want to continue, you must first
; C                     alter the error criterion to use positive values
; C                     for those components of ATOL corresponding to zero
; C                     solution components, then set INFO(1)=1 and call
; C                     the code again.
; C
; C             IDID = -4, the problem appears to be stiff.  It is very
; C                     inefficient to solve such problems with DDEABM.
; C                     The code DDEBDF in DEPAC handles this task
; C                     efficiently.  If you are absolutely sure you want
; C                     to continue with DDEABM, set INFO(1)=1 and call
; C                     the code again.
; C
; C             IDID = -5,-6,-7,..,-32  --- cannot occur with this code
; C                     but used by other members of DEPAC or possible
; C                     future extensions.
; C
; C                         *** Following A Terminated Task ***
; C         If
; C             IDID = -33, you cannot continue the solution of this
; C                     problem.  An attempt to do so will result in your
; C                     run being terminated.
; C
; C **********************************************************************
; C *Long Description:
; C
; C **********************************************************************
; C *             DEPAC Package Overview           *
; C **********************************************************************
; C
; C ....   You have a choice of three differential equation solvers from
; C ....   DEPAC. The following brief descriptions are meant to aid you in
; C ....   choosing the most appropriate code for your problem.
; C
; C ....   DDERKF is a fifth order Runge-Kutta code. It is the simplest of
; C ....   the three choices, both algorithmically and in the use of the
; C ....   code. DDERKF is primarily designed to solve non-stiff and
; C ....   mildly stiff differential equations when derivative evaluations
; C ....   are not expensive. It should generally not be used to get high
; C ....   accuracy results nor answers at a great many specific points.
; C ....   Because DDERKF has very low overhead costs, it will usually
; C ....   result in the least expensive integration when solving
; C ....   problems requiring a modest amount of accuracy and having
; C ....   equations that are not costly to evaluate. DDERKF attempts to
; C ....   discover when it is not suitable for the task posed.
; C
; C ....   DDEABM is a variable order (one through twelve) Adams code.
; C ....   Its complexity lies somewhere between that of DDERKF and
; C ....   DDEBDF.  DDEABM is primarily designed to solve non-stiff and
; C ....   mildly stiff differential equations when derivative evaluations
; C ....   are expensive, high accuracy results are needed or answers at
; C ....   many specific points are required. DDEABM attempts to discover
; C ....   when it is not suitable for the task posed.
; C
; C ....   DDEBDF is a variable order (one through five) backward
; C ....   differentiation formula code. it is the most complicated of
; C ....   the three choices. DDEBDF is primarily designed to solve stiff
; C ....   differential equations at crude to moderate tolerances.
; C ....   If the problem is very stiff at all, DDERKF and DDEABM will be
; C ....   quite inefficient compared to DDEBDF. However, DDEBDF will be
; C ....   inefficient compared to DDERKF and DDEABM on non-stiff problems
; C ....   because it uses much more storage, has a much larger overhead,
; C ....   and the low order formulas will not give high accuracies
; C ....   efficiently.
; C
; C ....   The concept of stiffness cannot be described in a few words.
; C ....   If you do not know the problem to be stiff, try either DDERKF
; C ....   or DDEABM. Both of these codes will inform you of stiffness
; C ....   when the cost of solving such problems becomes important.
; C
; C *********************************************************************
; C
; C***REFERENCES  L. F. Shampine and H. A. Watts, DEPAC - design of a user
; C                 oriented package of ODE solvers, Report SAND79-2374,
; C                 Sandia Laboratories, 1979.
; C***ROUTINES CALLED  DDEABM_DDES, XERMSG
; C***REVISION HISTORY  (YYMMDD)
; C   820301  DATE WRITTEN
; C   890531  Changed all specific intrinsics to generic.  (WRB)
; C   890831  Modified array declarations.  (WRB)
; C   891006  Cosmetic changes to prologue.  (WRB)
; C   891024  Changed references from DVNORM to DHVNRM.  (WRB)
; C   891024  REVISION DATE from Version 3.2
; C   891214  Prologue converted to Version 4.0 format.  (BAB)
; C   900510  Convert XERRWV calls to XERMSG calls.  (RWC)
; C   920501  Reformatted the REFERENCES section.  (WRB)
; C***END PROLOGUE  DDEABM
; C
;       INTEGER IALPHA, IBETA, IDELSN, IDID, IFOURU, IG, IHOLD,
;      1      INFO, IP, IPAR, IPHI, IPSI, ISIG, ITOLD, ITSTAR, ITWOU,
;      2      IV, IW, IWORK, IWT, IYP, IYPOUT, IYY, LIW, LRW, NEQ
;       DOUBLE PRECISION ATOL, RPAR, RTOL, RWORK, T, TOUT, Y
;       LOGICAL START,PHASE1,NORND,STIFF,INTOUT
; C
;       DIMENSION Y(*),INFO(15),RTOL(*),ATOL(*),RWORK(*),IWORK(*),
;      1          RPAR(*),IPAR(*)
; C
;       CHARACTER*8 XERN1
;       CHARACTER*16 XERN3
; C
;       EXTERNAL DF
; C
; C     CHECK FOR AN APPARENT INFINITE LOOP
; C

;; -----------  SOI start of IDL code ----------

pro ddeabm_dummy
  common ddeabm_func_common, ddeabm_nfev, ddeabm_funcerror
end

function ddeabm_func0n, func, a, y, private, _extra=fa
  common ddeabm_func_common, nfev, error
  nfev = nfev + 1
  dydx = call_function(func, a, y)
  if min(finite(dydx)) EQ 0 then error = -16
  return, dydx
end

function ddeabm_func1n, func, a, y, private, _extra=fa
  common ddeabm_func_common, nfev, error
  nfev = nfev + 1
  dydx = call_function(func, a, y, private)
  if min(finite(dydx)) EQ 0 then error = -16
  return, dydx
end

function ddeabm_func0e, func, a, y, private, _extra=fa
  common ddeabm_func_common, nfev, error
  nfev = nfev + 1
  dydx = call_function(func, a, y, _extra=fa)
  if min(finite(dydx)) EQ 0 then error = -16
  return, dydx
end

function ddeabm_func1e, func, a, y, private, _extra=fa
  common ddeabm_func_common, nfev, error
  nfev = nfev + 1
  dydx = call_function(func, a, y, private, _extra=fa)
  if min(finite(dydx)) EQ 0 then error = -16
  return, dydx
end

; *DECK DDEABM_DHSTRT
pro ddeabm_dhstrt, DF, NEQ, A, B, Y, YPRIME, ETOL, MORDER, SMALL, $
          BIG, SPY, PV, YP, SF, PRIVATE, FA, H, DFNAME
; C***BEGIN PROLOGUE  DDEABM_DHSTRT
; C***SUBSIDIARY
; C***PURPOSE  Subsidiary to DDEABM, DDEBDF and DDERKF
; C***LIBRARY   SLATEC
; C***TYPE      DOUBLE PRECISION (HSTART-S, DHSTRT-D)
; C***AUTHOR  Watts, H. A., (SNLA)
; C***DESCRIPTION
; C
; C   DDEABM_DHSTRT computes a starting step size to be used in solving initial
; C   value problems in ordinary differential equations.
; C
; C **********************************************************************
; C  ABSTRACT
; C
; C     Subroutine DDEABM_DHSTRT computes a starting step size to be used by an
; C     initial value method in solving ordinary differential equations.
; C     It is based on an estimate of the local Lipschitz constant for the
; C     differential equation   (lower bound on a norm of the Jacobian) ,
; C     a bound on the differential equation  (first derivative) , and
; C     a bound on the partial derivative of the equation with respect to
; C     the independent variable.
; C     (all approximated near the initial point A)
; C
; C     Subroutine DDEABM_DHSTRT uses a function subprogram DHVNRM for computing
; C     a vector norm. The maximum norm is presently utilized though it
; C     can easily be replaced by any other vector norm. It is presumed
; C     that any replacement norm routine would be carefully coded to
; C     prevent unnecessary underflows or overflows from occurring, and
; C     also, would not alter the vector or number of components.
; C
; C **********************************************************************
; C  On input you must provide the following
; C
; C      DF -- This is a subroutine of the form
; C                               DF(X,U,UPRIME,RPAR,IPAR)
; C             which defines the system of first order differential
; C             equations to be solved. For the given values of X and the
; C             vector  U(*)=(U(1),U(2),...,U(NEQ)) , the subroutine must
; C             evaluate the NEQ components of the system of differential
; C             equations  DU/DX=DF(X,U)  and store the derivatives in the
; C             array UPRIME(*), that is,  UPRIME(I) = * DU(I)/DX *  for
; C             equations I=1,...,NEQ.
; C
; C             Subroutine DF must not alter X or U(*). You must declare
; C             the name DF in an external statement in your program that
; C             calls DDEABM_DHSTRT. You must dimension U and UPRIME in DF.
; C
; C             RPAR and IPAR are DOUBLE PRECISION and INTEGER parameter
; C             arrays which you can use for communication between your
; C             program and subroutine DF. They are not used or altered by
; C             DDEABM_DHSTRT. If you do not need RPAR or IPAR, ignore these
; C             parameters by treating them as dummy arguments. If you do
; C             choose to use them, dimension them in your program and in
; C             DF as arrays of appropriate length.
; C
; C      NEQ -- This is the number of (first order) differential equations
; C             to be integrated.
; C
; C      A -- This is the initial point of integration.
; C
; C      B -- This is a value of the independent variable used to define
; C             the direction of integration. A reasonable choice is to
; C             set  B  to the first point at which a solution is desired.
; C             You can also use  B, if necessary, to restrict the length
; C             of the first integration step because the algorithm will
; C             not compute a starting step length which is bigger than
; C             ABS(B-A), unless  B  has been chosen too close to  A.
; C             (it is presumed that DDEABM_DHSTRT has been called with  B
; C             different from  A  on the machine being used. Also see the
; C             discussion about the parameter  SMALL.)
; C
; C      Y(*) -- This is the vector of initial values of the NEQ solution
; C             components at the initial point  A.
; C
; C      YPRIME(*) -- This is the vector of derivatives of the NEQ
; C             solution components at the initial point  A.
; C             (defined by the differential equations in subroutine DF)
; C
; C      ETOL -- This is the vector of error tolerances corresponding to
; C             the NEQ solution components. It is assumed that all
; C             elements are positive. Following the first integration
; C             step, the tolerances are expected to be used by the
; C             integrator in an error test which roughly requires that
; C                        ABS(LOCAL ERROR)  .LE.  ETOL
; C             for each vector component.
; C
; C      MORDER -- This is the order of the formula which will be used by
; C             the initial value method for taking the first integration
; C             step.
; C
; C      SMALL -- This is a small positive machine dependent constant
; C             which is used for protecting against computations with
; C             numbers which are too small relative to the precision of
; C             floating point arithmetic.  SMALL  should be set to
; C             (approximately) the smallest positive DOUBLE PRECISION
; C             number such that  (1.+SMALL) .GT. 1.  on the machine being
; C             used. The quantity  SMALL**(3/8)  is used in computing
; C             increments of variables for approximating derivatives by
; C             differences.  Also the algorithm will not compute a
; C             starting step length which is smaller than
; C             100*SMALL*ABS(A).
; C
; C      BIG -- This is a large positive machine dependent constant which
; C             is used for preventing machine overflows. A reasonable
; C             choice is to set big to (approximately) the square root of
; C             the largest DOUBLE PRECISION number which can be held in
; C             the machine.
; C
; C      SPY(*),PV(*),YP(*),SF(*) -- These are DOUBLE PRECISION work
; C             arrays of length NEQ which provide the routine with needed
; C             storage space.
; C
; C      RPAR,IPAR -- These are parameter arrays, of DOUBLE PRECISION and
; C             INTEGER type, respectively, which can be used for
; C             communication between your program and the DF subroutine.
; C             They are not used or altered by DDEABM_DHSTRT.
; C
; C **********************************************************************
; C  On Output  (after the return from DDEABM_DHSTRT),
; C
; C      H -- is an appropriate starting step size to be attempted by the
; C             differential equation method.
; C
; C           All parameters in the call list remain unchanged except for
; C           the working arrays SPY(*),PV(*),YP(*), and SF(*).
; C
; C **********************************************************************
; C
; C***SEE ALSO  DDEABM, DDEBDF, DDERKF
; C***ROUTINES CALLED  DHVNRM
; C***REVISION HISTORY  (YYMMDD)
; C   820301  DATE WRITTEN
; C   890531  Changed all specific intrinsics to generic.  (WRB)
; C   890831  Modified array declarations.  (WRB)
; C   890911  Removed unnecessary intrinsics.  (WRB)
; C   891024  Changed references from DVNORM to DHVNRM.  (WRB)
; C   891214  Prologue converted to Version 4.0 format.  (BAB)
; C   900328  Added TYPE section.  (WRB)
; C   910722  Updated AUTHOR section.  (ALS)
; C***END PROLOGUE  DDEABM_DHSTRT
; C
;       INTEGER IPAR, J, K, LK, MORDER, NEQ
;       DOUBLE PRECISION A, ABSDX, B, BIG, DA, DELF, DELY,
;      1      DFDUB, DFDXB, DHVNRM,
;      2      DX, DY, ETOL, FBND, H, PV, RELPER, RPAR, SF, SMALL, SPY,
;      3      SRYDPB, TOLEXP, TOLMIN, TOLP, TOLSUM, Y, YDPB, YP, YPRIME
;       DIMENSION Y(*),YPRIME(*),ETOL(*),SPY(*),PV(*),YP(*),
;      1          SF(*),RPAR(*),IPAR(*)
;       EXTERNAL DF
; C
; C     ..................................................................
; C
; C     BEGIN BLOCK PERMITTING ...EXITS TO 160
; C***FIRST EXECUTABLE STATEMENT  DDEABM_DHSTRT
         common ddeabm_func_common
         DX = B - A
         ABSDX = ABS(DX)
         RELPER = SMALL^0.375D0
; C
; C        ...............................................................
; C
; C             COMPUTE AN APPROXIMATE BOUND (DFDXB) ON THE PARTIAL
; C             DERIVATIVE OF THE EQUATION WITH RESPECT TO THE
; C             INDEPENDENT VARIABLE. PROTECT AGAINST AN OVERFLOW.
; C             ALSO COMPUTE A BOUND (FBND) ON THE FIRST DERIVATIVE
; C             LOCALLY.
; C
         DA = MAX([MIN([RELPER*ABS(A),ABSDX]), 100.0D0*SMALL*ABS(A)])
         DA = (DX GE 0)?(+DA):(-DA)
         IF (DA EQ 0.0D0) THEN DA = RELPER*DX

         SF = CALL_FUNCTION(DFNAME, DF, A+DA, Y, PRIVATE, _EXTRA=FA)
         if ddeabm_funcerror NE 0 then return

         YP = SF - YPRIME
         DELF = max(abs(YP))
         DFDXB = BIG
         IF (DELF LT BIG*ABS(DA)) THEN DFDXB = DELF/ABS(DA)
         FBND = max(abs(SF))
; C
; C        ...............................................................
; C
; C             COMPUTE AN ESTIMATE (DFDUB) OF THE LOCAL LIPSCHITZ
; C             CONSTANT FOR THE SYSTEM OF DIFFERENTIAL EQUATIONS. THIS
; C             ALSO REPRESENTS AN ESTIMATE OF THE NORM OF THE JACOBIAN
; C             LOCALLY.  THREE ITERATIONS (TWO WHEN NEQ=1) ARE USED TO
; C             ESTIMATE THE LIPSCHITZ CONSTANT BY NUMERICAL DIFFERENCES.
; C             THE FIRST PERTURBATION VECTOR IS BASED ON THE INITIAL
; C             DERIVATIVES AND DIRECTION OF INTEGRATION. THE SECOND
; C             PERTURBATION VECTOR IS FORMED USING ANOTHER EVALUATION OF
; C             THE DIFFERENTIAL EQUATION.  THE THIRD PERTURBATION VECTOR
; C             IS FORMED USING PERTURBATIONS BASED ONLY ON THE INITIAL
; C             VALUES. COMPONENTS THAT ARE ZERO ARE ALWAYS CHANGED TO
; C             NON-ZERO VALUES (EXCEPT ON THE FIRST ITERATION). WHEN
; C             INFORMATION IS AVAILABLE, CARE IS TAKEN TO ENSURE THAT
; C             COMPONENTS OF THE PERTURBATION VECTOR HAVE SIGNS WHICH ARE
; C             CONSISTENT WITH THE SLOPES OF LOCAL SOLUTION CURVES.
; C             ALSO CHOOSE THE LARGEST BOUND (FBND) FOR THE FIRST
; C             DERIVATIVE.
; C
; C                               PERTURBATION VECTOR SIZE IS HELD
; C                               CONSTANT FOR ALL ITERATIONS. COMPUTE
; C                               THIS CHANGE FROM THE
; C                                       SIZE OF THE VECTOR OF INITIAL
; C                                       VALUES.
         DELY = RELPER*max(abs(y))
         IF (DELY EQ 0.0D0) THEN DELY = RELPER
         DELY = (DX GE 0)?(+DELY):(-DELY)
         DELF = max(abs(YPRIME))
         FBND = MAX([FBND,DELF])
         IF (DELF NE 0.0D0) THEN BEGIN
; C           USE INITIAL DERIVATIVES FOR FIRST PERTURBATION
             SPY = YPRIME
             YP = YPRIME
         ENDIF ELSE BEGIN
; C           CANNOT HAVE A NULL PERTURBATION VECTOR
             SPY(*) = 0
             YP(*) = 1
             DELF = max(abs(yp))
         ENDELSE
; C
         DFDUB = 0.0D0
         LK = MIN([NEQ+1,3])
         FOR K = 1L, LK DO BEGIN

; C           DEFINE PERTURBED VECTOR OF INITIAL VALUES
             PV = Y + DELY*(YP/DELF)
             IF (K NE 2) THEN BEGIN
; C              EVALUATE DERIVATIVES ASSOCIATED WITH PERTURBED
; C              VECTOR  AND  COMPUTE CORRESPONDING DIFFERENCES

                 YP = CALL_FUNCTION(DFNAME, DF, A, PV, PRIVATE, _EXTRA=FA)
                 if ddeabm_funcerror NE 0 then return
                 PV = YP - YPRIME
            ENDIF ELSE BEGIN
; C              USE A SHIFTED VALUE OF THE INDEPENDENT VARIABLE
; C                                    IN COMPUTING ONE ESTIMATE

                YP = CALL_FUNCTION(DFNAME, DF, A+DA, PV, PRIVATE, _EXTRA=FA)
                if ddeabm_funcerror NE 0 then return
                PV = YP - SF
            ENDELSE
; C           CHOOSE LARGEST BOUNDS ON THE FIRST DERIVATIVE
; C                          AND A LOCAL LIPSCHITZ CONSTANT
            FBND = MAX([FBND,max(abs(yp))])
            DELF = max(abs(pv))
; C        ...EXIT
            IF (DELF GE BIG*ABS(DELY)) THEN GOTO, DHSTRT_150
            DFDUB = MAX([DFDUB,DELF/ABS(DELY)])
; C     ......EXIT
            IF (K EQ LK) THEN GOTO, DHSTRT_160
; C           CHOOSE NEXT PERTURBATION VECTOR
            IF (DELF EQ 0.0D0) THEN DELF = 1.0D0
            IF (K NE 2) THEN BEGIN
                DY = ABS(PV)
                wh = where(dy EQ 0, ct)
                if ct GT 0 then dy(wh) = DELF
            endif else begin
                DY = Y
                wh = where(dy EQ 0, ct)
                if ct GT 0 then dy(wh) = DELY/RELPER
            endelse
            wh = where(spy EQ 0, ct)
            if ct GT 0 then spy(wh) = yp(wh)
            wh = where(spy LT 0, ct)
            if ct GT 0 then dy(wh) = -dy(wh)
            yp(*) = dy
            DELF = max(abs(YP))
        ENDFOR
        DHSTRT_150:
; C
; C        PROTECT AGAINST AN OVERFLOW
        DFDUB = BIG
        DHSTRT_160:
; C
; C     ..................................................................
; C
; C          COMPUTE A BOUND (YDPB) ON THE NORM OF THE SECOND DERIVATIVE
; C
      YDPB = DFDXB + DFDUB*FBND
; C
; C     ..................................................................
; C
; C          DEFINE THE TOLERANCE PARAMETER UPON WHICH THE STARTING STEP
; C          SIZE IS TO BE BASED.  A VALUE IN THE MIDDLE OF THE ERROR
; C          TOLERANCE RANGE IS SELECTED.
; C
      TOLEXP = ALOG10(ETOL)
      TOLMIN = MIN(TOLEXP)
      TOLSUM = TOTAL(TOLEXP)
      TOLP = 10.0D0^(0.5D0*(TOLSUM/NEQ + TOLMIN)/(MORDER+1))
; C
; C     ..................................................................
; C
; C          COMPUTE A STARTING STEP SIZE BASED ON THE ABOVE FIRST AND
; C          SECOND DERIVATIVE INFORMATION
; C
; C                            RESTRICT THE STEP LENGTH TO BE NOT BIGGER
; C                            THAN ABS(B-A).   (UNLESS  B  IS TOO CLOSE
; C                            TO  A)
      H = ABSDX
; C
      IF (YDPB EQ 0.0D0 AND FBND EQ 0.0D0) THEN BEGIN
; GO TO 180
; C
; C        BOTH FIRST DERIVATIVE TERM (FBND) AND SECOND
; C                     DERIVATIVE TERM (YDPB) ARE ZERO
         IF (TOLP LT 1.0D0) THEN H = ABSDX*TOLP
     ENDIF ELSE IF (YDPB EQ 0.0D0) THEN BEGIN
; C        ONLY SECOND DERIVATIVE TERM (YDPB) IS ZERO
         IF (TOLP LT FBND*ABSDX) THEN H = TOLP/FBND
     ENDIF ELSE BEGIN
; C        SECOND DERIVATIVE TERM (YDPB) IS NON-ZERO
         SRYDPB = SQRT(0.5D0*YDPB)
         IF (TOLP LT SRYDPB*ABSDX) THEN H = TOLP/SRYDPB
     ENDELSE

; C     FURTHER RESTRICT THE STEP LENGTH TO BE NOT
; C                               BIGGER THAN  1/DFDUB
      IF (H*DFDUB GT 1.0D0) THEN H = 1.0D0/DFDUB

; C     FINALLY, RESTRICT THE STEP LENGTH TO BE NOT
; C     SMALLER THAN  100*SMALL*ABS(A).  HOWEVER, IF
; C     A=0. AND THE COMPUTED H UNDERFLOWED TO ZERO,
; C     THE ALGORITHM RETURNS  SMALL*ABS(B)  FOR THE
; C                                     STEP LENGTH.
      H = MAX([H,100.0D0*SMALL*ABS(A)])
      IF (H EQ 0.0D0) THEN H = SMALL*ABS(B)

; C     NOW SET DIRECTION OF INTEGRATION
      H = (DX GE 0)?(+H):(-H)

      RETURN
  END





; *DECK DDEABM_DDES

      pro ddeabm_ddes, DF, NEQ, T, Y, TOUT, INFO, RTOL, ATOL, IDID, $
        YPOUT, YP, YY, WT, P, PHI, ALPHA, BETA, PSI, V, W, SIG, G, GI, $
        H, EPS, X, XOLD, HOLD, TOLD, DELSGN, TSTOP, TWOU, FOURU, START, $
        PHASE1, NORND, STIFF, INTOUT, NS, KORD, KOLD, INIT, KSTEPS, $
        KLE4, IQUIT, KPREV, IVC, IV, KGI, PRIVATE, FA, dfname, $
                errmsg=errmsg, max_stepsize=max_stepsize
; C***BEGIN PROLOGUE  DDEABM_DDES
; C***SUBSIDIARY
; C***PURPOSE  Subsidiary to DDEABM
; C***LIBRARY   SLATEC
; C***TYPE      DOUBLE PRECISION (DES-S, DDES-D)
; C***AUTHOR  Watts, H. A., (SNLA)
; C***DESCRIPTION
; C
; C   DDEABM merely allocates storage for DDEABM_DDES to relieve the user of the
; C   inconvenience of a long call list.  Consequently  DDEABM_DDES  is used as
; C   described in the comments for  DDEABM .
; C
; C***SEE ALSO  DDEABM
; C***ROUTINES CALLED  D1MACH, DDEABM_DINTP, DDEABM_DSTEPS, XERMSG
; C***REVISION HISTORY  (YYMMDD)
; C   820301  DATE WRITTEN
; C   890531  Changed all specific intrinsics to generic.  (WRB)
; C   890831  Modified array declarations.  (WRB)
; C   891214  Prologue converted to Version 4.0 format.  (BAB)
; C   900328  Added TYPE section.  (WRB)
; C   900510  Convert XERRWV calls to XERMSG calls, cvt GOTOs to
; C           IF-THEN-ELSE.  (RWC)
; C   910722  Updated AUTHOR section.  (ALS)
; C***END PROLOGUE  DDEABM_DDES
; C
;       INTEGER IDID, INFO, INIT, IPAR, IQUIT, IV, IVC, K, KGI, KLE4,
;      1      KOLD, KORD, KPREV, KSTEPS, L, LTOL, MAXNUM, NATOLP, NEQ,
;      2      NRTOLP, NS
;       DOUBLE PRECISION A, ABSDEL, ALPHA, ATOL, BETA, D1MACH,
;      1      DEL, DELSGN, DT, EPS, FOURU, G, GI, H,
;      2      HA, HOLD, P, PHI, PSI, RPAR, RTOL, SIG, T, TOLD, TOUT,
;      3      TSTOP, TWOU, U, V, W, WT, X, XOLD, Y, YP, YPOUT, YY
;       LOGICAL STIFF,CRASH,START,PHASE1,NORND,INTOUT
; C
;      DIMENSION Y(*),YY(*),WT(*),PHI(NEQ,16),P(*),YP(*),
;     1  YPOUT(*),PSI(12),ALPHA(12),BETA(12),SIG(13),V(12),W(12),G(13),
;     2  GI(11),IV(10),INFO(15),RTOL(*),ATOL(*),RPAR(*),IPAR(*)
;      CHARACTER*8 XERN1
;      CHARACTER*16 XERN3, XERN4
; C
;       EXTERNAL DF
; C
; C.......................................................................
; C
; C  THE EXPENSE OF SOLVING THE PROBLEM IS MONITORED BY COUNTING THE
; C  NUMBER OF  STEPS ATTEMPTED. WHEN THIS EXCEEDS  MAXNUM, THE COUNTER
; C  IS RESET TO ZERO AND THE USER IS INFORMED ABOUT POSSIBLE EXCESSIVE
; C  WORK.
; C
      MAXNUM = 500L
; C
; C.......................................................................
; C
; C***FIRST EXECUTABLE STATEMENT  DDEABM_DDES
      common ddeabm_func_common
      IF (INFO(1-1) EQ 0) THEN BEGIN
; C
; C ON THE FIRST CALL , PERFORM INITIALIZATION --
; C        DEFINE THE MACHINE UNIT ROUNDOFF QUANTITY  U  BY CALLING THE
; C        FUNCTION ROUTINE  D1MACH. THE USER MUST MAKE SURE THAT THE
; C        VALUES SET IN D1MACH ARE RELEVANT TO THE COMPUTER BEING USED.
; C
         U=(machar(/double)).eps  ;; XXX
; C                       -- SET ASSOCIATED MACHINE DEPENDENT PARAMETERS
         TWOU=2.D0*U
         FOURU=4.D0*U
; C                       -- SET TERMINATION FLAG
         IQUIT=0L
; C                       -- SET INITIALIZATION INDICATOR
         INIT=0L
; C                       -- SET COUNTER FOR ATTEMPTED STEPS
         KSTEPS=0L
; C                       -- SET INDICATOR FOR INTERMEDIATE-OUTPUT
         INTOUT= 0L
; C                       -- SET INDICATOR FOR STIFFNESS DETECTION
         STIFF= 0L
; C                       -- SET STEP COUNTER FOR STIFFNESS DETECTION
         KLE4=0L
; C                       -- SET INDICATORS FOR STEPS CODE
         START= 1L
         PHASE1= 1L
         NORND= 1L
; C                       -- RESET INFO(1) FOR SUBSEQUENT CALLS
         INFO(1-1)=1L
      ENDIF
; C
; C.......................................................................
; C
; C      CHECK VALIDITY OF INPUT PARAMETERS ON EACH ENTRY
; C
      IF (INFO(1-1) NE 0  AND  INFO(1-1) NE 1) THEN BEGIN
          errmsg = 'IN DDEABM, INFO(1-1) MUST BE '+ $
            'SET TO 0 FOR THE START OF A NEW PROBLEM, AND MUST BE '+ $
            'SET TO 1 FOLLOWING AN INTERRUPTED TASK.  YOU ARE '+ $
            'ATTEMPTING TO CONTINUE THE INTEGRATION ILLEGALLY BY '+ $
            'CALLING THE CODE WITH INFO(1-1) = ' + strtrim(info(1-1),2)
          IDID=-33L
      ENDIF

      IF (INFO(2-1) NE 0  AND  INFO(2-1) NE 1) THEN BEGIN
          errmsg = 'IN DDEABM, INFO(2-1) MUST BE '+ $
            '0 OR 1 INDICATING SCALAR AND VECTOR ERROR TOLERANCES, '+ $
            'RESPECTIVELY.  YOU HAVE CALLED THE CODE WITH INFO(2-1) = '+ $
            strtrim(info(2-1),2)
         IDID=-33L
      ENDIF

      IF (INFO(3-1) NE 0  AND  INFO(3-1) NE 1) THEN BEGIN
          errmsg = 'IN DDEABM, INFO(3-1) MUST BE '+ $
            '0 OR 1 INDICATING THE INTERVAL OR INTERMEDIATE-OUTPUT '+ $
            'MODE OF INTEGRATION, RESPECTIVELY.  YOU HAVE CALLED '+ $
            'THE CODE WITH  INFO(3-1) = '+strtrim(info(3-1),2)
         IDID=-33L
      ENDIF

      IF (INFO(4-1) NE 0  AND  INFO(4-1) NE 1) THEN BEGIN
         errmsg = 'IN DDEABM, INFO(4-1) MUST BE '+ $
           '0 OR 1 INDICATING WHETHER OR NOT THE INTEGRATION '+ $
           'INTERVAL IS TO BE RESTRICTED BY A POINT TSTOP.  YOU '+ $
           'HAVE CALLED THE CODE WITH INFO(4-1) = '+strtrim(info(4-1),2)
         IDID=-33L
      ENDIF

      IF (NEQ LT 1) THEN BEGIN
         errmsg = 'IN DDEABM,  THE NUMBER OF '+ $
           'EQUATIONS NEQ MUST BE A POSITIVE INTEGER.  YOU HAVE '+ $
           'CALLED THE CODE WITH  NEQ = '+strtrim(neq,2)
         IDID=-33L
      ENDIF

      whr = where(rtol LT 0, nrtolp)
      wha = where(atol LT 0, natolp)

      NRTOLP = total( (RTOL LT 0) ) NE 0
      NATOLP = total( (ATOL LT 0) ) NE 0

      if nrtolp NE 0 then begin
          errmsg = 'IN DDEABM, THE RELATIVE '+ $
            'ERROR TOLERANCES RTOL MUST BE NON-NEGATIVE.  YOU '+ $
            'HAVE CALLED THE CODE WITH  RTOL('+strtrim(whr(0),2)+') = '+ $
            strtrim(rtol(whr(0)),2)+ $
            '. IN THE CASE OF VECTOR ERROR TOLERANCES, '+$
            'NO FURTHER CHECKING OF RTOL COMPONENTS IS DONE.'
          IDID = -33L
      endif else if natolp NE 0 then begin
          errmsg = 'IN DDEABM, THE ABSOLUTE '+ $
            'ERROR TOLERANCES ATOL MUST BE NON-NEGATIVE.  YOU '+ $
            'HAVE CALLED THE CODE WITH  ATOL('+strtrim(wha(0),2)+'-1) = '+ $
            strtrim(atol(wha(0)),2)+ $
            '.  IN THE CASE OF VECTOR ERROR TOLERANCES, '+ $
            'NO FURTHER CHECKING OF ATOL COMPONENTS IS DONE.'
          IDID = -33L
      ENDIF

     DDES_100:
     IF (INFO(4-1) EQ 1) THEN BEGIN
         IF ( (TOUT-T)*(TSTOP-T) LT 0 $
              OR ABS(TOUT-T) GT ABS(TSTOP-T)) THEN BEGIN
            errmsg = 'IN DDEABM, YOU HAVE '+ $
              'CALLED THE CODE WITH  TOUT = '+strtrim(tout,2)+' BUT '+ $
              'YOU HAVE ALSO TOLD THE CODE (INFO(4-1) = 1) NOT TO '+ $
              'INTEGRATE PAST THE POINT TSTOP = '+strtrim(tstop,2)+ $
              ' THESE INSTRUCTIONS CONFLICT.'
            IDID=-33L
         ENDIF
      ENDIF
; C
; C     CHECK SOME CONTINUATION POSSIBILITIES
; C
      IF (INIT NE 0) THEN BEGIN
         IF (T EQ TOUT) THEN BEGIN
            errmsg = 'IN DDEABM, YOU HAVE '+ $
              'CALLED THE CODE WITH  T = TOUT = '+strtrim(T,2)+ $
              '.  THIS IS NOT ALLOWED ON CONTINUATION CALLS.'
            IDID=-33L
         ENDIF

         IF (T NE TOLD) THEN BEGIN
            errmsg = 'IN DDEABM, YOU HAVE '+ $
              'CHANGED THE VALUE OF T FROM '+strtrim(told,2)+' TO '+ $
              strtrim(t,2)+'  THIS IS NOT ALLOWED ON CONTINUATION CALLS.'
            IDID=-33L
         ENDIF

         IF (INIT NE 1) THEN BEGIN
            IF (DELSGN*(TOUT-T) LT 0.D0) THEN BEGIN
               errmsg = 'IN DDEABM, BY '+ $
                 'CALLING THE CODE WITH TOUT = '+strtrim(tout,2)+ $
                 ' YOU ARE ATTEMPTING TO CHANGE THE DIRECTION OF '+ $
                 'INTEGRATION.  THIS IS NOT ALLOWED WITHOUT '+ $
                 'RESTARTING.'
               IDID=-33L
            ENDIF
         ENDIF
      ENDIF
; C
; C     INVALID INPUT DETECTED
; C
      IF (IDID EQ (-33)) THEN BEGIN
         IF (IQUIT NE (-33)) THEN BEGIN
            IQUIT = -33L
            INFO(1-1) = -1L
         ENDIF ELSE BEGIN
             errmsg = 'IN DDEABM, INVALID '+ $
               'INPUT WAS DETECTED ON SUCCESSIVE ENTRIES.  IT IS '+ $
               'IMPOSSIBLE TO PROCEED BECAUSE YOU HAVE NOT '+ $
               'CORRECTED THE PROBLEM, SO EXECUTION IS BEING '+ $
               'TERMINATED.'
         ENDELSE
         RETURN
      ENDIF
; C
; C.......................................................................
; C
; C     RTOL = ATOL = 0. IS ALLOWED AS VALID INPUT AND INTERPRETED AS
; C     ASKING FOR THE MOST ACCURATE SOLUTION POSSIBLE. IN THIS CASE,
; C     THE RELATIVE ERROR TOLERANCE RTOL IS RESET TO THE SMALLEST VALUE
; C     FOURU WHICH IS LIKELY TO BE REASONABLE FOR THIS METHOD AND MACHINE
; C
      wh = where(rtol+atol EQ 0, ct)
      if ct GT 0 then begin
          ;; Expand RTOL if necessary to be per-vector
          if n_elements(rtol) LT n_elements(atol) then $
            rtol = rtol + atol*0
          rtol(wh) = fouru
          IDID = -2L
      endif

      DDES_190:
      IF (IDID EQ (-2)) THEN BEGIN
; C                       RTOL=ATOL=0 ON INPUT, SO RTOL IS CHANGED TO A
; C                                                SMALL POSITIVE VALUE
          INFO(1-1)=-1
          RETURN
      ENDIF
; C
; C     BRANCH ON STATUS OF INITIALIZATION INDICATOR
; C            INIT=0 MEANS INITIAL DERIVATIVES AND NOMINAL STEP SIZE
; C                   AND DIRECTION NOT YET SET
; C            INIT=1 MEANS NOMINAL STEP SIZE AND DIRECTION NOT YET SET
; C            INIT=2 MEANS NO FURTHER INITIALIZATION REQUIRED
; C
      IF (INIT EQ 0) THEN GOTO, DDES_210
      IF (INIT EQ 1) THEN GOTO, DDES_220
      GOTO, DDES_240
; C
; C.......................................................................
; C
; C     MORE INITIALIZATION --
; C                         -- EVALUATE INITIAL DERIVATIVES
; C
      DDES_210:
      INIT=1L
      A=T
      YP = CALL_FUNCTION(DFNAME, DF, A, Y, PRIVATE, _EXTRA=FA)
      if ddeabm_funcerror NE 0 then return
      

      IF (T EQ TOUT) THEN BEGIN
          IDID=2L
          YPOUT = YP
          TOLD=T
          RETURN
      ENDIF
; C
; C                         -- SET INDEPENDENT AND DEPENDENT VARIABLES
; C                                              X AND YY(*) FOR STEPS
; C                         -- SET SIGN OF INTEGRATION DIRECTION
; C                         -- INITIALIZE THE STEP SIZE
; C
      DDES_220:
      INIT = 2L
      X = T
      YY = Y
      DELSGN = (TOUT GE T)?(+1):(-1)
      DELSGNX = (TOUT GE X)?(+1):(-1)
      H = MAX([FOURU*ABS(X),ABS(TOUT-X)])
      if n_elements(max_stepsize) GT 0 then h = h < max_stepsize(0)
      H = H * DELSGNX
; C
; C.......................................................................
; C
; C   ON EACH CALL SET INFORMATION WHICH DETERMINES THE ALLOWED INTERVAL
; C   OF INTEGRATION BEFORE RETURNING WITH AN ANSWER AT TOUT
; C
      DDES_240:
      DEL = TOUT - T
      ABSDEL = ABS(DEL)
; C
; C.......................................................................
; C
; C   IF ALREADY PAST OUTPUT POINT, INTERPOLATE AND RETURN
; C
      DDES_250:
      IF (ABS(X-T) GE ABSDEL) THEN BEGIN
          DDEABM_DINTP, X,YY,TOUT,Y,YPOUT,NEQ,KOLD,PHI,IVC,IV,KGI,GI, $
                                             ALPHA,G,W,XOLD,P
          IDID = 3L
          IF (X EQ TOUT) THEN BEGIN
              IDID = 2L
              INTOUT = 0L
          ENDIF
          T = TOUT
          TOLD = T
          RETURN
      ENDIF
; C
; C   IF CANNOT GO PAST TSTOP AND SUFFICIENTLY CLOSE,
; C   EXTRAPOLATE AND RETURN
; C
      IF ((INFO(4-1) EQ 1) AND $
          (ABS(TSTOP-X) LT FOURU*ABS(X))) THEN BEGIN
          DT = TOUT - X
          Y = YY + DT*YP

          YPOUT = CALL_FUNCTION(DFNAME, DF, TOUT, Y, PRIVATE, _EXTRA=FA)
          if ddeabm_funcerror NE 0 then return
          IDID = 3L
          T = TOUT
          TOLD = T
          RETURN
      ENDIF

      IF (INFO(3-1) EQ 0  OR  INTOUT EQ 0) EQ 0 THEN BEGIN
; C
; C   INTERMEDIATE-OUTPUT MODE
; C
          IDID = 1L
          Y = YY
          YPOUT = YP
          T = X
          TOLD = T
          INTOUT = 0L
          RETURN
      ENDIF
; C
; C.......................................................................
; C
; C     MONITOR NUMBER OF STEPS ATTEMPTED
; C
      IF (KSTEPS GT MAXNUM) THEN BEGIN
; C
; C                       A SIGNIFICANT AMOUNT OF WORK HAS BEEN EXPENDED
          IDID=-1L
          KSTEPS=0L
          IF (STIFF) THEN BEGIN
; C
; C                       PROBLEM APPEARS TO BE STIFF
              IDID=-4L
              STIFF= 0L
              KLE4=0L
          ENDIF

          Y = YY
          YPOUT = YP
          T = X
          TOLD = T
          INFO(1-1) = -1
          INTOUT = 0L
          RETURN
      ENDIF
; C
; C.......................................................................
; C
; C   LIMIT STEP SIZE, SET WEIGHT VECTOR AND TAKE A STEP
; C
      HA = ABS(H)
      IF (INFO(4-1) EQ 1) THEN BEGIN
          HA = MIN([HA,ABS(TSTOP-X)])
      ENDIF
      H = (H GE 0)?(HA):(-HA)
      EPS = 1.0D0
      LTOL = 1L

      wt = rtol*abs(yy) + atol
      wh = where(wt LE 0, ct)

; C
; C                       RELATIVE ERROR CRITERION INAPPROPRIATE
      if ct GT 0 then begin
          IDID = -3L
          Y = YY
          YPOUT = YP
          T = X
          TOLD = T
          INFO(1-1) = -1L
          INTOUT = 0L
          RETURN
      endif

      DDES_380:
      DDEABM_DSTEPS, DF,NEQ,YY,X,H,EPS,WT,START,HOLD,KORD,KOLD,CRASH,PHI,P, $
                YP,PSI,ALPHA,BETA,SIG,V,W,G,PHASE1,NS,NORND,KSTEPS, $
                TWOU,FOURU,XOLD,KPREV,IVC,IV,KGI,GI,PRIVATE,FA, DFNAME, $
                max_stepsize=max_stepsize
      if ddeabm_funcerror NE 0 then return

; C
; C.......................................................................
; C
      IF (CRASH) THEN BEGIN
; C
; C                       TOLERANCES TOO SMALL
          IDID = -2L
          RTOL = EPS*RTOL
          ATOL = EPS*ATOL
          Y = YY
          YPOUT = YP
          T = X
          TOLD = T
          INFO(1-1) = -1L
          INTOUT = 0L
          RETURN
      ENDIF
; C
; C   (STIFFNESS TEST) COUNT NUMBER OF CONSECUTIVE STEPS TAKEN WITH THE
; C   ORDER OF THE METHOD BEING LESS OR EQUAL TO FOUR
; C
  DDES_420:
      KLE4 = KLE4 + 1
      IF(KOLD GT 4) THEN KLE4 = 0L
      IF(KLE4 GE 50) THEN STIFF = 1L
      INTOUT = 1L
      GOTO, DDES_250
    END


; *DECK DDEABM_DINTP
pro DDEABM_DINTP, X, Y, XOUT, YOUT, YPOUT, NEQN, KOLD, PHI, IVC, $
       IV, KGI, GI, ALPHA, OG, OW, OX, OY
; C***BEGIN PROLOGUE  DDEABM_DINTP
; C***PURPOSE  Approximate the solution at XOUT by evaluating the
; C            polynomial computed in DDEABM_DSTEPS at XOUT.  Must be used in
; C            conjunction with DDEABM_DSTEPS.
; C***LIBRARY   SLATEC (DEPAC)
; C***CATEGORY  I1A1B
; C***TYPE      DOUBLE PRECISION (SINTRP-S, DINTP-D)
; C***KEYWORDS  ADAMS METHOD, DEPAC, INITIAL VALUE PROBLEMS, ODE,
; C             ORDINARY DIFFERENTIAL EQUATIONS, PREDICTOR-CORRECTOR,
; C             SMOOTH INTERPOLANT
; C***AUTHOR  Watts, H. A., (SNLA)
; C***DESCRIPTION
; C
; C   The methods in subroutine  DDEABM_DSTEPS  approximate the solution near  X
; C   by a polynomial.  Subroutine  DDEABM_DINTP  approximates the solution at
; C   XOUT  by evaluating the polynomial there.  Information defining this
; C   polynomial is passed from  DDEABM_DSTEPS  so  DDEABM_DINTP  cannot be used alone.
; C
; C   Subroutine DDEABM_DSTEPS is completely explained and documented in the text
; C   "Computer Solution of Ordinary Differential Equations, the Initial
; C   Value Problem"  by L. F. Shampine and M. K. Gordon.
; C
; C   Input to DDEABM_DINTP --
; C
; C   The user provides storage in the calling program for the arrays in
; C   the call list
; C      DIMENSION Y(NEQN),YOUT(NEQN),YPOUT(NEQN),PHI(NEQN,16),OY(NEQN)
; C                AND ALPHA(12),OG(13),OW(12),GI(11),IV(10)
; C   and defines
; C      XOUT -- point at which solution is desired.
; C   The remaining parameters are defined in  DDEABM_DSTEPS  and passed to
; C   DDEABM_DINTP  from that subroutine
; C
; C   Output from  DDEABM_DINTP --
; C
; C      YOUT(*) -- solution at  XOUT
; C      YPOUT(*) -- derivative of solution at  XOUT
; C   The remaining parameters are returned unaltered from their input
; C   values.  Integration with  DDEABM_DSTEPS  may be continued.
; C
; C***REFERENCES  H. A. Watts, A smoother interpolant for DE/STEP, INTRP
; C                 II, Report SAND84-0293, Sandia Laboratories, 1984.
; C***ROUTINES CALLED  (NONE)
; C***REVISION HISTORY  (YYMMDD)
; C   840201  DATE WRITTEN
; C   890831  Modified array declarations.  (WRB)
; C   890831  REVISION DATE from Version 3.2
; C   891214  Prologue converted to Version 4.0 format.  (BAB)
; C   920501  Reformatted the REFERENCES section.  (WRB)
; C***END PROLOGUE  DDEABM_DINTP
; C
;       INTEGER I, IQ, IV, IVC, IW, J, JQ, KGI, KOLD, KP1, KP2,
;      1        L, M, NEQN
;       DOUBLE PRECISION ALP, ALPHA, C, G, GDI, GDIF, GI, GAMMA, H, HI,
;      1       HMU, OG, OW, OX, OY, PHI, RMU, SIGMA, TEMP1, TEMP2, TEMP3,
;      2       W, X, XI, XIM1, XIQ, XOUT, Y, YOUT, YPOUT
; C
;       DIMENSION Y(*),YOUT(*),YPOUT(*),PHI(NEQN,16),OY(*)
;       DIMENSION G(13),C(13),W(13),OG(13),OW(12),ALPHA(12),GI(11),IV(10)
; C
; C***FIRST EXECUTABLE STATEMENT  DDEABM_DINTP
      KP1 = KOLD + 1
      KP2 = KOLD + 2

      HI = XOUT - OX
      H = X - OX
      XI = HI/H
      XIM1 = XI - 1.D0

      G = DBLARR(13) & C = G & W = G
; C
; C   INITIALIZE W(*) FOR COMPUTING G(*)
; C
      XIQ = XI
      IQ = dindgen(kp1)+1
      XIQ = XI^(IQ+1)
      W(0:KP1-1) = XIQ/(IQ*(IQ+1))

; C
; C   COMPUTE THE DOUBLE INTEGRAL TERM GDI
; C
      IF (KOLD LE KGI) THEN BEGIN
          GDI = GI(KOLD-1)
          GOTO, DINTP_60
      ENDIF
          
      IF (IVC LE 0) THEN BEGIN
          GDI = 1.0D0/(KP1*(KP1+1))
          M = 2L
      ENDIF ELSE BEGIN
          IW = IV(IVC-1)
          GDI = OW(IW-1)
          M = KOLD - IW + 3
      ENDELSE

      IF (M LE KOLD) THEN BEGIN
;;        XXX:    (M>1) is a kludge
          FOR I = (M>1), KOLD DO $
            GDI = OW(KP2-I-1) - ALPHA(I-1)*GDI
      ENDIF

; C
; C   COMPUTE G(*) AND C(*)
; C
      DINTP_60:
      G(0:1) = [XI, XI^2/2]
      C(0:1) = [1d, XI]
      IF (KOLD GE 2) THEN BEGIN
          FOR I = 2L, KOLD DO BEGIN
              ALP = ALPHA(I-1)
              GAMMA = 1.0D0 + XIM1*ALP
              L = KP2 - I
              W(0:L-1) = GAMMA*W(0:L-1) - ALP*W(1:L)
              G(I) = W(0)
              C(I) = GAMMA*C(I-1)
          ENDFOR
      ENDIF
; C
; C   DEFINE INTERPOLATION PARAMETERS
; C
      SIGMA = (W(1) - XIM1*W(0))/GDI
      RMU = XIM1*C(KOLD)/GDI  ;; *** NOTE: KP1-1 is KOLD
      HMU = RMU/H

; C   INTERPOLATE FOR THE SOLUTION -- YOUT
; C   AND FOR THE DERIVATIVE OF THE SOLUTION -- YPOUT

      YOUT(*) = 0
      YPOUT(*) = 0
      J = lindgen(KOLD)+1
      I = KP2 - J - 1 ;; *** NOTE: -1 is here
      GDIF = OG(I) - OG(I-1)
      TEMP2 = (G(I) - G(I-1)) - SIGMA*GDIF
      TEMP3 = (C(I) - C(I-1)) + RMU*GDIF

      FOR J = 0L, KOLD-1 DO BEGIN
        YOUT  = YOUT  + TEMP2(j)*PHI(*,I(j))
        YPOUT = YPOUT + TEMP3(j)*PHI(*,I(j))
      ENDFOR

      YOUT = ( ((1.0D0 - SIGMA)*OY + SIGMA*Y) + $
               H*(YOUT + (G(0) - SIGMA*OG(0))*PHI(*,0)) )
      YPOUT = ( HMU*(OY - Y) + $
                (YPOUT + (C(0) + RMU*OG(0))*PHI(*,0)) )

      RETURN
  END

; *DECK DDEABM_DSTEPS
  pro DDEABM_DSTEPS, DF, NEQN, Y, X, H, EPS, WT, START, HOLD, K, $
              KOLD, CRASH, PHI, P, YP, PSI, ALPHA, BETA, SIG, V, W, G, $
              PHASE1, NS, NORND, KSTEPS, TWOU, FOURU, XOLD, KPREV, IVC, IV, $
              KGI, GI, PRIVATE, FA, dfname, max_stepsize=max_stepsize
; C***BEGIN PROLOGUE  DDEABM_DSTEPS
; C***PURPOSE  Integrate a system of first order ordinary differential
; C            equations one step.
; C***LIBRARY   SLATEC (DEPAC)
; C***CATEGORY  I1A1B
; C***TYPE      DOUBLE PRECISION (STEPS-S, DSTEPS-D)
; C***KEYWORDS  ADAMS METHOD, DEPAC, INITIAL VALUE PROBLEMS, ODE,
; C             ORDINARY DIFFERENTIAL EQUATIONS, PREDICTOR-CORRECTOR
; C***AUTHOR  Shampine, L. F., (SNLA)
; C           Gordon, M. K., (SNLA)
; C             MODIFIED BY H.A. WATTS
; C***DESCRIPTION
; C
; C   Written by L. F. Shampine and M. K. Gordon
; C
; C   Abstract
; C
; C   Subroutine  DDEABM_DSTEPS  is normally used indirectly through subroutine
; C   DDEABM .  Because  DDEABM  suffices for most problems and is much
; C   easier to use, using it should be considered before using  DDEABM_DSTEPS
; C   alone.
; C
; C   Subroutine DDEABM_DSTEPS integrates a system of  NEQN  first order ordinary
; C   differential equations one step, normally from X to X+H, using a
; C   modified divided difference form of the Adams Pece formulas.  Local
; C   extrapolation is used to improve absolute stability and accuracy.
; C   The code adjusts its order and step size to control the local error
; C   per unit step in a generalized sense.  Special devices are included
; C   to control roundoff error and to detect when the user is requesting
; C   too much accuracy.
; C
; C   This code is completely explained and documented in the text,
; C   Computer Solution of Ordinary Differential Equations, The Initial
; C   Value Problem  by L. F. Shampine and M. K. Gordon.
; C   Further details on use of this code are available in "Solving
; C   Ordinary Differential Equations with ODE, STEP, and INTRP",
; C   by L. F. Shampine and M. K. Gordon, SLA-73-1060.
; C
; C
; C   The parameters represent --
; C      DF -- subroutine to evaluate derivatives
; C      NEQN -- number of equations to be integrated
; C      Y(*) -- solution vector at X
; C      X -- independent variable
; C      H -- appropriate step size for next step.  Normally determined by
; C           code
; C      EPS -- local error tolerance
; C      WT(*) -- vector of weights for error criterion
; C      START -- logical variable set .TRUE. for first step,  .FALSE.
; C           otherwise
; C      HOLD -- step size used for last successful step
; C      K -- appropriate order for next step (determined by code)
; C      KOLD -- order used for last successful step
; C      CRASH -- logical variable set .TRUE. when no step can be taken,
; C           .FALSE. otherwise.
; C      YP(*) -- derivative of solution vector at  X  after successful
; C           step
; C      KSTEPS -- counter on attempted steps
; C      TWOU -- 2.*U where U is machine unit roundoff quantity
; C      FOURU -- 4.*U where U is machine unit roundoff quantity
; C      RPAR,IPAR -- parameter arrays which you may choose to use
; C            for communication between your program and subroutine F.
; C            They are not altered or used by DDEABM_DSTEPS.
; C   The variables X,XOLD,KOLD,KGI and IVC and the arrays Y,PHI,ALPHA,G,
; C   W,P,IV and GI are required for the interpolation subroutine SINTRP.
; C   The remaining variables and arrays are included in the call list
; C   only to eliminate local retention of variables between calls.
; C
; C   Input to DDEABM_DSTEPS
; C
; C      First call --
; C
; C   The user must provide storage in his calling program for all arrays
; C   in the call list, namely
; C
; C     DIMENSION Y(NEQN),WT(NEQN),PHI(NEQN,16),P(NEQN),YP(NEQN),PSI(12),
; C    1  ALPHA(12),BETA(12),SIG(13),V(12),W(12),G(13),GI(11),IV(10),
; C    2  RPAR(*),IPAR(*)
; C
; C    **Note**
; C
; C   The user must also declare  START ,  CRASH ,  PHASE1  and  NORND
; C   logical variables and  DF  an EXTERNAL subroutine, supply the
; C   subroutine  DF(X,Y,YP)  to evaluate
; C      DY(I)/DX = YP(I) = DF(X,Y(1),Y(2),...,Y(NEQN))
; C   and initialize only the following parameters.
; C      NEQN -- number of equations to be integrated
; C      Y(*) -- vector of initial values of dependent variables
; C      X -- initial value of the independent variable
; C      H -- nominal step size indicating direction of integration
; C           and maximum size of step.  Must be variable
; C      EPS -- local error tolerance per step.  Must be variable
; C      WT(*) -- vector of non-zero weights for error criterion
; C      START -- .TRUE.
; C      YP(*) -- vector of initial derivative values
; C      KSTEPS -- set KSTEPS to zero
; C      TWOU -- 2.*U where U is machine unit roundoff quantity
; C      FOURU -- 4.*U where U is machine unit roundoff quantity
; C   Define U to be the machine unit roundoff quantity by calling
; C   the function routine  D1MACH,  U = D1MACH(4), or by
; C   computing U so that U is the smallest positive number such
; C   that 1.0+U .GT. 1.0.
; C
; C   DDEABM_DSTEPS  requires that the L2 norm of the vector with components
; C   LOCAL ERROR(L)/WT(L)  be less than  EPS  for a successful step.  The
; C   array  WT  allows the user to specify an error test appropriate
; C   for his problem.  For example,
; C      WT(L) = 1.0  specifies absolute error,
; C            = ABS(Y(L))  error relative to the most recent value of the
; C                 L-th component of the solution,
; C            = ABS(YP(L))  error relative to the most recent value of
; C                 the L-th component of the derivative,
; C            = MAX(WT(L),ABS(Y(L)))  error relative to the largest
; C                 magnitude of L-th component obtained so far,
; C            = ABS(Y(L))*RELERR/EPS + ABSERR/EPS  specifies a mixed
; C                 relative-absolute test where  RELERR  is relative
; C                 error,  ABSERR  is absolute error and  EPS =
; C                 MAX(RELERR,ABSERR) .
; C
; C      Subsequent calls --
; C
; C   Subroutine  DDEABM_DSTEPS  is designed so that all information needed to
; C   continue the integration, including the step size  H  and the order
; C   K , is returned with each step.  With the exception of the step
; C   size, the error tolerance, and the weights, none of the parameters
; C   should be altered.  The array  WT  must be updated after each step
; C   to maintain relative error tests like those above.  Normally the
; C   integration is continued just beyond the desired endpoint and the
; C   solution interpolated there with subroutine  SINTRP .  If it is
; C   impossible to integrate beyond the endpoint, the step size may be
; C   reduced to hit the endpoint since the code will not take a step
; C   larger than the  H  input.  Changing the direction of integration,
; C   i.e., the sign of  H , requires the user set  START = .TRUE. before
; C   calling  DDEABM_DSTEPS  again.  This is the only situation in which  START
; C   should be altered.
; C
; C   Output from DDEABM_DSTEPS
; C
; C      Successful Step --
; C
; C   The subroutine returns after each successful step with  START  and
; C   CRASH  set .FALSE. .  X  represents the independent variable
; C   advanced one step of length  HOLD  from its value on input and  Y
; C   the solution vector at the new value of  X .  All other parameters
; C   represent information corresponding to the new  X  needed to
; C   continue the integration.
; C
; C      Unsuccessful Step --
; C
; C   When the error tolerance is too small for the machine precision,
; C   the subroutine returns without taking a step and  CRASH = .TRUE. .
; C   An appropriate step size and error tolerance for continuing are
; C   estimated and all other information is restored as upon input
; C   before returning.  To continue with the larger tolerance, the user
; C   just calls the code again.  A restart is neither required nor
; C   desirable.
; C
; C***REFERENCES  L. F. Shampine and M. K. Gordon, Solving ordinary
; C                 differential equations with ODE, STEP, and INTRP,
; C                 Report SLA-73-1060, Sandia Laboratories, 1973.
; C***ROUTINES CALLED  D1MACH, DDEABM_DHSTRT
; C***REVISION HISTORY  (YYMMDD)
; C   740101  DATE WRITTEN
; C   890531  Changed all specific intrinsics to generic.  (WRB)
; C   890831  Modified array declarations.  (WRB)
; C   890831  REVISION DATE from Version 3.2
; C   891214  Prologue converted to Version 4.0 format.  (BAB)
; C   920501  Reformatted the REFERENCES section.  (WRB)
; C***END PROLOGUE  DDEABM_DSTEPS
; C
;       INTEGER I, IFAIL, IM1, IP1, IPAR, IQ, J, K, KM1, KM2, KNEW,
;      1      KOLD, KP1, KP2, KSTEPS, L, LIMIT1, LIMIT2, NEQN, NS, NSM2,
;      2      NSP1, NSP2
;       DOUBLE PRECISION ABSH, ALPHA, BETA, BIG, D1MACH,
;      1      EPS, ERK, ERKM1, ERKM2, ERKP1, ERR,
;      2      FOURU, G, GI, GSTR, H, HNEW, HOLD, P, P5EPS, PHI, PSI, R,
;      3      REALI, REALNS, RHO, ROUND, RPAR, SIG, TAU, TEMP1,
;      4      TEMP2, TEMP3, TEMP4, TEMP5, TEMP6, TWO, TWOU, U, V, W, WT,
;      5      X, XOLD, Y, YP
;       LOGICAL START,CRASH,PHASE1,NORND
;       DIMENSION Y(*),WT(*),PHI(NEQN,16),P(*),YP(*),PSI(12),
;      1  ALPHA(12),BETA(12),SIG(13),V(12),W(12),G(13),GI(11),IV(10),
;      2  RPAR(*),IPAR(*)
;       DIMENSION TWO(13),GSTR(13)
;       EXTERNAL DF
;       SAVE TWO, GSTR

;       DATA TWO(1),TWO(2),TWO(3),TWO(4),TWO(5),TWO(6),TWO(7),TWO(8),
;      1     TWO(9),TWO(10),TWO(11),TWO(12),TWO(13)
;      2     /2.0D0,4.0D0,8.0D0,16.0D0,32.0D0,64.0D0,128.0D0,256.0D0,
;      3      512.0D0,1024.0D0,2048.0D0,4096.0D0,8192.0D0/
;       DATA GSTR(1),GSTR(2),GSTR(3),GSTR(4),GSTR(5),GSTR(6),GSTR(7),
;      1     GSTR(8),GSTR(9),GSTR(10),GSTR(11),GSTR(12),GSTR(13)
;      2     /0.5D0,0.0833D0,0.0417D0,0.0264D0,0.0188D0,0.0143D0,0.0114D0,
;      3      0.00936D0,0.00789D0,0.00679D0,0.00592D0,0.00524D0,0.00468D0/

    common ddeabm_func_common
    TWO = 2d^(dindgen(13)+1)
    GSTR = [ 0.5D0,0.0833D0,0.0417D0,0.0264D0,0.0188D0,0.0143D0,0.0114D0, $
             0.00936D0,0.00789D0,0.00679D0,0.00592D0,0.00524D0,0.00468D0 ]

; C
; C       ***     BEGIN BLOCK 0     ***
; C   CHECK IF STEP SIZE OR ERROR TOLERANCE IS TOO SMALL FOR MACHINE
; C   PRECISION.  IF FIRST STEP, INITIALIZE PHI ARRAY AND ESTIMATE A
; C   STARTING STEP SIZE.
; C                   ***
; C
; C   IF STEP SIZE IS TOO SMALL, DETERMINE AN ACCEPTABLE ONE
; C
; C***FIRST EXECUTABLE STATEMENT  DDEABM_DSTEPS
      CRASH = 1L
      IF (ABS(H) LT FOURU*ABS(X)) THEN BEGIN
          H = (FOURU*ABS(X)) * ( (H GE 0)?(+1):(-1) )
          RETURN
      ENDIF

      P5EPS = 0.5D0*EPS
; C
; C   IF ERROR TOLERANCE IS TOO SMALL, INCREASE IT TO AN ACCEPTABLE VALUE
; C
      ROUND = TOTAL( (Y/WT)^2 )
      ROUND = TWOU*SQRT(ROUND)

      IF (P5EPS LT ROUND) THEN BEGIN
          EPS = 2.0D0*ROUND*(1.0D0 + FOURU)
          RETURN
      ENDIF

      CRASH = 0L
      G(0) = 1.0D0
      G(1) = 0.5D0
      SIG(0) = 1.0D0
      IF (NOT START) THEN GOTO, DSTEPS_99
; C
; C   INITIALIZE.  COMPUTE APPROPRIATE STEP SIZE FOR FIRST STEP
; C
; C     CALL DF(X,Y,YP,RPAR,IPAR)
; C     SUM = 0.0
      PHI(*,0) = YP
      PHI(*,1) = 0
; C20     SUM = SUM + (YP(L-1)/WT(L-1))**2
; C     SUM = SQRT(SUM)
; C     ABSH = ABS(H)
; C     IF(EPS .LT. 16.0*SUM*H*H) ABSH = 0.25*SQRT(EPS/SUM)
; C     H = SIGN(MAX(ABSH,FOURU*ABS(X)),H)
; C
      U = (machar(/double)).eps          ;; XXX
      BIG = SQRT((machar(/double)).xmax) ;; XXX


      ;; Save and restore values from PHI
      phi3 = phi(*,2) & phi4 = phi(*,3)
      phi5 = phi(*,4) & phi6 = phi(*,5)
      DDEABM_DHSTRT, DF,NEQN,X,X+H,Y,YP,WT,1,U,BIG, $
        phi3, phi4, phi5, phi6, private, fa, h, dfname
      phi(*,2) = phi3 & phi(*,3) = phi4 
      phi(*,4) = phi5 & phi(*,5) = phi6
      if ddeabm_funcerror NE 0 then return

      HOLD = 0.0D0
      K = 1L
      KOLD = 0L
      KPREV = 0L
      START = 0L
      PHASE1 = 1L
      NORND = 1L
      IF (P5EPS LE 100.0D0*ROUND) THEN BEGIN
          NORND = 0L
          PHI(*,14) = 0
      ENDIF

      DSTEPS_99:
      IFAIL = 0L
; C       ***     END BLOCK 0     ***
; C
; C       ***     BEGIN BLOCK 1     ***
; C   COMPUTE COEFFICIENTS OF FORMULAS FOR THIS STEP.  AVOID COMPUTING
; C   THOSE QUANTITIES NOT CHANGED WHEN STEP SIZE IS NOT CHANGED.
; C                   ***
; C
      DSTEPS_100:
      KP1 = K+1
      KP2 = K+2
      KM1 = K-1
      KM2 = K-2
; C
; C   NS IS THE NUMBER OF DSTEPS TAKEN WITH SIZE H, INCLUDING THE CURRENT
; C   ONE.  WHEN K.LT.NS, NO COEFFICIENTS CHANGE
; C
      IF (H NE HOLD) THEN NS = 0L
      IF (NS LE KOLD) THEN NS = NS+1
      NSP1 = NS+1
      IF (K LT NS) THEN GOTO, DSTEPS_199
; C
; C   COMPUTE THOSE COMPONENTS OF ALPHA(*),BETA(*),PSI(*),SIG(*) WHICH
; C   ARE CHANGED
; C
      BETA(NS-1) = 1.0D0
      ALPHA(NS-1) = 1.0D0/NS
      TEMP1 = H*NS
      SIG(NSP1-1) = 1.0D0
      IF (K GE NSP1) THEN BEGIN
          FOR I = NSP1, K DO BEGIN
              IM1 = I-1-1  ;; *** Note IM1-1 here!
              II = I-1
              TEMP2 = PSI(IM1)
              PSI(IM1) = TEMP1
              BETA(II) = BETA(IM1)*PSI(IM1)/TEMP2
              TEMP1 = TEMP2 + H
              ALPHA(II) = H/TEMP1
              SIG(I) = I*ALPHA(II)*SIG(II)
          ENDFOR
      ENDIF
      PSI(K-1) = TEMP1
; C
; C   COMPUTE COEFFICIENTS G(*)
; C
; C   INITIALIZE V(*) AND SET W(*).
; C
      IF (NS LE 1) THEN BEGIN
          KK = dindgen(K)+1
          V(0:K-1) = 1.0D0/(KK*(KK+1))
          W(0:K-1) = V(0:K-1)
          IVC = 0L
          KGI = 0L
          IF (K EQ 1) THEN GOTO, DSTEPS_140
          KGI = 1L
          GI(0) = W(1)
          GOTO, DSTEPS_140
      ENDIF
; C
; C   IF ORDER WAS RAISED, UPDATE DIAGONAL PART OF V(*)
; C
      IF (K LE KPREV) THEN GOTO, DSTEPS_130
      IF (IVC NE 0) THEN BEGIN 
          JV = KP1 - IV(IVC-1)
          IVC = IVC - 1
      ENDIF ELSE BEGIN
          JV = 1L
          TEMP4 = K*KP1
          V(K-1) = 1.0D0/TEMP4
          W(K-1) = V(K-1)
          IF (K EQ 2) THEN BEGIN
              KGI = 1L
              GI(0) = W(1)
          ENDIF
      ENDELSE
      NSM2 = NS-2
      IF (NSM2 GE JV) THEN BEGIN
          
          FOR J = JV, NSM2 DO BEGIN
              I = K-J-1  ;; *** NOTE: I-1 here!
              V(I) = V(I) - ALPHA(J)*V(I+1)
              W(I) = V(I)
          ENDFOR

          IF (I EQ 2) THEN BEGIN
              KGI = NS - 1
              GI(KGI-1) = W(1)
          ENDIF
      ENDIF
; C
; C   UPDATE V(*) AND SET W(*)
; C
      DSTEPS_130:
      LIMIT1 = KP1 - NS
      TEMP5 = ALPHA(NS-1)
      V(0:LIMIT1-1) = V(0:LIMIT1-1) - TEMP5*V(1:LIMIT1)
      W(0:LIMIT1-1) = V(0:LIMIT1-1)
      G(NSP1-1) = W(0)
      IF (LIMIT1 NE 1) THEN BEGIN
          KGI = NS
          GI(KGI-1) = W(1)
      ENDIF
      W(LIMIT1) = V(LIMIT1)
      IF (K LT KOLD) THEN BEGIN 
          IVC = IVC + 1
          IV(IVC-1) = LIMIT1 + 2
      ENDIF
; C
; C   COMPUTE THE G(*) IN THE WORK VECTOR W(*)
; C
      DSTEPS_140:
      NSP2 = NS + 2
      KPREV = K
      IF (KP1 GE NSP2) THEN BEGIN
          FOR I = NSP2, KP1 DO BEGIN
              LIMIT2 = KP2 - I
              TEMP6 = ALPHA(I-2)
              W(0:LIMIT2-1) = W(0:LIMIT2-1) - TEMP6*W(1:LIMIT2)
              G(I-1) = W(0)
          ENDFOR
      ENDIF

      DSTEPS_199:
; C       ***     END BLOCK 1     ***
; C
; C       ***     BEGIN BLOCK 2     ***
; C   PREDICT A SOLUTION P(*), EVALUATE DERIVATIVES USING PREDICTED
; C   SOLUTION, ESTIMATE LOCAL ERROR AT ORDER K AND ERRORS AT ORDERS K,
; C   K-1, K-2 AS IF CONSTANT STEP SIZE WERE USED.
; C                   ***
; C
; C   INCREMENT COUNTER ON ATTEMPTED DSTEPS
; C
      KSTEPS = KSTEPS + 1
; C
; C   CHANGE PHI TO PHI STAR
; C
      IF (K GE NSP1) THEN BEGIN
          FOR I = NSP1, K DO BEGIN
              TEMP1 = BETA(I-1)
              PHI(*,I-1) = TEMP1*PHI(*,I-1)
          ENDFOR
      ENDIF
; C
; C   PREDICT SOLUTION AND DIFFERENCES
; C
      PHI(*,KP2-1) = PHI(*,KP1-1)
      PHI(*,KP1-1) = 0
      P(*) = 0
      FOR J = 1L, K DO BEGIN
        I = KP1 - J - 1 ;; *** NOTE: I-1 here!
        TEMP2 = G(I)
        P = P + TEMP2*PHI(*,I)
        PHI(*,I) = PHI(*,I) + PHI(*,I+1)
      ENDFOR
      IF NOT (NORND) THEN BEGIN
          TAU = H*P - PHI(*,14)
          P = Y + TAU
          PHI(*,15) = (P - Y) - TAU
      ENDIF ELSE BEGIN
          P = Y + H*P
      ENDELSE
      XOLD = X

      X = X + H
      ABSH = ABS(H)

      YP = CALL_FUNCTION(DFNAME, DF, X, P, PRIVATE, _EXTRA=FA)
      if ddeabm_funcerror NE 0 then return

; C
; C   ESTIMATE ERRORS AT ORDERS K,K-1,K-2
; C
      ERKM2 = 0.0D0
      ERKM1 = 0.0D0
      TEMP3 = 1.0D0/WT
      TEMP4 = YP - PHI(*,0)
      ERK = total((temp4*temp3)^2)
      IF (KM2 GT 0) THEN $
        ERKM2 = TOTAL( ((PHI(*,KM1-1)+TEMP4)*TEMP3)^2 )
      IF (KM2 GE 0) THEN $
        ERKM1 = TOTAL( ((PHI(*,K-1)+TEMP4)*TEMP3)^2 )

      IF (KM2 GT 0) THEN $
        ERKM2 = ABSH*SIG(KM1-1)*GSTR(KM2-1)*SQRT(ERKM2)
      IF (KM2 GE 0) THEN $
        ERKM1 = ABSH*SIG(K-1)*GSTR(KM1-1)*SQRT(ERKM1)
      TEMP5 = ABSH*SQRT(ERK)
      ERR = TEMP5*(G(K-1)-G(KP1-1))
      ERK = TEMP5*SIG(KP1-1)*GSTR(K-1)
      KNEW = K
; C
; C   TEST IF ORDER SHOULD BE LOWERED
; C
      IF (KM2 GT 0) THEN BEGIN
          IF(MAX([ERKM1,ERKM2]) LE ERK) THEN KNEW = KM1
      ENDIF ELSE IF (KM2 EQ 0) THEN BEGIN
          IF(ERKM1 LE 0.5D0*ERK) THEN KNEW = KM1
      ENDIF
; C
; C   TEST IF STEP SUCCESSFUL
; C
      IF (ERR LE EPS) THEN GOTO, DSTEPS_400
; C       ***     END BLOCK 2     ***
; C
; C       ***     BEGIN BLOCK 3     ***
; C   THE STEP IS UNSUCCESSFUL.  RESTORE  X, PHI(*,*), PSI(*) .
; C   IF THIRD CONSECUTIVE FAILURE, SET ORDER TO ONE.  IF STEP FAILS MORE
; C   THAN THREE TIMES, CONSIDER AN OPTIMAL STEP SIZE.  DOUBLE ERROR
; C   TOLERANCE AND RETURN IF ESTIMATED STEP SIZE IS TOO SMALL FOR MACHINE
; C   PRECISION.
; C                   ***
; C
; C   RESTORE X, PHI(*,*) AND PSI(*)
; C
      PHASE1 = 0L
      X = XOLD
      FOR I = 0L, K-1 DO BEGIN
        PHI(*,I) = (PHI(*,I) - PHI(*,I+1))/BETA(I)
      ENDFOR
      IF (K GE 2) THEN BEGIN 
          PSI(0:K-2) = PSI(1:K-1) - H
      ENDIF

; C
; C   ON THIRD FAILURE, SET ORDER TO ONE.  THEREAFTER, USE OPTIMAL STEP
; C   SIZE
; C
      IFAIL = IFAIL + 1
      TEMP2 = 0.5D0
      IF (IFAIL - 3 GT 0) THEN BEGIN
          IF (P5EPS LT 0.25D0*ERK) THEN TEMP2 = SQRT(P5EPS/ERK)
      ENDIF
      IF (IFAIL - 3 GE 0) THEN KNEW = 1L
      H = TEMP2*H
      K = KNEW
      NS = 0L
      IF (ABS(H) LT FOURU*ABS(X)) THEN BEGIN
          CRASH = 1L
          H = (FOURU*ABS(X))*( (H GE 0)?(+1):(-1) )
          EPS = EPS + EPS
          RETURN
      ENDIF

      GOTO, DSTEPS_100
; C       ***     END BLOCK 3     ***
; C
; C       ***     BEGIN BLOCK 4     ***
; C   THE STEP IS SUCCESSFUL.  CORRECT THE PREDICTED SOLUTION, EVALUATE
; C   THE DERIVATIVES USING THE CORRECTED SOLUTION AND UPDATE THE
; C   DIFFERENCES.  DETERMINE BEST ORDER AND STEP SIZE FOR NEXT STEP.
; C                   ***
      DSTEPS_400:
      KOLD = K
      HOLD = H
; C
; C   CORRECT AND EVALUATE
; C
      TEMP1 = H*G(KP1-1)
      IF NOT (NORND) THEN BEGIN
          
          TEMP3 = Y
          RHO = TEMP1*(YP - PHI(*,0)) - PHI(*,15)
          Y = P + RHO
          PHI(*,14) = (Y - P) - RHO
          P = TEMP3
      ENDIF ELSE BEGIN
          TEMP3 = Y
          Y = P + TEMP1*(YP - PHI(*,0))
          P = TEMP3
      ENDELSE

      YP = CALL_FUNCTION(DFNAME, DF, X, Y, PRIVATE, _EXTRA=FA)
      if ddeabm_funcerror NE 0 then return


; C
; C   UPDATE DIFFERENCES FOR NEXT STEP
; C
      PHI(*,KP1-1) = YP - PHI(*,0)
      PHI(*,KP2-1) = PHI(*,KP1-1) - PHI(*,KP2-1)
      FOR I = 0L, K-1 DO BEGIN
          PHI(*,I) = PHI(*,I) + PHI(*,KP1-1)
      ENDFOR
; C
; C   ESTIMATE ERROR AT ORDER K+1 UNLESS:
; C     IN FIRST PHASE WHEN ALWAYS RAISE ORDER,
; C     ALREADY DECIDED TO LOWER ORDER,
; C     STEP SIZE NOT CONSTANT SO ESTIMATE UNRELIABLE
; C
      ERKP1 = 0.0D0
      IF (KNEW EQ KM1  OR  K EQ 12) THEN PHASE1 = 0L
      IF (PHASE1) THEN GOTO, DSTEPS_450
      IF (KNEW EQ KM1) THEN GOTO, DSTEPS_455
      IF (KP1 GT NS) THEN GOTO, DSTEPS_460
      ERKP1 = TOTAL( (PHI(*,KP2-1)/WT)^2 )
      ERKP1 = ABSH*GSTR(KP1-1)*SQRT(ERKP1)
; C
; C   USING ESTIMATED ERROR AT ORDER K+1, DETERMINE APPROPRIATE ORDER
; C   FOR NEXT STEP
; C
      IF (K LE 1) THEN BEGIN
          IF (ERKP1 GE 0.5D0*ERK) THEN GOTO, DSTEPS_460
      ENDIF ELSE BEGIN
          IF (ERKM1 LE MIN([ERK,ERKP1])) THEN GOTO, DSTEPS_455
          IF (ERKP1 GE ERK  OR  K EQ 12) THEN GOTO, DSTEPS_460
      ENDELSE
; C
; C   HERE ERKP1 .LT. ERK .LT. MAX(ERKM1,ERKM2) ELSE ORDER WOULD HAVE
; C   BEEN LOWERED IN BLOCK 2.  THUS ORDER IS TO BE RAISED
; C
; C   RAISE ORDER
; C
      DSTEPS_450:
      K = KP1
      ERK = ERKP1
      GOTO, DSTEPS_460
; C
; C   LOWER ORDER
; C
      DSTEPS_455:
      K = KM1
      ERK = ERKM1
; C
; C   WITH NEW ORDER DETERMINE APPROPRIATE STEP SIZE FOR NEXT STEP
; C
      DSTEPS_460:
      HNEW = H + H
      IF NOT ( (PHASE1) OR $
               (P5EPS GE ERK*TWO(K)) ) THEN BEGIN
          HNEW = H
          IF (P5EPS LT ERK) THEN BEGIN
              TEMP2 = K+1
              R = (P5EPS/ERK)^(1.0D0/TEMP2)
              HNEW = ABSH*MAX([0.5D0,MIN([0.9D0,R])])
              HNEW = MAX([HNEW,FOURU*ABS(X)])
              HNEW = (H GE 0)?(+HNEW):(-HNEW)
          ENDIF
      ENDIF

      if n_elements(max_stepsize) GT 0 then begin
          HNEW = HNEW < max_stepsize(0) > (-max_stepsize(0))
      endif

      H = HNEW
      RETURN
; C       ***     END BLOCK 4     ***
  END


; ------------------------------------------------------------------------

pro DDEABM, DF, T, Y, TOUT0, PRIVATE, FUNCTARGS=fa, STATE=state, $
            CONTROL=control, $
            init=init0, intermediate=intermediate, tstop=TSTOP0, $
            epsrel=RTOL, epsabs=ATOL, status=IDID, $
            TGRID=tgrid, YGRID=ygrid, YPGRID=ypgrid, $
            NGRID=ngrid0, NOUTGRID=nsamp, $
            TIMPULSE=timpulse, YIMPULSE=yimpulse, $
            MAX_STEPSIZE=max_stepsize, $
            NFEV=nfev, errmsg=errmsg, dostatusline=dostatusline

      common ddeabm_func_common, ddeabm_nfev, ddeabm_funcerror

      IDID = -33
      errmsg = ''
      if n_params() EQ 0 then begin
          message, 'USAGE:', /info
          message, '  DDEABM, FUNCNAME, T0, Y0, TOUT, STATE, PRIVATE, '+$
            'FUNCTARGS=fa, INIT=init, [EPSREL=epsrel, EPSABS=epsabs, '+$
            'STATUS=status, /INTERMEDIATE, ...]', /info
          return
      endif

; C***FIRST EXECUTABLE STATEMENT  DDEABM
      NEQ = N_ELEMENTS(Y)

      IF NEQ LT 1 THEN BEGIN
          errmsg = 'The number of equations, NEQ, must be greater than '+$
            'or equal to 1'
          idid = -33L
          RETURN
      ENDIF

      ;; Initialize the number of function evaluations
      ddeabm_nfev = 0L
      nfev = 0L

      ;; Construct the wrapper function to be used
      dfname = 'ddeabm_func'
      dfname = dfname + ((n_elements(private) GT 0)?'1':'0')
      dfname = dfname + ((n_elements(fa) GT 0     )?'e':'n')
      
      ;; If either of the tolerances are undefined, then define with
      ;; default tolerances and with the same number of elements
      if n_elements(rtol) EQ 0 AND n_elements(atol) EQ 0 then begin
          rtol = 1d-4
          atol = 1d-6
      endif else if n_elements(rtol) GT 0 AND n_elements(atol) EQ 0 then begin
          atol = rtol*0d
      endif else if n_elements(atol) GT 0 AND n_elements(rtol) EQ 0 then begin
          rtol = atol*0D
      endif

      ;; Compare to be sure the same number 
      if n_elements(rtol) NE n_elements(atol) $
        OR (n_elements(y) NE n_elements(rtol) AND n_elements(rtol) NE 1) $
        then begin
          errmsg = 'The number of absolute and relative tolerance values '+ $
            'must match the number of equations being solved.'
          idid = -33L
      endif

      ;; Be sure to initialize if there is no state variable
      ;; NOTE: DDEABM uses INIT=0 to mean initialize; INIT=1 means
      ;;       don't initialize, which is the opposite sense from the
      ;;       input keyword.
      userinit = 1-keyword_set(init0)
      if n_elements(state) EQ 0 then userinit = 0L

      ;; Construct the INFO array from keywords
      INFO = [ userinit, n_elements(rtol) GT 1, $
               keyword_set(intermediate), n_elements(tstop0) GT 0]

      ;; Construct the STATE array if this is the first pass
      IF ( INFO(1-1) EQ 0 ) OR N_ELEMENTS(STATE) EQ 0 THEN BEGIN
          STATE = {YPOUT: dblarr(neq), TSTAR: 0D, YP: dblarr(NEQ), $
                   YY: dblarr(NEQ), WT: dblarr(NEQ), P: dblarr(NEQ), $
                   PHI: dblarr(NEQ,16), ALPHA: dblarr(12), BETA: dblarr(12), $
                   PSI: dblarr(12), V: dblarr(12), W: dblarr(12), $
                   SIG: dblarr(13), G: dblarr(13), GI: dblarr(11), $
                   XOLD: 0D, HOLD: 0D, TOLD: 0D, DELSN: 0D, TWOU: 0D, $
                   FOURU: 0D, H: 0D, EPS: 0D, X: 0D, TSTOP: 0D, $
                   START: 0L, PHASE1: 0L, NORND: 0L, STIFF: 0L, $
                   INTOUT: 0L, NS: 0L, KORD: 0L, KOLD: 0L, INTERNAL_INIT: 0L, $
                   KSTEPS: 0L, KLE4: 0L, IQUIT: 0L, KPREV: 0L, IVC: 0L, $
                   IV: lonarr(10), KGI: 0L, NEQ: NEQ, COUNT: 0L}
      ENDIF

      if n_elements(tstop0) GT 0 then $
        state.tstop = tstop0(0)

      IF (STATE.COUNT GE 5) THEN BEGIN
         IF (T EQ STATE.TSTAR) THEN BEGIN
             errmsg = 'AN APPARENT INFINITE LOOP HAS BEEN DETECTED.  '+ $
               'YOU HAVE MADE REPEATED CALLS AT T = '+strtrim(t,2)+ $
               ' AND THE INTEGRATION HAS NOT ADVANCED.  CHECK THE '+ $
               'WAY YOU HAVE SET PARAMETERS FOR THE CALL TO THE '+ $
               'CODE, PARTICULARLY INFO(1-1).'
             RETURN
         ENDIF
     ENDIF

     IF NEQ NE STATE.NEQ THEN BEGIN
         errmsg = 'You have initialized DDEABM with a different number '+$
           'of equations, NEQ, than this call has provided.'
         RETURN
     ENDIF
; C
; C     CHECK LRW AND LIW FOR SUFFICIENT STORAGE ALLOCATION
; C
      IDID=0L
; C
; C     COMPUTE THE INDICES FOR THE ARRAYS TO BE STORED IN THE WORK ARRAY
; C
     YPOUT = STATE.YPOUT
     TSTAR = STATE.TSTAR
     YP    = STATE.YP   
     YY    = STATE.YY   
     WT    = STATE.WT   
     P     = STATE.P    
     PHI   = STATE.PHI  
     ALPHA = STATE.ALPHA
     BETA  = STATE.BETA 
     PSI   = STATE.PSI  
     V     = STATE.V    
     W     = STATE.W    
     SIG   = STATE.SIG  
     G     = STATE.G    
     GI    = STATE.GI   
     XOLD  = STATE.XOLD 
     HOLD  = STATE.HOLD 
     TOLD  = STATE.TOLD 
     DELSN = STATE.DELSN
     TWOU  = STATE.TWOU 
     FOURU = STATE.FOURU

     H     = STATE.H
     EPS   = STATE.EPS
     X     = STATE.X
     TSTOP = STATE.TSTOP


     STATE.TSTAR = T

      IF (INFO(1-1) NE 0) THEN BEGIN
          START = STATE.START
          PHASE1 = STATE.PHASE1
          NORND = STATE.NORND
          STIFF = STATE.STIFF
          INTOUT = STATE.INTOUT
      ENDIF
      
     NS     = STATE.NS
     KORD   = STATE.KORD
     KOLD   = STATE.KOLD
     INTERNAL_INIT   = STATE.INTERNAL_INIT
     KSTEPS = STATE.KSTEPS
     KLE4   = STATE.KLE4
     IQUIT  = STATE.IQUIT
     KPREV  = STATE.KPREV
     IVC    = STATE.IVC
     IV     = STATE.IV
     KGI    = STATE.KGI

     if n_elements(ngrid0) GT 0 then begin
         if NOT keyword_set(intermediate) then begin
             errmsg = 'ERROR: NGRID and /INTERMEDIATE must be specified '+$
               'together'
             return
         endif
         ngrid = round(ngrid0(0))
     endif else begin
         ngrid = n_elements(tout0)
     endelse

     tgrid = dblarr(ngrid)
     ygrid = dblarr(neq, ngrid)
     ypgrid = dblarr(neq, ngrid)
     forward = tout0(0) GT t     ;; 1=FORWARD; 0=BACKWARD

     ki = 1L
     nimpulse = n_elements(timpulse)
     if nimpulse GT 0 then begin
         if nimpulse NE n_elements(yimpulse)/neq then begin
             errmsg = 'ERROR: TIMPULSE and YIMPULSE must have the same '+$
               'number of samples'
             return
         endif

         if forward then begin
             wh = where(timpulse GT tout0(0), ct)
             if ct EQ 0 then ki = nimpulse else ki = min(wh)
         endif else begin
             wh = where(timpulse LT tout0(0), ct)
             if ct EQ 0 then ki = 0L else ki = max(wh)
         endelse

     endif

     ;; Initialize the user function
     if info(1-1) EQ 0 AND keyword_set(control) then begin
         ddeabm_funcerror = call_function(dfname, df, t, $
                                          CONTROL={message: 'INITIALIZE'}, $
                                          y, private, _EXTRA=fa)
         if ddeabm_funcerror LT 0 then begin
             errmsg = 'ERROR: user function failed to initialize'
             goto, FINISH_INTEGRATION
         endif
     endif


     i = 0L   ;; Output grid position counter
     nsamp = 0L
     while (i LT ngrid) do begin
         if keyword_set(dostatusline) then $
           statusline, string(i, ngrid, format='(I8,"/",I8)'), 0, /left
         doimpulse = 0  ;; Signal to process an impulse (0=no; 1=yes; 2=both)
         if keyword_set(intermediate) then begin
             TOUT = TOUT0(0)
         endif else begin
             TOUT = TOUT0(i)

             if (ki GE 0) AND (ki LT nimpulse) then begin
                 if abs(timpulse(ki)-t) LE abs(tout-t) then begin
                     doimpulse = 1
                     if TIMPULSE(ki) EQ TOUT then doimpulse = 2
                     TOUT = TIMPULSE(ki)
                 endif
             endif
         endelse

         ddeabm_funcerror = 0
         DDEABM_DDES, DF,NEQ,T,Y,TOUT,INFO,RTOL,ATOL,IDID,YPOUT, $
           YP,YY,WT,P,PHI, $
           ALPHA,BETA,PSI,V, $
           W,SIG,G,GI,H, $
           EPS,X,XOLD,HOLD, $
           TOLD,DELSN,TSTOP,TWOU, $
           FOURU,START,PHASE1,NORND,STIFF,INTOUT, NS, KORD, KOLD, INTERNAL_INIT, $
           KSTEPS, KLE4, IQUIT, KPREV, IVC, IV, KGI, PRIVATE, FA, dfname, $
           ERRMSG=errmsg, max_stepsize=max_stepsize

         if ddeabm_funcerror NE 0 then begin
             case ddeabm_funcerror of
                 -16:  errmsg = 'ERROR: user function returned non-finite values'
                 else: errmsg = 'ERROR: unknown internal error occurred'
             endcase
             goto, FINISH_INTEGRATION
         endif

         if IDID GT 0 then begin

             ;; === Store the result
             if (doimpulse EQ 0) OR (doimpulse EQ 2) then begin
                 ;; This was not an impulse-only stopping point...
                 ;; Store the result
                 tgrid(i) = T
                 ygrid(*,i) = Y
                 ypgrid(*,i) = YPOUT
                 nsamp = nsamp + 1
                 i = i + 1
             endif

             ;; === Handle any impulse changes
             if (doimpulse GT 0) then begin    
                 ;; Apply an impulse
                 if forward then begin
                     Y = Y + YIMPULSE(*,ki)
                     ki = ki + 1
                 endif else begin
                     Y = Y - YIMPULSE(*,ki)
                     ki = ki - 1
                 endelse

                 ;; Special case: the same TOUT can be listed twice
                 ;;   for before and after an impulse.  In that case
                 ;;   store the same values after we have incremented.
                 if (doimpulse EQ 2) then if (TOUT0(i) EQ TOUT0(i-1)) then begin
                     tgrid(i) = T
                     ygrid(*,i) = Y
                     ;; We have to re-call the function since we
                     ;; crossed the discontinuity.
                     ypgrid(*,i) = call_function(dfname,df, t, y, private,_EXTRA=fa)
                     nsamp = nsamp + 1
                     i = i + 1
                 endif

                 doimpulse = 0
                 info(1-1) = 0L  ;; NOTE that INIT=0 means initialize!!!
             endif

             ;; Reset KSTEPS since we successfully integrated this step
             ksteps = 0L

             ;; End if we reach the stopping point early
             if keyword_set(intermediate) then begin
                 if idid EQ 2 OR idid EQ 3 OR t EQ tout(0) then $
                   goto, FINISH_INTEGRATION
             endif

         endif else begin
             goto, FINISH_INTEGRATION
         endelse

     endwhile

     FINISH_INTEGRATION:

     STATE.YPOUT = YPOUT
     STATE.YP    = YP   
     STATE.YY    = YY   
     STATE.WT    = WT   
     STATE.P     = P    
     STATE.PHI   = PHI  
     STATE.ALPHA = ALPHA
     STATE.BETA  = BETA 
     STATE.PSI   = PSI  
     STATE.V     = V    
     STATE.W     = W    
     STATE.SIG   = SIG  
     STATE.G     = G    
     STATE.GI    = GI   
     STATE.XOLD  = XOLD 
     STATE.HOLD  = HOLD 
     STATE.TOLD  = TOLD 
     STATE.DELSN = DELSN
     STATE.TWOU  = TWOU 
     STATE.FOURU = FOURU
      
     STATE.H     = H
     STATE.EPS   = EPS
     STATE.X     = X
     STATE.TSTOP = TSTOP

     STATE.NS     = NS
     STATE.KORD   = KORD
     STATE.KOLD   = KOLD
     STATE.INTERNAL_INIT   = INTERNAL_INIT
     STATE.KSTEPS = KSTEPS
     STATE.KLE4   = KLE4
     STATE.IQUIT  = IQUIT
     STATE.KPREV  = KPREV
     STATE.IVC    = IVC
     STATE.IV     = IV
     STATE.KGI    = KGI

     STATE.START  = START
     STATE.PHASE1 = PHASE1
     STATE.NORND  = NORND
     STATE.STIFF  = STIFF
     STATE.INTOUT = INTOUT
     NFEV = DDEABM_NFEV

     ;; Pass back INIT to user, and remember to invert the sense
     ;; between the internal variable and the external variable.
     INIT0 = (info(0) EQ 0)
     TSTOP0 = STATE.TSTOP
     ;; XXX what to do about interrupted case where INFO(1-1) is
     ;; negative, and the user must reset it?
     ;; Answer: Add RESUME keyword, and enforce the behavior that
     ;; RESUME must only be set after an interruption.  Must save one
     ;; more variable in STATE with the previous INFO(1-1) value.

      IF (IDID NE (-2)) THEN STATE.COUNT = STATE.COUNT + 1L
      IF (T NE STATE.TSTAR) THEN STATE.COUNT = 0L

      if keyword_set(dostatusline) then begin
          statusline, /close
      endif

      RETURN
      END