#! /bin/sh

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

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.46  2001/02/05 08:32:31  geraint
## typo
##
## Revision 1.45  2001/02/05 04:32:35  geraint
## Octave version 2.1.x compatability and #ifdef statements for standalone rep
##
## Revision 1.46  2001/01/08 06:21:59  geraint
## #ifdef STANDALONE stuff
##
## Revision 1.45  2001/01/07 01:25:49  geraint

cat <<EOF  > $filename
#include <octave/oct.h>

#include <octave/${feval_header}>
#include <octave/ov-struct.h>
#include <octave/oct-map.h>
#include <octave/lo-mappers.h>
#include <octave/variables.h>

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

#ifdef STANDALONE
#define DECLARE(name) extern octave_value_list F##name (const octave_value_list &, int);
DECLARE(mtt_euler)
DECLARE(mtt_implicit)
DECLARE(mtt_write)
DECLARE(${sys}_${ode})
DECLARE(${sys}_${odeo})
DECLARE(${sys}_input)
DECLARE(${sys}_numpar)
DECLARE(${sys}_simpar)
DECLARE(${sys}_smxa)
DECLARE(${sys}_smxax)
DECLARE(${sys}_state)
DECLARE(${sys}_logic)
#endif // STANDALONE

ColumnVector
mtt_${ode} (ColumnVector x, ColumnVector u, double t, ColumnVector par)
inline ColumnVector
mtt_${ode} (const ColumnVector &x,
	    const ColumnVector &u,
	    const double &t,
	    const ColumnVector &par)
{
  octave_value_list args, f;
  args (0) = octave_value (x);
  args (1) = octave_value (u);
  args (2) = octave_value (t);
  args (3) = octave_value (par);
#ifdef STANDALONE
  f = F${sys}_${ode} (args, 1);
#else
  f = feval ("${sys}_${ode}", args, 1);
#endif
  return f(0).${vector_value} ();
}

ColumnVector
mtt_${odeo} (ColumnVector x, ColumnVector u, double t, ColumnVector par)
inline ColumnVector
mtt_${odeo} (const ColumnVector &x,
	     const ColumnVector &u,
	     const double &t,
	     const ColumnVector &par)
{
  octave_value_list args, f;
  args (0) = octave_value (x);
  args (1) = octave_value (u);
  args (2) = octave_value (t);
  args (3) = octave_value (par);
#ifdef STANDALONE
  f = F${sys}_${odeo} (args, 1);
#else
  f = feval ("${sys}_${odeo}", args, 1);
#endif
  return f(0).${vector_value} ();
}

// #define mtt_implicit(x,dx,AA,AAx,ddt,nx,open) call_mtt_implicit((x),(dx),(AA),(AAx),(ddt),(nx),(open))
ColumnVector
mtt_implicit (ColumnVector x,
		   ColumnVector dx,
		   Matrix AA,
		   ColumnVector AAx,
		   double ddt,
		   int nx,
		   ColumnVector open_switches)
inline ColumnVector
mtt_implicit (const ColumnVector &x,
	      const ColumnVector &dx,
	      const Matrix &AA,
	      const ColumnVector &AAx,
	      const double &ddt,
	      const int &nx,
	      const ColumnVector &open_switches)
{
  octave_value_list args, f;
  args (0) = octave_value (x);
  args (1) = octave_value (dx);
  args (2) = octave_value (AA);
  args (3) = octave_value (AAx);
  args (4) = octave_value (ddt);
  args (5) = octave_value ((double)nx);
  args (6) = octave_value (open_switches);
#ifdef STANDALONE
  f = Fmtt_implicit (args, 1);
#else
  f = feval ("mtt_implicit", args, 1);
#endif
  return f(0).${vector_value} ();
}

ColumnVector
mtt_euler (ColumnVector x,
		   ColumnVector dx,
		   double ddt,
		   int nx,
		   ColumnVector open_switches)
inline ColumnVector
mtt_euler (const ColumnVector &x,
	   const ColumnVector &dx,
	   const double &ddt,
	   const int &nx,
	   const ColumnVector &open_switches)
{
  octave_value_list args, f;
  args (0) = octave_value (x);
  args (1) = octave_value (dx);
  args (2) = octave_value (ddt);
  args (3) = octave_value ((double)nx);
  args (4) = octave_value (open_switches);
#ifdef STANDALONE
  f = Fmtt_euler (args, 1);
#else
  f = feval ("mtt_euler", args, 1);
#endif
  return f(0).${vector_value} ();
}

