/*
Copyright (C) 2000 Marc Compere, compere@mail.utexas.edu

All Octave library functions and compiler macros were copied from
examples found in lsode.cc.
The one-line copyright notice for lsode.h is included below:

Copyright (C) 1996, 1997 John W. Eaton

The file sdirk4.h is intended for use with Octave.

This program is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2, or (at your option) any
later version.

This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
for more details.

You should have received a copy of the GNU General Public License
along with Octave; see the file COPYING.  If not, write to the Free
Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

*/

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <string>

#include <iostream.h>

#include "oct.h"

#include "ov-fcn.h"
#include "lo-mappers.h"

#include "defun-dld.h"
#include "error.h"
#include "gripes.h"
#include "oct-obj.h"
#include "ov-fcn.h"
#include "pager.h"
#include "unwind-prot.h"
#include "utils.h"
#include "variables.h"

#include <cfloat>
#include <cmath>

#include "f77-fcn.h"
#include "lo-error.h"



extern "C"
{
  int solout (int*, double*, double*,
               double*, int*, int*);

  int mas (int*, double*, int*);

  int jac (int*, double*, double*, double*, int*);

  int sdirk4_user_function (int*, double*, double*, double*);

  int F77_FCN (sdirk4, SDIRK4) (int*, int (*)(int*, double*, double*, double*),
                                 double*, double*, double*, double*,
                                 double*, double*, int*,
                                 int(*)(int*, double*, double*, double*, int*),
                                 int*, int*, int*,
                                 int(*)(int*, double*, int*),
                                 int*, int*, int*,
                                 int(*)(int*, double*, double*, double*, int*, int*),
                                 int*, double*, int*, int*,
                                 int*, int*, int*);

//      SUBROUTINE SDIRK4(N,FCN,
//                        X,Y,XEND,H,
//                        RTOL,ATOL,ITOL,
//                        JAC,
//                        IJAC,MLJAC,MUJAC,
//                        MAS,
//                        IMAS,MLMAS,MUMAS,
//                        SOLOUT,
//                        IOUT,WORK,LWORK,IWORK,
//                        LIWORK,LRCONT,IDID)
//
//                    SUBROUTINE FCN(N,X,Y,F)
//                    REAL*8 X,Y(N),F(N)
//
//                    SUBROUTINE JAC(N,X,Y,DFY,LDFY)
//                    REAL*8 X,Y(N),DFY(LDFY,N)
//
//                    SUBROUTINE MAS(N,AM,LMAS)
//                    REAL*8 AM(LMAS,N)
//
//                    SUBROUTINE SOLOUT (NR,XOLD,X,Y,N,IRTRN)
//                    REAL*8 X,Y(N)
}

int n;               // N (int, DIMENSION OF THE SYSTEM)
double h;            // INITIAL STEP SIZE GUESS, IF H=0.D0, THE CODE PUTS H=1.D-6
int itol;            // ITOL (int, ITOL=0: BOTH RTOL AND ATOL ARE SCALARS)
int ijac;            // IJAC (int, IJAC=0: JACOBIAN IS COMPUTED INTERNALLY)
int mljac;           // MLJAC (int, MLJAC=N: JACOBIAN IS A FULL MATRIX)
int mujac;           // MUJAC (int, NEED NOT BE DEFINED IF MLJAC=N.)
int imas;            // IMAS (int, IF IMAS=0, THIS MATRIX IS ASSUMED TO BE THE IDENTITY)
int mlmas;           // MLMAS (int, MLMAS=N: THE FULL MATRIX CASE)
int mumas;           // MUMAS (int, NEED NOT BE DEFINED IF MLMAS=N)
int iout;            // IOUT (int, IOUT=1: SUBROUTINE IS AVAILABLE FOR OUTPUT)
Array<double> work;  // WORK (double*, ARRAY OF WORKING SPACE OF LENGTH "LWORK".)
int lwork;           // LWORK (int, LWORK = 2*N*N+12*N+7)
Array<int> iwork;    // IWORK (int*, INTEGER WORKING SPACE OF LENGHT "LIWORK".)
int liwork;          // LIWORK (int, 2*N+4)
int lrcont;          // LRCONT (int, 5*N+2)
int idid;            // IDID (REPORTS ON SUCCESSFULNESS UPON RETURN,
                        //       IDID=1  COMPUTATION SUCCESSFUL,
                        //       IDID=-1 COMPUTATION UNSUCCESSFUL)


double t_start;
double t_stop;
double rel_tol;
double abs_tol;
double h_initial;


// Global pointer for user defined function required by sdirk4
static octave_function *sdirk4_fcn;

// Is this a recursive call?
static int call_depth = 0;

// global var's
octave_value_list args;
octave_value_list fout;
Matrix statehistory;
ColumnVector timehistory;
bool abscissa_trace;
int skip_step;
int step_counter;