# Copyright (C) 2001 by Peter J. Gawthrop

###############################################################
## Version control history
###############################################################
## $Header$
## $Log$
## Revision 1.326  2001/10/26 01:01:47  geraint
## fixcc when rdae_is_dae (-cr).
##
## Revision 1.325  2001/10/11 03:15:55  geraint
## Fixed make dependencies for c++ simulation code
## - re-simulation no longer required for different output format.
##
## Revision 1.324  2001/10/05 23:39:43  geraint
## Made odes.dat2 .PRECIOUS to prevent deletion when simulation receives ^C.
##

		    *)
			echo $1 is an unknown integration method - use dassl, euler, rk4 or implicit;
                        exit;;
		esac;;
	-ae )
                mtt_switches="$mtt_switches $1";
		case $2 in
		    dassl)
			mtt_switches="$mtt_switches $2";
			algebraic_solver=Dassl_Solver;
			shift;
			;;
		    fsolve | hybrd)
			mtt_switches="$mtt_switches $2";
			algebraic_solver=Hybrd_Solver;
			shift;
			;;
		    hooke)
			mtt_switches="$mtt_switches $2";

    echo '         -c  c-code generation'
    echo '         -cc C++ code generation'
    echo '         -cr Use cr before resolving equations'
    echo '         -d  <dir>  use directory <dir>'
    echo '         -dc Maximise derivative (not integral) causality'
    echo '         -dc Maximise derivative (not integral) causality'
    echo '         -i <implicit|euler|rk4|dassl>   Use implicit, euler, rk4 or dassl integration'
    echo '         -ae <reduce|hybrd|hooke>   Solve algebraic equations with Reduce, hybrd (fsolve) or Hooke and Jeeves'
    echo '         -ae <reduce|hybrd|dassl|hooke>   Solve algebraic equations with specified solver'
    echo '         -o ode is same as dae'
    echo '         -oct use oct files in place of m files where appropriate'
    echo '         -opt optimise code generation'
    echo '         -p  print environment variables'
    echo '         -partition partition hierachical system'
    echo '         -pdf generate pdf in place of ps'
    echo '         -r  reset time stamp on representation'

#include "mtt_Dassl_Solver.hh"


// used by "-ae dassl" NOT "-i dassl"


MTT::Dassl_Solver *MTT::Dassl_Solver::static_ptr;

ColumnVector
MTT::Dassl_Solver::f_dassl (const ColumnVector &tryUi,
			    const ColumnVector &tryUidot,
			    double t, int &ires)
{
  static MTT::Dassl_Solver *p = MTT::Dassl_Solver::static_ptr;

  ColumnVector residual = p->eval(tryUi);
  ColumnVector uidoterr = (tryUi - p->_ui) - tryUidot;  
  p->_yz = residual + uidoterr;
  return p->_yz;
}

void
MTT::Dassl_Solver::Solve (void)
{
  const double t0 = 0.0;
  const double t1 = 1.0;
  DAEFunc fcn(&MTT::Dassl_Solver::f_dassl);
  DASSL   eqn(AlgebraicSolver::_ui,t0,fcn);
  AlgebraicSolver::_ui = eqn.do_integrate (t1);
}

#ifndef MTT_DASSLSOLVER
#define MTT_DASSLSOLVER


#include "mtt_AlgebraicSolver.hh"
#include <octave/DASSL.h>


#ifdef  OCTAVE_DEV
#include <octave/parse.h>
#define VECTOR_VALUE column_vector_value
#else   // !OCTAVE_DEV
#include <octave/toplev.h>
#define VECTOR_VALUE vector_value
#endif  // OCTAVE_DEV


// -ae dassl


namespace MTT
{
  class Dassl_Solver : public MTT::AlgebraicSolver {
    
    // used only when called because of "-ae dassl"
    // this is not used when called by "-i dassl"
    
  public:
    
    Dassl_Solver (const int npar,
		  const int nu,
		  const int nx,
		  const int ny,
		  const int nyz)
      : MTT::AlgebraicSolver (npar, nu, nx, ny, nyz)
    {
      static_ptr = this;
    }
    
    static ColumnVector
    f_dassl (const ColumnVector &tryUi,
	     const ColumnVector &tryUidot,
	     double t, int &ires);
    
    ~Dassl_Solver (void) {};
    
  protected:
    
    void
    Solve (void);
    
  public:
    
    static Dassl_Solver *static_ptr;
    
  };
}


#endif // MTT_DASSLSOLVER