#! /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
## 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

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  

if [ "$method" = "implicit" ]; then
    ode=csex
    odeo=cseo
    algorithm="mtt_implicit(x,dx,AA,AAx,ddt,$Nx,open_switches)"
else
    ode=ode
    odeo=odeo
    algorithm="mtt_euler(x,dx,ddt,$Nx,open_switches)"
fi

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
    [x]  =  ${sys}_state(par); 
    [x] = mtt_state_update(x);
  endif

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

## 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
    [dx] = ${sys}_$ode(x,u,t,par);	# State derivative
EOF

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

  t = MTT_data(:,1);
  Y = MTT_data(:,2);
  X  = MTT_data(:,4);

endfunction

EOF
} # make_m

make_cc() 
{
# get octave version
octave_development=`octave --version | awk '{print $4}' | awk -F\. '{print $2}'`
if [ $octave_development ]; then
    vector_value=column_vector_value
    feval_header=parse.h
else
    vector_value=vector_value
    feval_header=toplev.h
fi

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(${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

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} ();
}

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} ();
}

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} ();
}

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} ();
}

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} ();
}

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} ();
}

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["last"]		= F${sys}_simpar (args, 1)(0).map_value ()["last"];
  f["stepfactor"]     	= F${sys}_simpar (args, 1)(0).map_value ()["stepfactor"];
  f["wmin"]		= F${sys}_simpar (args, 1)(0).map_value ()["wmin"];
  f["wmax"]		= F${sys}_simpar (args, 1)(0).map_value ()["wmax"];
  f["wsteps"]		= F${sys}_simpar (args, 1)(0).map_value ()["wsteps"];
  f["input"]		= F${sys}_simpar (args, 1)(0).map_value ()["input"];
#else
  f["first"]		= feval ("${sys}_simpar", args, 1)(0).map_value ()["first"];
  f["dt"]		= feval ("${sys}_simpar", args, 1)(0).map_value ()["dt"];
  f["last"]		= feval ("${sys}_simpar", args, 1)(0).map_value ()["last"];
  f["stepfactor"]     	= feval ("${sys}_simpar", args, 1)(0).map_value ()["stepfactor"];
  f["wmin"]		= feval ("${sys}_simpar", args, 1)(0).map_value ()["wmin"];
  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);
}

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 ();
}

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} ();
}

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} ();
}

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} ();
}


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
    cout << data;
#else // ! STANDALONE
  set_global_value ("MTT_data", data);
#endif

}


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:
      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 ();
      par    		= args (1).${vector_value} ();
      x      		= args (0).${vector_value} ();
      break;
    case 2:
      first		= mtt_simpar ()["first"].double_value ();
      dt		= mtt_simpar ()["dt"].double_value ();
      last		= mtt_simpar ()["last"].double_value ();
      stepfactor     	= mtt_simpar ()["stepfactor"].double_value ();
      par    		= args (1).${vector_value} ();
      x      		= args (0).${vector_value} ();
      break;
    case 1:
      first		= mtt_simpar ()["first"].double_value ();
      dt		= mtt_simpar ()["dt"].double_value ();
      last		= mtt_simpar ()["last"].double_value ();
      stepfactor     	= mtt_simpar ()["stepfactor"].double_value ();
      par    		= mtt_numpar ();
      x      		= args (0).${vector_value} ();
      break;
    case 0:
      first		= mtt_simpar ()["first"].double_value ();
      dt		= mtt_simpar ()["dt"].double_value ();
      last		= mtt_simpar ()["last"].double_value ();
      stepfactor     	= mtt_simpar ()["stepfactor"].double_value ();
      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	= 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 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, 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> (stepfactor)) ? 0 : j;
    }

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

#ifdef STANDALONE
int main (void)
{
  octave_value_list args;
  F${sys}_ode2odes (args, 3);
  return 0;
}
#endif
EOF
}

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