ColumnVector
mtt_input (ColumnVector x, ColumnVector y, const double t, ColumnVector par)
inline ColumnVector
mtt_input (const ColumnVector &x,
	   const ColumnVector &y,
	   const double &t,
	   const ColumnVector &par)
{
  octave_value_list args;
  args (0) = octave_value (x);
  args (1) = octave_value (y);
  args (2) = octave_value (t);
  args (3) = octave_value (par);
  octave_value_list f;
#ifdef STANDALONE
  f = F${sys}_input (args, 1);
#else
  f = feval ("${sys}_input", args, 1);
#endif
  return f(0).${vector_value} ();
}

ColumnVector
inline ColumnVector
mtt_numpar (void)
{
  octave_value_list args;
  octave_value_list f;
#ifdef STANDALONE
  f = F${sys}_numpar (args, 1);
#else
  f = feval ("${sys}_numpar", args, 1);
#endif
  return f(0).${vector_value} ();
}

Octave_map
inline Octave_map
mtt_simpar (void)
{
  octave_value_list args;
  Octave_map f;
#ifdef STANDALONE
  f["first"]		= F${sys}_simpar (args, 1)(0).map_value ()["first"];
  f["dt"]		= F${sys}_simpar (args, 1)(0).map_value ()["dt"];

  f["wmax"]		= feval ("${sys}_simpar", args, 1)(0).map_value ()["wmax"];
  f["wsteps"]		= feval ("${sys}_simpar", args, 1)(0).map_value ()["wsteps"];
  f["input"]		= feval ("${sys}_simpar", args, 1)(0).map_value ()["input"];
#endif
  return (f);
}

Matrix
mtt_smxa (ColumnVector x, ColumnVector u, double t, ColumnVector par)
inline Matrix
mtt_smxa (const ColumnVector &x,
	  const ColumnVector &u,
	  const double &t,
	  const ColumnVector &par)
{
  octave_value_list args;
  args (0) = octave_value (x);
  args (1) = octave_value (u);
  args (2) = octave_value (t);
  args (3) = octave_value (par);
  octave_value_list f;
#ifdef STANDALONE
  f = F${sys}_smxa (args, 1);
#else
  f = feval ("${sys}_smxa", args, 1);
#endif
  return f(0).matrix_value ();
}

ColumnVector
mtt_smxax (ColumnVector x, ColumnVector u, double t, ColumnVector par)
inline ColumnVector
mtt_smxax (const ColumnVector &x,
	   const ColumnVector &u,
	   const double &t,
	   const ColumnVector &par)
{
  octave_value_list args;
  args (0) = octave_value (x);
  args (1) = octave_value (u);
  args (2) = octave_value (t);
  args (3) = octave_value (par);
  octave_value_list f;
#ifdef STANDALONE
  f = F${sys}_smxax (args, 1);
#else
  f = feval ("${sys}_smxax", args, 1);
#endif
  return f(0).${vector_value} ();
}

ColumnVector
mtt_state (ColumnVector x)
inline ColumnVector
mtt_state (const ColumnVector &x)
{
  octave_value_list args;
  args (0) = octave_value (x);
  octave_value_list f;
#ifdef STANDALONE
  f = F${sys}_state (args, 1);
#else
  f = feval ("${sys}_state", args, 1);
#endif
  return f(0).${vector_value} ();
}

ColumnVector
mtt_logic (ColumnVector x, ColumnVector u, double t, ColumnVector par)
inline ColumnVector
mtt_logic (const ColumnVector &x,
	   const ColumnVector &u,
	   const double &t,
	   const ColumnVector &par)
{
  octave_value_list args;
  args (0) = octave_value (x);
  args (1) = octave_value (u);
  args (2) = octave_value (t);
  args (3) = octave_value (par);

  octave_value_list f;
#ifdef STANDALONE
  f = F${sys}_logic (args, 1);
#else
  f = feval ("${sys}_logic", args, 1);
#endif
  return f(0).${vector_value} ();
}


