#! /bin/sh

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

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.70  2002/04/28 18:41:27  geraint
## Fixed [ 549658 ] awk should be gawk.
## Replaced calls to awk with call to gawk.
##
## Revision 1.69  2002/04/17 13:46:58  geraint
## #include <fstream> for -oct as well as -cc.
##
## Revision 1.68  2002/04/15 10:54:31  geraint
## Statically declare outputs and initialise to zero.
## This is necessary to prevent spurious values from being output when no assignments are made (i.e. when "y(i) := 0 for all u"  (Reduce:see NERO)).
##
## Revision 1.67  2002/04/09 12:04:21  geraint
## Replaced ios:: with std::ios:: for g++-3.0 compatability.
##
## Revision 1.66  2002/03/26 11:58:58  geraint
## Added cputime monitoring.
##
## Revision 1.65  2001/11/15 06:24:11  geraint
## Updated (-i dassl) residual function to use new DAEFunc (octave-2.1.35).
## YZ residual dependency on Ui still requires some work.
##
## Revision 1.64  2001/08/08 02:15:00  geraint
## Rationalisation of solver code, beginning with algebraic solvers.
##
## Revision 1.63  2001/08/07 04:39:24  geraint
## Consolidated dassl and residual functions.
##
## Revision 1.62  2001/08/01 22:14:32  geraint
## Bug fix for dassl.
##
## Revision 1.61  2001/08/01 04:06:07  geraint
## Added -i dassl for -cc and -oct.
##
## Revision 1.60  2001/07/16 22:23:00  geraint
## Fixed misleading variable name in .cc rep.
##
## Revision 1.59  2001/07/13 04:54:04  geraint
## Branch merge: numerical-algebraic-solution back to main.
##
## Revision 1.58  2001/07/13 00:51:39  geraint
## Fixed generation of odes.sg from .m and .oct simulations.
## .cc, .m and .oct simulations now all write mtt_data (lower case).
##
## Revision 1.57.2.5  2001/07/13 04:02:31  geraint
## Implemented numerical algebraic solution for _ode2odes.oct.
##
## Revision 1.57.2.4  2001/07/02 00:34:56  geraint
## gcc-3.0 compatibility.
##
## Revision 1.57.2.3  2001/06/25 23:28:29  geraint
## Generic mtt_rate and mtt_output - allows method independent calls.
##
## Revision 1.57.2.2  2001/06/05 03:20:40  geraint
## added -ae option to select algebraic equation solution method.
##
## Revision 1.57.2.1  2001/05/04 04:07:24  geraint
## Numerical solution of algebraic equations.
## sys_ae.cc written for unsolved inputs.
## Solution of equations using hybrd from MINPACK (as used by Octave fsolve).
##
## 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
##
## Revision 1.55  2001/03/27 13:21:59  geraint
## Octave version compatibility for save_ascii_data(_for_plotting).
##
## Revision 1.54  2001/03/27 01:14:27  geraint
## Improved determination of Octave version.
##
## Revision 1.53  2001/03/21 03:24:59  geraint
## Calculate inputs before outputs (.cc).
##
## Revision 1.52  2001/03/19 02:28:52  geraint
## Branch merge: merging-ode2odes-exe back to MAIN.
##
## Revision 1.51.2.7  2001/03/17 09:51:07  geraint
## Implemented Runge-Kutta IV fixed-step method (-i rk4).
##
## Revision 1.51.2.6  2001/03/16 03:56:13  geraint
## Removed psignal/siginfo.h - problematic and unnecessary.
##
## Revision 1.51.2.5  2001/03/12 23:16:37  geraint
## Minor improvements to signal handling (.exe).
##
## Revision 1.51.2.4  2001/03/12 03:59:30  geraint
## SIGINT  (C-c C-c) now causes simulation data to be dumped to MTT.core.
## SIGQUIT (C-c C-\) as for SIGINT, then raises default SIGQUIT.
## SIGFPE  as for SIGINT, then raises default SIGABRT.
##
## Revision 1.51.2.3  2001/03/07 04:06:55  geraint
## Irix: catch SIGFPE and write data before aborting (.exe).
## GNU/Linux: nada.
##
## Revision 1.51.2.2  2001/03/02 00:45:21  geraint
## Separated Euler and Implicit methods in .cc code and dependencies.
##
## Revision 1.51.2.1  2001/03/01 05:05:53  geraint
## Minor revisions.
##
## Revision 1.51  2001/02/19 06:33:19  geraint
## Removed operation form loop.
##
## Revision 1.50  2001/02/18 09:18:49  geraint
## Removed temporary Matrices from mtt_implicit.cc
##
## Revision 1.49  2001/02/14 06:06:34  geraint
## Removed octave_value_list wrappers from standalone.exe - speed improvements
##
## Revision 1.48  2001/02/11 07:08:59  geraint
## Static declarations of octave_value_lists: small .exe speed improvement
##
## Revision 1.47  2001/02/11 05:25:52  geraint
## Reduced number of matrix operations during .oct simulation data write
##
## 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
## Compatibility with Octave 2.1.33
##
## Revision 1.44  2000/12/05 12:11:45  peterg
## Changed function name to name()
##
## Revision 1.43  2000/12/04 10:59:40  peterg
## *** empty log message ***
##
## Revision 1.42  2000/11/10 14:19:50  peterg
## Corrected the csex and cseo functions
##
## Revision 1.41  2000/11/09 17:06:39  peterg
## Now does euler for cc
##
## Revision 1.40  2000/10/17 09:55:00  peterg
## Replaced switchopen by logic
##
## Revision 1.39  2000/10/14 08:04:40  peterg
## Changed arguments to _inout for consistency
##
## Revision 1.38  2000/10/11 09:08:08  peterg
## cse --> csex
##
## Revision 1.37  2000/08/01 12:25:06  peterg
## Now includes euler
##
## Revision 1.36  2000/05/19 17:48:16  peterg
## Argument to state
##
## Revision 1.35  2000/05/18 18:59:40  peterg
## Removed the First time stuff
##
## Revision 1.34  2000/05/16 18:56:14  peterg
## *** empty log message ***
##
## Revision 1.33  2000/05/11 19:33:18  peterg
## Uniform version for _sim.m
##
## Revision 1.32  2000/05/11 08:30:00  peterg

