#! /bin/sh

     ###################################### 
     ##### Model Transformation Tools #####
    ######################################

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$





## Revision 1.57  2001/04/01 03:38:54  geraint
## Reset row to zero after write to file, ready for subsequent runs.
## Eliminates SIGSEGV in Octave when _ode2odes called multiple times.
##
## Revision 1.56  2001/03/30 15:13:58  gawthrop
## Rationalised simulation modes to each return mtt_data
##

filename=${sys}_ode2odes.${lang}

if [ -n "$3" ]; then
  method=$3    
else
  method=implicit  
fi







echo Creating $filename with $method integration method

# Find system constants
Nx=`mtt_getsize $sys x` # States
Nu=`mtt_getsize $sys u` # Inputs 
Ny=`mtt_getsize $sys y` # Inputs  

	;;
esac

cat <<EOF  > $filename
#include <octave/oct.h>
#include <octave/load-save.h>
#include <octave/lo-mappers.h>
#include <octave/NLEqn.h>
#include <octave/ov-struct.h>
#include <octave/variables.h>

#ifndef STANDALONE
#include <octave/${feval_header}>
#endif

#include "${sys}_def.h"
#include "${sys}_sympar.h"

#ifdef STANDALONE
#include <csignal>
#include <fstream>



extern ColumnVector F${sys}_ae (
	ColumnVector &x,
	ColumnVector &u,
	const double &t,
	ColumnVector &par);

fi
cat <<EOF >> $filename

void set_signal_handlers (void);

#endif // STANDALONE

class Solver
{
public:
  Solver (void) { 
    solve_errors.input		= 0;
    solve_errors.user		= 0;
    solve_errors.converge	= 0;
    solve_errors.progress	= 0;
    solve_errors.limit		= 0;
    solve_errors.unknown	= 0;
 }
  ColumnVector solve (const ColumnVector &x,
		      const ColumnVector &u,
		      const double &t,
		      const ColumnVector &par)
  {
    initialise(x,u,t,par);
    NLFunc fcn(&Solver::evaluate);
    NLEqn  eqn(Ui,fcn);
    Ui = eqn.solve(info);
    switch (info)
    {
      case -2:
	solve_errors.input++;
	write_stats();
        break;
      case -1:
	solve_errors.user++;
	write_stats();
        break;
      case 1:
	solve_errors.converge++;
        break;
      case 3:
	solve_errors.progress++;
	write_stats();
        break;
      case 4:
	solve_errors.limit++;
	write_stats();
	break;
      default:
	solve_errors.unknown++;
	write_stats();
	break;;
    }
    U.insert (Ui,MTTNU);
    return U;
  }
  void write_stats (void)
  {
    cerr << "input error (" << solve_errors.input << ") "
	 << ",  user (" << solve_errors.user << ") "
	 << ",  converge (" << solve_errors.converge << ") "
	 << ",  progress (" << solve_errors.progress << ") "
	 << ",  limit (" << solve_errors.limit << ") "
	 << ",  unknown (" << solve_errors.unknown << ") "
	 << endl;
  }
private:
  static ColumnVector evaluate (const ColumnVector &tryUi)
  {
    U.insert(tryUi,MTTNU);
    return F${sys}_ae(X,U,T,P);
  }
  void initialise (const ColumnVector &x,
		   const ColumnVector &u,
		   const double &t,
		   const ColumnVector &par)
  {
    X = x;
    T = t;
    P = par;
    U.insert(u,0);
    Ui = ColumnVector(MTTNYZ,0.0);
  }
private:
  struct {
    long int input;
    long int user;
    long int converge;
    long int progress;
    long int limit;
    long int unknown; 
  } solve_errors;
  static double T;
  static ColumnVector X;
  static ColumnVector P;
  static ColumnVector U;
  static ColumnVector Ui;
  int info;
};
double Solver::T = 0.0;
ColumnVector Solver::U(MTTNU+MTTNYZ,0.0);
ColumnVector Solver::X(MTTNX,0.0);
ColumnVector Solver::P(MTTNPAR,0.0);
ColumnVector Solver::Ui(MTTNYZ,0.0);

inline ColumnVector
mtt_input (ColumnVector &x,
	   ColumnVector &y,
	   const double &t,
	   ColumnVector &par)
{
#ifdef STANDALONE
  static ColumnVector u  (MTTNU);
  static ColumnVector ui (MTTNYZ);
  static ColumnVector U  (MTTNU+MTTNYZ);
  static Solver ae;


  u = F${sys}_input (x, y, t, par);
  if (MTTNYZ == 0)
    {
      return u;
    }
  else