void
mtt_write (double t, ColumnVector x, ColumnVector y){
  octave_value_list args;
  args (0) = octave_value (t);
  args (1) = octave_value (x);
  args (2) = octave_value (y);
inline void
mtt_write (const double &t,
	   const ColumnVector &x,
	   const ColumnVector &y,
	   const int &nrows)
{
  static Matrix data;
  static int row;
  const int nx = x.length (), ny = y.length ();
  register int col = 0;

  if (0 == row)
    data = Matrix (nrows, 1+ny+1+nx, 0.0);

  data (row, col) = t;
  for (register int i = 0; i < ny; i++)
    data (row, ++col) = y (i);
  data (row, ++col) = t;
  for (register int i = 0; i < nx; i++)
    data (row, ++col) = x (i);

  row++;

  if (nrows == row)
#ifdef STANDALONE
  Fmtt_write (args, 1);
#else
    cout << data;
#else // ! STANDALONE
  feval ("mtt_write", args, 1);
  set_global_value ("MTT_data", data);
#endif
}


//void
//mtt_write (double t, ColumnVector x, ColumnVector y, int nx, int ny)
//{
//  register int i;
//  cout.precision (5);		// this should be passed in as an argument
//  cout.width (12);		// as should this (instead of nx, ny)
//  cout << t;
//  for (i = 0; i < y.length (); i++)
//    {
//      cout.width (12);
//      cout << '\t' << y (i);
//    }
//  cout.width (12);
//  cout << "\t\t" << t;
//  for (i = 0; i < x.length (); i++)
//    {
//      cout.width (12);
//      cout << '\t' << x (i);
//    }
//  cout << endl;
//}

ColumnVector nozeros (const ColumnVector v0, const double tol = 0.0)
{
  ColumnVector v (v0.length ());
  register int j;
  for (register int i = j = 0; i < v.length (); i++)
    {
      if (tol <= abs(v0 (i)))
	{
	  v (j) = v0 (i);
	  j++;
	}
    }
  return (j)
    ? v.extract (0, --j)
    : ColumnVector ();
}


DEFUN_DLD (${sys}_ode2odes, args, ,
"Octave ode2odes representation of system 
Usage: ${sys}_ode2odes (x, par, simpar)
")
{
  octave_value_list retval;

  ColumnVector	x;
  ColumnVector	par;
  Octave_map	simpar;

  double
    first	= 0.0,
    dt		= 0.0,
    last	= 0.0,
    stepfactor	= 0.0;

  int nargin = args.length ();
  switch (nargin)
    {
    case 3:
      simpar["first"]		= args (2).map_value ()["first"];
      simpar["dt"]		= args (2).map_value ()["dt"];
      simpar["last"]		= args (2).map_value ()["last"];
      simpar["stepfactor"]     	= args (2).map_value ()["stepfactor"];
      first		= args (2).map_value ()["first"].double_value ();
      dt		= args (2).map_value ()["dt"].double_value ();
      last		= args (2).map_value ()["last"].double_value ();
      stepfactor     	= args (2).map_value ()["stepfactor"].double_value ();
      simpar["wmin"]		= args (2).map_value ()["wmin"];
      simpar["wmax"]		= args (2).map_value ()["wmax"];
      simpar["wsteps"]		= args (2).map_value ()["wsteps"];
      simpar["input"]		= args (2).map_value ()["input"];
      par    = args (1).${vector_value} ();
      x      = args (0).${vector_value} ();
      par    		= args (1).${vector_value} ();
      x      		= args (0).${vector_value} ();
      break;
    case 2:
      simpar["first"]		= mtt_simpar ()["first"];
      simpar["dt"]		= mtt_simpar ()["dt"];
      simpar["last"]		= mtt_simpar ()["last"];
      simpar["stepfactor"]     	= mtt_simpar ()["stepfactor"];
      first		= mtt_simpar ()["first"].double_value ();
      dt		= mtt_simpar ()["dt"].double_value ();
      last		= mtt_simpar ()["last"].double_value ();
      stepfactor     	= mtt_simpar ()["stepfactor"].double_value ();
      simpar["wmin"]		= mtt_simpar ()["wmin"];
      simpar["wmax"]		= mtt_simpar ()["wmax"];
      simpar["wsteps"]		= mtt_simpar ()["wsteps"];
      simpar["input"]		= mtt_simpar ()["input"];
      par    = args (1).${vector_value} ();
      x      = args (0).${vector_value} ();
      par    		= args (1).${vector_value} ();
      x      		= args (0).${vector_value} ();
      break;
    case 1:
      simpar["first"]		= mtt_simpar ()["first"];
      simpar["dt"]		= mtt_simpar ()["dt"];
      simpar["last"]		= mtt_simpar ()["last"];
      simpar["stepfactor"]     	= mtt_simpar ()["stepfactor"];
      first		= mtt_simpar ()["first"].double_value ();
      dt		= mtt_simpar ()["dt"].double_value ();
      last		= mtt_simpar ()["last"].double_value ();
      stepfactor     	= mtt_simpar ()["stepfactor"].double_value ();
      simpar["wmin"]		= mtt_simpar ()["wmin"];
      simpar["wmax"]		= mtt_simpar ()["wmax"];
      simpar["wsteps"]		= mtt_simpar ()["wsteps"];
      simpar["input"]		= mtt_simpar ()["input"];
      par    = mtt_numpar ();
      x      = args (0).${vector_value} ();
      par    		= mtt_numpar ();
      x      		= args (0).${vector_value} ();
      break;
    case 0:
      simpar["first"]		= mtt_simpar ()["first"];
      simpar["dt"]		= mtt_simpar ()["dt"];
      simpar["last"]		= mtt_simpar ()["last"];
      simpar["stepfactor"]     	= mtt_simpar ()["stepfactor"];
      first		= mtt_simpar ()["first"].double_value ();
      dt		= mtt_simpar ()["dt"].double_value ();
      last		= mtt_simpar ()["last"].double_value ();
      stepfactor     	= mtt_simpar ()["stepfactor"].double_value ();
      simpar["wmin"]		= mtt_simpar ()["wmin"];
      simpar["wmax"]		= mtt_simpar ()["wmax"];
      simpar["wsteps"]		= mtt_simpar ()["wsteps"];
      simpar["input"]		= mtt_simpar ()["input"];
      par    = mtt_numpar ();
      x      = mtt_state (par);
      par    		= mtt_numpar ();
      x      		= mtt_state (par);
      break;
    default:
      usage("${sys}_ode2odes (x par simpar)", nargin);
      error("aborting.");
    }

  ColumnVector	dx (MTTNX);
  ColumnVector	u (MTTNU);
  ColumnVector	y (MTTNY);

  Matrix	AA (MTTNX, MTTNX);
  ColumnVector	AAx (MTTNX);

  ColumnVector	open_switches (MTTNX);

  register double t	= 0.0;

  const double	ddt	= simpar ["dt"].double_value () / simpar ["stepfactor"].double_value ();
  const int	ilast	= static_cast<int> (round (simpar ["last"].double_value () / ddt)) + 1;

  const double	ddt	= dt / stepfactor;
  const int	ilast	= static_cast<int> (round ((last - first) / ddt)) + 1;
  const int	nrows	= static_cast<int> (round ((last - first) / dt)) + 1;
  // cse translation
  // LSODE will need ODEFUNC

  for (register int j = 0, i = 1; i <= ilast; i++)
    {
      y	= mtt_${odeo} (x, u, t, par);
      u	= mtt_input (x, y, t, par);
      if (0 == j)
	{
	  //mtt_write (t, x, y, MTTNX, MTTNY);
           mtt_write (t, x, y);
           mtt_write (t, x, y, nrows);
	}
      dx = mtt_${ode} (x, u, t, par);
EOF

if [ "$method" = "implicit" ]; then
cat <<EOF >> $filename

      AA = mtt_smxa (x, u, ddt, par);
      AAx = mtt_smxax (x, u, ddt, par);
EOF
fi

## Common stuff
cat <<EOF >> $filename
      open_switches = mtt_logic (x, u, t, par);
      x =  $algorithm; 
      t += ddt;
      j++;
      j = (j == static_cast<int> (simpar ["stepfactor"].double_value ())) ? 0 : j;
      j = (j == static_cast<int> (stepfactor)) ? 0 : j;
    }

  retval (0) = octave_value (y);
  retval (1) = octave_value (x);
  retval (2) = octave_value (t);
  return (retval);
}

#include <octave/oct.h>
#include <octave/variables.h>

#ifdef STANDALONE
static Matrix MTT_data;
#endif

DEFUN_DLD (mtt_write, args, ,
	   "append current data to output")
{
  const double		t	= args(0).double_value ();
#ifdef OCTAVE_DEV
  const ColumnVector	x	= args(1).column_vector_value ();
  const ColumnVector	y	= args(2).column_vector_value ();
#else
  const ColumnVector	x	= args(1).vvector_value ();
  const ColumnVector	y	= args(2).vector_value ();
#endif

  const int		nx	= x.length ();
  const int		ny	= y.length ();


  ColumnVector Output (2+nx+ny, 0.0);
  
  Output (0) = Output (1+nx) = t;
  Output.insert (x.transpose (), 1);
  Output.insert (y.transpose (), 2+nx);

  Matrix data;

  if (0.0 == t)
    {
      data = static_cast<Matrix> (Output.transpose ());
    }
  else
    {
#ifdef STANDALONE
      data = MTT_data.transpose ();
#else
      data = get_global_value ("MTT_data").matrix_value ().transpose ();
#endif
      data = data.append (Output);
    }
  data = data.transpose ();
#ifdef STANDALONE
  MTT_data = data;
  cout << Output.transpose () << endl;
#else
  set_global_value ("MTT_data", data);
#endif
  return data;
}