##
## Revision 1.31  2000/05/10 18:33:25  peterg
## Use smxa and smxax in place of smx
##
## Revision 1.30  2000/04/18 11:24:19  peterg
## Removed _numpar.
##
## Revision 1.29  2000/04/07 19:10:57  peterg
## *** empty log message ***
##
## Revision 1.28  1999/12/08 05:56:52  peterg
## Reordered the writing of the input and output.
## Note that last value now discarded.
##
## Revision 1.27  1999/11/15 22:47:53  peterg
## Generates method-specific code.
##
## Revision 1.26  1999/10/20 01:31:43  peterg
## *** empty log message ***
##
## Revision 1.25  1999/08/29 06:55:26  peterg
## Removed [MTTu]  = zero_input($Nu);	# Zero the input
## to avoide the p2c bug ????
##
## Revision 1.24  1999/08/27 06:02:16  peterg
## removed zero_input to avoid p2c bug
##
## Revision 1.23  1999/08/02 13:39:19  peterg
## Replaced zero_vector by zero_input
##
## Revision 1.22  1999/04/20 06:16:07  peterg
## Removed initialisation of AA and AAx
## Remove _switch calls -- uses _switchopen exclusively
##
## Revision 1.21  1999/04/02 06:29:25  peterg
## New implicit method - solves numerical prob with ISW
##
## Revision 1.20  1999/04/02 02:13:58  peterg
## Back to RCS
##
## Revision 1.19  1999/03/30 21:39:25  peterg
## In implicit approach, set derivatives to zero (when switch is off)
## before update. This seems to stop numerical leakage though non-return
## switches.
##
## Revision 1.18  1999/03/15 01:17:07  peterg
## Removed some spurious debugging code
##
## Revision 1.17  1999/03/15 01:09:15  peterg
## Fixed bugs when Nx=0 (no state)
##
## Revision 1.16  1999/03/06 02:28:38  peterg
## Rearranged evaluation to: state - input - output - write
##
## Revision 1.15  1999/03/06 02:19:43  peterg
## Changed args to _input
##
## Revision 1.14  1998/10/01 16:02:01  peterg
## Integration with switches handled separately fro Euler and Implicit.
##
## Revision 1.13  1998/09/30 17:41:24  peterg
## Implicit method now allows for switches via _switchA
##
## Revision 1.12  1998/08/27 08:55:18  peterg
## Mods to integration methods
##
## Revision 1.11  1998/08/25 12:28:31  peterg
## Move initila switch to after initial input
##
## Revision 1.10  1998/08/25 12:22:45  peterg
## Put _switch after update and also at initilisation
##
## Revision 1.9  1998/08/15 13:46:59  peterg
## New versions of integration routines
##
## Revision 1.8  1998/08/11 13:28:03  peterg
## Lowercase mttLAST etc
##
## Revision 1.7  1998/07/30 11:29:54  peterg
## Added implicit integration stuff
##
## Revision 1.6  1998/07/30 10:44:37  peterg
## INcluded othe integration methods.
##
## Revision 1.5  1998/07/26 11:02:20  peterg
## Put mtt or MTT  in front of variable names to avoid clashes with
## globals
##
## Revision 1.4  1998/07/25 20:14:00  peterg
## update code added for flexibility and octave efficiency
##
###############################################################


# Bourne shell script: make_ode2odes

# Copyright (c) P.J.Gawthrop July 1998.

# Tell user
sys=$1
lang=$2
filename=${sys}_ode2odes.${lang}

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

if [ -n "$4" ]; then
    algebraic_solver=$4
else
  algebraic_solver="Reduce_Solver"
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` # Outputs  

case "$method" in
    "implicit")
	ode=csex
	odeo=cseo
	algorithm="mtt_implicit(x,dx,AA,AAx,ddt,$Nx,open_switches)"
	;;
    "dassl")
	ode=ode
	odeo=odeo
	algorithm="mtt_dassl(x,u,t,par,dx,ddt,MTTNX,MTTNYZ,open_switches)"
	;;
    "euler" | "rk4" | *)
	ode=ode
	odeo=odeo
	algorithm="mtt_euler(x,dx,ddt,$Nx,open_switches)"
	;;
esac

make_m()
{
#lang_header  $1 ode2odes m  'x,par,simpar' '[Y,X,t]' > $filename
mtt_header ${sys} ode2odes m > $filename
cat <<EOF >> $filename 
global MTT_data;

  if nargin<3
    simpar = ${sys}_simpar(); 
    [simpar.dt] = mtt_simpar_update;
  endif
  if nargin<2
    par = ${sys}_numpar(); 
    [par] = mtt_numpar_update(par);
  endif
  if nargin<1
    [x0]  =  ${sys}_state(par); 
    [x0] = mtt_state_update(x);
  endif

  ## Initialise
  t = 0.0;
  ddt = simpar.dt/simpar.stepfactor;
  ilast = round(simpar.last/ddt)+1; # Total number of steps
  x = x0;

## Following removed due to p2c bug
##  [u]  = zero_input($Nu);		# Zero the input
  for MTTi=1:$Ny
    y(MTTi) = 0;
  endfor;

  mttj = 0;
    for it = 1:ilast		#Integration loop
    [u] = ${sys}_input(x,y,t,par);	# Input
    [y] = ${sys}_$odeo(x,u,t,par);  # Output 
    if mttj==0
     mtt_write(t,x,y,$Nx,$Ny);   # Write it out
    endif
EOF

if [ "$method" = "rk4" ]; then
cat << EOF >> $filename
    [k1] = ddt * ${sys}_${ode}(x,u,t,par);
    [k2] = ddt * ${sys}_${ode}(x+k1/2,u,t+ddt/2,par);
    [k3] = ddt * ${sys}_${ode}(x+k2/2,u,t+ddt/2,par);
    [k4] = ddt * ${sys}_${ode}(x+k3,u,t+ddt,par);
    [dx] = [k1 + 2.0 * [k2 + k3] + k4] / (6.0 * ddt);
EOF
else
cat << EOF >> $filename
    [dx] = ${sys}_$ode(x,u,t,par);	# State derivative
EOF
fi

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

    [AA] = ${sys}_smxa(x,u,ddt,par);	# (I-Adt) and (I-Adt)x
    [AAx] = ${sys}_smxax(x,u,ddt,par); # (I-Adt) and (I-Adt)x

EOF
fi

cat <<EOF >> $filename
    [open_switches] = ${sys}_logic(x,u,t,par); # Switch logic
    [x] = $algorithm; # Integration update
    t = t + ddt;   # Time update
    mttj = mttj+1;    # Increment counter
    if mttj==simpar.stepfactor
       mttj = 0;      # Reset counter
    endif

  endfor;			# Integration loop

## Create the output data
    mtt_data = MTT_data;
endfunction

EOF
} # make_m

make_cc() 
{
# get octave version
case `$MATRIX --version | gawk -F\. '{print $2}'` in
    0) # stable
	vector_value=vector_value
	feval_header=toplev.h
	save_ascii_data_function=save_ascii_data
	;;
    1) # development
	vector_value=column_vector_value
	feval_header=parse.h
	save_ascii_data_function=save_ascii_data_for_plotting
	;;
    *)
	vector_value=column_vector_value
	feval_header=parse.h
	save_ascii_data_function=save_ascii_data_for_plotting
	;;
esac

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

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

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

#include "mtt_${algebraic_solver}.hh"

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

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

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

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

extern ColumnVector F${sys}_numpar (
	void);

extern ColumnVector F${sys}_simpar (
     	void);

extern ColumnVector F${sys}_state (
	ColumnVector &par);

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

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

EOF
case "$method" in
    "implicit")
	cat <<EOF >> $filename
extern ColumnVector Fmtt_implicit (
	ColumnVector &x,
	ColumnVector &dx,
	Matrix &AA,
	ColumnVector &AAx,
	const double &ddt,
	const int &nx,
	const ColumnVector &open_switches);

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

extern ColumnVector F${sys}_smxax (
	ColumnVector &x,
	ColumnVector &u,
	const double &t,
	ColumnVector &par);
 
EOF
    ;;
    "dassl")
	cat <<EOF >> $filename
extern ColumnVector Fmtt_dassl (
	ColumnVector &x,
	const ColumnVector &u,
	const double &t,
	const ColumnVector &par,
	const ColumnVector &dx,
	const double &ddt,
	const int nx,
	const int nyz,
	const ColumnVector &open_switches);

EOF
    ;;
    "euler" | "rk4" | *)
	cat <<EOF >> $filename
extern ColumnVector Fmtt_euler (
	ColumnVector &x,
	const ColumnVector &dx,
	const double &ddt,
	const int &nx,
	const ColumnVector &open_switches);
 
EOF
    ;;
esac
cat <<EOF >> $filename

void set_signal_handlers (void);

#endif // STANDALONE

ColumnVector
mtt_ae (ColumnVector &x,
	ColumnVector &u,
	const double &t,
	ColumnVector &par)
{
#ifdef STANDALONE
    return F${sys}_ae(x,u,t,par);
#else
    static 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);
    f = feval ("${sys}_ae", args, 1);
    return f(0).${vector_value} ();
#endif
}

inline ColumnVector
mtt_input (ColumnVector &x,
	   ColumnVector &y,
	   const double &t,
	   ColumnVector &par)
{
  static MTT::${algebraic_solver} ae(MTTNPAR,MTTNU,MTTNX,MTTNY,MTTNYZ);
  static ColumnVector u  (MTTNU);

#ifdef STANDALONE
  u = F${sys}_input (x, y, t, par);
#else
  static octave_value_list args, f;
  args (0) = octave_value (x);
  args (1) = octave_value (y);
  args (2) = octave_value (t);
  args (3) = octave_value (par);
  f = feval ("${sys}_input", args, 1);
  u = f(0).${vector_value} ();
#endif
  if (MTTNYZ == 0)
    {
      return u;
    }
  else
    {
      return ae.solve(x,u,t,par);
    }
}

inline ColumnVector
mtt_logic (ColumnVector &x,
	   ColumnVector &u,
	   const double &t,
	   ColumnVector &par)
{
#ifdef STANDALONE
  return F${sys}_logic (x, u, t, par);
#else
  static 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);
  f = feval ("${sys}_logic", args, 1);
  return f(0).${vector_value} ();
#endif
}

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

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

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

inline ColumnVector
mtt_rate (ColumnVector &x,
	  ColumnVector &u,
	  const double &t,
	  ColumnVector &par)
{
#ifdef STANDALONE
  return F${sys}_${ode} (x, u, t, par);
#else
  static 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);
  f = feval ("${sys}_${ode}", args, 1);
  return f(0).${vector_value} ();
#endif
}

inline ColumnVector
mtt_output (ColumnVector &x,
	     ColumnVector &u,
	     const double &t,
	     ColumnVector &par)
{
#ifdef STANDALONE
  return F${sys}_${odeo} (x, u, t, par);
#else
  static 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);
  f = feval ("${sys}_${odeo}", args, 1);
  return f(0).${vector_value} ();
#endif
}

EOF
case "$method" in
    "implicit")
	cat <<EOF >> $filename
inline ColumnVector
mtt_implicit (ColumnVector &x,
	      ColumnVector &dx,
	      Matrix &AA,
	      ColumnVector &AAx,
	      const double &ddt,
	      const int &nx,
	      const ColumnVector &open_switches)
{
#ifdef STANDALONE
  return Fmtt_implicit (x, dx, AA, AAx, ddt, nx, open_switches);
#else
  static 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);
  f = feval ("mtt_implicit", args, 1);
  return f(0).${vector_value} ();
#endif
}

inline Matrix
mtt_smxa (ColumnVector &x,
	  ColumnVector &u,
	  const double &t,
	  ColumnVector &par)
{
#ifdef STANDALONE
  return F${sys}_smxa (x, u, t, par);
#else
  static 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);
  f = feval ("${sys}_smxa", args, 1);
  return f(0).matrix_value ();
#endif
}

inline ColumnVector
mtt_smxax (ColumnVector &x,
	   ColumnVector &u,
	   const double &t,
	   ColumnVector &par)
{
#ifdef STANDALONE
  return F${sys}_smxax (x, u, t, par);
#else
  static 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);
  f = feval ("${sys}_smxax", args, 1);
  return f(0).${vector_value} ();
#endif
}

EOF
    ;;
    "dassl")
cat <<EOF >> $filename
inline ColumnVector
mtt_dassl (ColumnVector &x,
	   const ColumnVector &u,
	   const double &t,
	   const ColumnVector &par,
	   const ColumnVector &dx,
	   const double &ddt,
	   const int &nx,
	   const int &nyz,
	   const ColumnVector &open_switches)
{
#ifdef STANDALONE
   return Fmtt_dassl (x, u, t, par, dx, ddt, nx, nyz, open_switches);
#else
   static 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);
    args (4) = octave_value (dx);
    args (5) = octave_value (ddt);
    args (6) = octave_value (static_cast<double> (nx));
    args (7) = octave_value (static_cast<double> (nyz));
    args (8) = octave_value (open_switches);
    f = feval ("mtt_dassl", args, 1);
    return f(0).${vector_value} ();
#endif
}

#ifdef STANDALONE
ColumnVector
Fmtt_residual (const ColumnVector &X,
	      const ColumnVector &DX,
	      double t,
	      int &ires)
{
#else // !STANDALONE
DEFUN_DLD (mtt_residual, args, , "mtt_residual")
{
    static ColumnVector	X  (MTTNX+MTTNYZ);
    static ColumnVector	DX (MTTNX+MTTNYZ);
    static double      	t;
    static int		&ires;

    X  = args(0).${vector_value} ();
    DX = args(1).${vector_value} ();
    t  = args(2).double_value ();
    ires = static_cast<int>(args(3).double_value ());
#endif // STANDALONE

    static ColumnVector residual (MTTNX+MTTNYZ);
    static ColumnVector U (MTTNU+MTTNYZ);
    static ColumnVector u (MTTNU);
    static ColumnVector y (MTTNY,0.0);
    static ColumnVector par (MTTNPAR);
    static ColumnVector dx(MTTNX);
    static ColumnVector yz(MTTNYZ);

    static ColumnVector x (MTTNX);
    static ColumnVector ui (MTTNYZ);

    static octave_value_list new_args;

    x = X.extract (0,MTTNX-1);
    if (MTTNYZ > 0)
      ui = X.extract (MTTNX,MTTNX+MTTNYZ-1);    

#ifdef STANDALONE
    par = F${sys}_numpar();
    u = F${sys}_input(x,y,t,par);
#else
    par = feval ("${sys}_numpar",  new_args, 1)(0).${vector_value} ();
    new_args(0) = octave_value (x);
    new_args(1) = octave_value (u);
    new_args(2) = octave_value (t);
    new_args(3) = octave_value (par);
    u = feval ("${sys}_input", new_args, 1)(0).${vector_value} ();
#endif

    U.insert (u,0);
    if (MTTNYZ > 0)
      U.insert (ui,MTTNU);

#ifdef STANDALONE
    dx = F${sys}_${ode} (x,U,t,par);
    yz = F${sys}_ae (x,U,t,par);
#else
    new_args(1) = octave_value (U);
    dx = feval ("${sys}_${ode}", new_args, 1)(0).${vector_value} ();
    yz = feval ("${sys}_ae", new_args, 1)(0).${vector_value} ();
#endif

    for (register int i = 0; i < MTTNX; i++)
      residual (i) = dx(i) - DX(i);
    
    if (MTTNYZ > 0)
      {
	residual.insert (yz,MTTNX);

	// XXX:
	// DASSL needs residual to be dependent on Ui and Uidot
	// mtt_dassl always sets the initial Ui to zero, so
	// Ui - h*Uidot should equal zero BUT, we don't know h
	static double t_old;
	double step;
	if (t != t_old)
	  {
	    step = t - t_old;
	    t = t_old;
	  }
	else
	  step = t;	
	for (register int i = MTTNX; i < MTTNX+MTTNYZ; i++)
	  residual(i) += X(i) - DX(i)*step;
      }

#ifdef STANDALONE
    return residual;
#else // !STANDALONE
    return octave_value (residual);
#endif // STANDALONE
}
    
EOF
    ;;
    "euler" | "rk4" | *)
cat <<EOF >> $filename
inline ColumnVector
mtt_euler (ColumnVector &x,
	   const ColumnVector &dx,
	   const double &ddt,
	   const int &nx,
	   const ColumnVector &open_switches)
{
#ifdef STANDALONE
  return Fmtt_euler (x, dx, ddt, nx, open_switches);
#else
  static 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);
  f = feval ("mtt_euler", args, 1);
  return f(0).${vector_value} ();
#endif
}

EOF
    ;;
esac
cat <<EOF >> $filename

inline void
mtt_write (const double &t,
	   ColumnVector &x,
	   ColumnVector &y,
	   const int &nrows,
	   const bool dump_data = false,
	   std::ostream &file = std::cout)
{
  static Matrix data;
  static int row;

  if (dump_data)
    {
       if (row > 0)
         {
           Matrix written_data = data.extract (0, 0, row-1, data.cols ()-1);
           $save_ascii_data_function (file, written_data, "mtt_data");
	 }
       return;
    }

  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.elem (row, col) = t;
  for (register int i = 0; i < ny; i++)
    data.elem (row, ++col) = y.elem (i);
  data.elem (row, ++col) = t;
  for (register int i = 0; i < nx; i++)
    data.elem (row, ++col) = x.elem (i);

  row++;

  static std::fstream fcputime ("MTT.cputime", std::ios::out | std::ios::trunc | std::ios::app);
  static clock_t cputime0 = clock();
  static clock_t cputime1 = cputime0;
  clock_t cputime = clock();
  fcputime << t << '\t'
	   << static_cast <double> (cputime - cputime0) / CLOCKS_PER_SEC << '\t'
	   << static_cast <double> (cputime - cputime1) / CLOCKS_PER_SEC << std::endl;
  cputime1 = cputime;

  if (nrows == row)
    {
#ifdef STANDALONE
      $save_ascii_data_function (file, data, "mtt_data");
//    std::cout << data << std::endl;
#else // ! STANDALONE
  set_global_value ("MTT_data", data);
#endif
    row = 0;
    fcputime.close();
    }
}

#ifdef STANDALONE
void dump_data (std::ostream &file)
{
  ColumnVector null (0);
  mtt_write (0.0, null, null, 0, true, file);
}

void handle_signal (int signum)
{
  // handle some signals to ensure data is written.
  std::cerr << "# Writing data to MTT.core (signal " << signum << ")" << std::endl;
  std::ofstream corefile ("MTT.core");
  dump_data (corefile);
  switch (signum)
    {
      case SIGFPE:
        // Intel chips do not raise SIGFPE for DIVZERO :-(
//        raise (SIGABRT);
        break;
      case SIGINT:
	break;
      case SIGQUIT:
	signal (SIGQUIT, SIG_DFL);
	raise (SIGQUIT);
	break;
      default:
        std::cerr << "# Warning: make_ode2odes needs updating!" << std::endl;
	signal (signum, SIG_DFL);
	raise (signum);
        break;
    }
  corefile.close ();
  set_signal_handlers ();
}

void set_signal_handlers (void)
{
  signal (SIGFPE,	handle_signal);
  signal (SIGINT,	handle_signal);
  signal (SIGQUIT,	handle_signal);
}

int main (void) {
  set_signal_handlers ();
#else
DEFUN_DLD (${sys}_ode2odes, args, ,
"Octave ode2odes representation of system with $method integration method\nUsage: mtt_data = ${sys}_ode2odes (x0, par, simpar)\n")
{
  static octave_value_list retval;
#endif // STANDALONE
  static ColumnVector	x0;
  static ColumnVector	par;
  static ColumnVector	simpar;

  static double
    first	= 0.0,
    dt		= 0.0,
    last	= 0.0,
    stepfactor	= 0.0;
#ifndef STANDALONE
  int nargin = args.length ();
  switch (nargin)
    {
    case 3:
      first		= args (2).${vector_value} ()(0);
      dt		= args (2).${vector_value} ()(1);
      last		= args (2).${vector_value} ()(2);
      stepfactor     	= args (2).${vector_value} ()(3);
      par    		= args (1).${vector_value} ();
      x0      		= args (0).${vector_value} ();
      break;
    case 2:
      first		= mtt_simpar ()(0);
      dt		= mtt_simpar ()(1);
      last		= mtt_simpar ()(2);
      stepfactor     	= mtt_simpar ()(3);
      par    		= args (1).${vector_value} ();
      x0      		= args (0).${vector_value} ();
      break;
    case 1:
      first		= mtt_simpar ()(0);
      dt		= mtt_simpar ()(1);
      last		= mtt_simpar ()(2);
      stepfactor     	= mtt_simpar ()(3);
      par    		= mtt_numpar ();
      x0      		= args (0).${vector_value} ();
      break;
    case 0:
#endif // ! STANDALONE
      first		= mtt_simpar ()(0);
      dt		= mtt_simpar ()(1);
      last		= mtt_simpar ()(2);
      stepfactor     	= mtt_simpar ()(3);
      par    		= mtt_numpar ();
      x0      		= mtt_state (par);
#ifndef STANDALONE
      break;
    default:
      usage("${sys}_ode2odes (x par simpar)", nargin);
      error("aborting.");
    }
#endif // STANDALONE
  static ColumnVector	dx (MTTNX, 0.0);
  static ColumnVector	x (MTTNX, 0.0);
  static ColumnVector	u (MTTNU, 0.0);
  static ColumnVector	y (MTTNY, 0.0);

  static Matrix	AA (MTTNX, MTTNX, 0.0);
  static ColumnVector	AAx (MTTNX, 0.0);
  
  static ColumnVector	open_switches (MTTNX, 0.0);

  register double t	= 0.0;

  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;

  for (register int i = 0; i < MTTNY; i++)
    {
      y (i) = 0.0;
    }

  for (register int i = 0; i < MTTNX; i++)
    {
      x (i) = x0 (i);
    }

  for (register int j = 0, i = 1; i <= ilast; i++)
    {
      u	= mtt_input (x, y, t, par);
      y	= mtt_output (x, u, t, par);
      if (0 == j)
	{
           mtt_write (t, x, y, nrows);
	}
EOF
case "$method" in
    "rk4")
	cat << EOF >> $filename
      {
        static ColumnVector
          k1 (MTTNX,0.0),
          k2 (MTTNX,0.0),
          k3 (MTTNX,0.0),
          k4 (MTTNX,0.0);

        const double
          t1 = t + ddt/2.0,
          t2 = t + ddt;

        ColumnVector
          x1 (x),
          x2 (x),
          x3 (x);

        k1 = ddt * mtt_rate (x , u, t , par); x1 += k1 * 0.5;
        k2 = ddt * mtt_rate (x1, u, t1, par); x2 += k2 * 0.5;
        k3 = ddt * mtt_rate (x2, u, t1, par); x3 += k3;
        k4 = ddt * mtt_rate (x3, u, t2, par);
        dx = (k1 + 2.0 * (k2 + k3) + k4) / (6.0 * ddt);
      }
EOF
    ;;
    "dassl")
    ;;
    "implicit")
	cat << EOF >> $filename
      dx = mtt_rate (x, u, t, par);
      AA = mtt_smxa (x, u, ddt, par);
      AAx = mtt_smxax (x, u, ddt, par);
EOF
    ;;
    "euler" | *)
	cat << EOF >> $filename
      dx = mtt_rate (x, u, t, par);
EOF
    ;;
esac

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

#ifdef STANDALONE
  return 0;
#else
/*  
  retval (0) = octave_value (y);
  retval (1) = octave_value (x);
  retval (2) = octave_value (t);
*/
  retval = octave_value (get_global_value ("MTT_data"));
  return (retval);
#endif // STANDALONE
}
EOF
}

case ${lang} in
    m)
        make_m	
	;;
    cc)
	make_cc
	;;
    *)
	echo Language ${lang} is not supported
esac
