.POSIX:

ifeq (0,1)
CC=g++

CFLAGS+=-static
OCTAVE_LIBS_PATH=-L/usr/local/lib/octave
OCTAVE_LIBRARIES=-loctave -lcruft -lm -lncurses -ldl -lstdc++ -lc -lkpathsea -lreadline -lf2c

%.exe: %.cc
	$(CC) $< $(CFLAGS) $(OCTAVE_LIBS_PATH) $(OCTAVE_LIBRARIES) -o $@

%.oct: %.cc
	mkoctfile $<
endif



all:	$(sys)_$(rep).cc

parser: parse_m2cc.cc
	g++ $< $(CFLAGS) -o parser

ifeq ($(rep),ode2odes)
$(sys)_ode2odes.cc: ode2odes.cc
	cat ode2odes.cc | sed 's/\$$/$(sys)/' > $(sys)_ode2odes.cc
else
$(sys)_$(rep).cc: $(sys)_def.h $(sys)_sympar.h $(sys)_def.r $(sys)_$(rep).m
	mtt_m2cc.sh ${sys} ${rep} parser
endif

$(sys)_def.h: $(sys)_def.m $(sys)_sympar.txt
	def_m2h.sh $(sys)

$(sys)_sympar.h: $(sys)_sympar.txt
	sympar_txt2h.sh $(sys)

$(sys)_def.m:
	mtt $(sys) def m

# getsize needs def.r
$(sys)_def.r:
	mtt $(sys) def r

$(sys)_sympar.txt:
	mtt $(sys) sympar txt

#! /bin/sh
# $Id$
# $Log$
# Revision 1.3  2000/12/05 12:13:52  peterg
# Changed function name to name()
#
# Revision 1.2  2000/12/04 12:04:46  peterg
# Changed $() to `` for sh compatibility -- geraint
#
# Revision 1.1  2000/12/04 12:02:23  peterg
# Initial revision
#
# Revision 1.1  2000/10/31 04:32:28  geraint
# Initial revision
#

SYS=$1
IN=${SYS}_def.m
SYM=${SYS}_sympar.txt
OUT=${SYS}_def.h

get_array_size ()
{
vec=$1
awk -v vec=${vec} '($1 == vec && $2 == "=") { print $3 }' | sed s/\;//
}

 echo "// ${SYS}_def.h, generated by MTT on `date`"			>  ${OUT}
 echo ""								>> ${OUT}
 echo "const int MTTNU   = `cat ${IN} | get_array_size nu`;"		>> ${OUT}
 echo "const int MTTNX   = `cat ${IN} | get_array_size nx`;"		>> ${OUT}
 echo "const int MTTNY   = `cat ${IN} | get_array_size ny`;"		>> ${OUT}
 echo "const int MTTNZ   = `cat ${IN} | get_array_size nz`;"		>> ${OUT}
 echo "const int MTTNYZ  = `cat ${IN} | get_array_size nyz`;"		>> ${OUT}
 echo "const int MTTNPAR = `wc -l ${SYM} | awk '{ print $1 }'`;"       	>> ${OUT}

cat <<EOF >> ${OUT}

// typedefs won't work because it is illegal to initialise ColumnVector in typedef
// use "ColumnVector mttx (MTTNX);" until the proper classes are ready


#if 0 // NB: These classes are not ready for use yet!
class InputVector : public ColumnVector
{
public:
  InputVector (void) : ColumnVector (MTTNU) { ; }
};
class StateVector : public ColumnVector
{
public:
  StateVector (void) : ColumnVector (MTTNX) { ; }
};
class OutputVector : public ColumnVector
{
public:
  OutputVector (void) : ColumnVector (MTTNY) { ; }
};
class ParameterVector : public ColumnVector
{
public:
  ParameterVector (void) : ColumnVector (MTTNPAR) { ; }
};
class StateMatrix : public Matrix
{
public:
  StateMatrix (void) : Matrix (MTTNX, MTTNX) { ; }
};
#endif

EOF

/* $Id$
 * $Log$
 * Revision 1.1  2000/10/31 04:29:11  geraint
 * Initial revision
 *
 */



#include <iostream>
#include <string>



class Bracket
{
public:
  friend ostream &operator<< (ostream &str, Bracket *b);
  friend string  &operator<< (string  &str, Bracket *b);
  friend string  &operator+= (string  &str, Bracket *b);
  friend string   operator+  (string  &str, Bracket *b);
  virtual int get_nesting_depth (void)		= 0;

protected:
  Bracket (char left, char right)	  	{ _l = left; _r = right; }
  ostream &open  (ostream &str)			{ this->increment_nesting (); return str << this->_l; }
  string  &open  (string  &str)			{ this->increment_nesting (); return str += this->_l; }
  ostream &close (ostream &str)			{ this->decrement_nesting (); return str << this->_r; }
  string  &close (string  &str)			{ this->decrement_nesting (); return str += this->_r; }

protected:
  virtual int increment_nesting (void)		= 0;
  virtual int decrement_nesting (void)		= 0;
  virtual ostream &display (ostream &str)	= 0;
  virtual string  &display (string  &str)	= 0;

private:
  char _l, _r;
};



class Brace : public Bracket
{
public:
  int get_nesting_depth (void)			{ return   (this->nesting_depth);}
  string indentation (int n = nesting_depth)
  {
    string s = "";
    int i;
    for (i = 0; i < n; i++)
      s += "   ";
    return s;
  }

protected:
  Brace () : Bracket ('{', '}')			{ ; }
  int increment_nesting (void)			{ return ++(this->nesting_depth); }
  int decrement_nesting (void)			{ return --(this->nesting_depth); }
  ostream &open (ostream &str)
  {
    str << endl << this->indentation ();
    this->Bracket::open (str);
    str << endl << this->indentation ();
    return str;
  }
  string &open (string &str)
  {
    str += '\n' + this->indentation ();
    this->Bracket::open (str);
    str += '\n' + this->indentation ();
    return str;
  }
  ostream &close (ostream &str)
  {
    str << endl << this->indentation (nesting_depth - 1);
    this->Bracket::close (str);
    str << endl << this->indentation ();
    return str;
  }
  string &close (string &str)
  {
    str += '\n' + this->indentation (nesting_depth - 1);
    this->Bracket::close (str);
    str += '\n' + this->indentation ();
    return str;
  }

private:
  static int nesting_depth;
};
int Brace::nesting_depth;



class Paren : public Bracket
{
public:
  int get_nesting_depth (void)			{ return   (this->nesting_depth);}

protected:
  Paren () : Bracket ('(', ')')			{ ; }
  int increment_nesting (void)			{ return ++(this->nesting_depth); }
  int decrement_nesting (void)			{ return --(this->nesting_depth); }

private:
  static int nesting_depth;
};
int Paren::nesting_depth;



class SquareBracket : public Bracket
{
public:
  int get_nesting_depth (void)			{ return   (this->nesting_depth);}

protected:
  SquareBracket () : Bracket ('[', ']')	       	{ ; }
  int increment_nesting (void)			{ return ++(this->nesting_depth); }
  int decrement_nesting (void)			{ return --(this->nesting_depth); }

private:
  static int nesting_depth;
};
int SquareBracket::nesting_depth;



class LeftBrace : public Brace
{
public:
  ostream &display (ostream &str)		{ return this->open  (str); }
  string  &display (string  &str)		{ return this->open  (str); }
};



class RightBrace : public Brace
{
public:
  ostream &display (ostream &str)		{ return this->close (str); }
  string  &display (string  &str)		{ return this->close (str); }
};



class LeftParen : public Paren
{
public:
  ostream &display (ostream &str)		{ return this->open  (str); }
  string  &display (string  &str)		{ return this->open  (str); }
};



class RightParen : public Paren
{
public:
  ostream &display (ostream &str)		{ return this->close (str); }
  string  &display (string  &str)		{ return this->close (str); }
};



class LeftSquareBracket : public SquareBracket
{
public:
  ostream &display (ostream &str)		{ return this->open  (str); }
  string  &display (string  &str)		{ return this->open  (str); }
};



class Rightbrace : public SquareBracket
{
public:
  ostream &display (ostream &str)		{ return this->close (str); }
  string  &display (string  &str)		{ return this->close (str); }
};



ostream &operator<< (ostream &str, Bracket *b)
{
  b->display (str);
  return str;
}



string  &operator<< (string  &str, Bracket *b)
{
  string s;
  b->display (s);
  str += s;
  return str;
}



string  &operator+= (string  &str, Bracket *b)
{
  b->display (str);
  return str;
}



string   operator+  (string  &str, Bracket *b)
{
  string s;
  b->display (s);
  return (str + s);
}

// $Id$
// $Log$
// Revision 1.1  2000/11/28 04:50:29  geraint
// Initial revision
//

inline Matrix
ones (const int r = 1, const int c = 1)
{
  Matrix m (r, c);
  register int i, j;
  for (i = 0; i < r; i++)
    for (j = 0; j < c; j++)
      m (i, j) = 1.0;
  return m;
}

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

inline ColumnVector
zeros (const int r)
{
  ColumnVector v (r);
  register int i;
  for (i = 0; i < r; i++)
    v (i) = 0.0;
  return v;
}

inline Matrix
zeros (const int r, const int c)
{
  Matrix m (r, c);
  register int i, j;
  for (i = 0; i < r; i++)
    for (j = 0; j < c; j++)
      m (i, j) = 0.0;
  return m;
}

template <class T>
inline int
sign (T x)
{
  return
    (0 < x) ? +1 :
    (0 > x) ? -1 :
    0;
}

#! /bin/sh

SYS=$1
REP=$2
PARSER=${3:-indent}

IN=${SYS}_${REP}.m
OUT=${SYS}_${REP}.cc
TMP=${SYS}_${REP}_m2cc.tmp

rep_declarations ()
{
(case ${REP} in
    simpar)
	cat <<EOF
  struct {
    double first;
    double dt;
    double last;
    int    stepfactor;
    int    wmin;
    int    wmax;
    int    wsteps;
    int    input;
  } mttsimpar;

EOF
	;;
    *)
	;;
esac)
echo ""
};

rep_footer ()
{
(case ${REP} in
    simpar)
	cat <<EOF
  mttsimpar_map ["first"]      = (double) mttsimpar.first;
  mttsimpar_map ["dt"]         = (double) mttsimpar.dt;
  mttsimpar_map ["last"]       = (double) mttsimpar.last;
  mttsimpar_map ["stepfactor"] = (double) mttsimpar.stepfactor;
  mttsimpar_map ["wmin"]       = (double) mttsimpar.wmin;
  mttsimpar_map ["wmax"]       = (double) mttsimpar.wmax;
  mttsimpar_map ["wsteps"]     = (double) mttsimpar.wsteps;
  mttsimpar_map ["input"]      = (double) mttsimpar.input;
EOF
	;;
    *)
	;;
esac)
};

find_code ()
{
    file_in=${1:-${IN}}
    portion=${2:-"body"}
    head=`cat ${file_in} | awk '($2 == "BEGIN" && $3 == "Code") { print NR }'`
    foot=`cat ${file_in} | awk '($2 == "END"   && $3 == "Code") { print NR }'`
    case ${portion} in
	head)
	    start=0
	    end=${head}
	    ;;
	body)
	    start=${head}
	    end=${foot}
	    ;;
	foot)
	    start=${foot}
	    end=end
	    ;;
	*)
	    echo "Error in find_code: portion unknown"
	    return -1
	    ;;
    esac
    cat ${file_in} |\
    awk --assign start=${start} --assign end=${end} '
	(start < NR && NR < end) { print $0 }'
};


strip_junk ()
{
case ${REP} in
    numpar)
	grep -v "mttpar = \[\]"
	;;
    switchopen)
	grep -v "\[open\] = zero"
	;;
    *)
	cat
	;;
esac
};

fix_comment_delimiter ()
{
    # it would be preferable if we didn't use '%' as a delimiter
    # (a % b) gives the remainder of ((int)a / (int)b) in C/C++
    sed 's/[%#]/\/\//g' | sed 's/\/\/\/\//\/\//g'
    
};

decrement_indices ()
{    
    # first section appends '-1' to container indices
    # to convert from FORTRAN-type numbering to C-type numbering
    sed 's/mtta(\([0-9][0-9]*\),\([0-9][0-9]*\))/mtta(\1-1,\2-1)/g'	|\
    sed 's/mttax(\([0-9][0-9]*\))/mttax(\1-1)/g'			|\
    sed 's/mttdx(\([0-9][0-9]*\))/mttdx(\1-1)/g'			|\
    sed 's/mttedx(\([0-9][0-9]*\))/mttedx(\1-1)/g'			|\
    sed 's/mttpar(\([0-9][0-9]*\))/mttpar(\1-1)/g'			|\
    sed 's/mttu(\([0-9][0-9]*\))/mttu(\1-1)/g'				|\
    sed 's/mttx(\([0-9][0-9]*\))/mttx(\1-1)/g'				|\
    sed 's/mtty(\([0-9][0-9]*\))/mtty(\1-1)/g'				|\
    sed 's/mttyz(\([0-9][0-9]*\))/mttyz(\1-1)/g'			|\
    sed 's/mttz(\([0-9][0-9]*\))/mttz(\1-1)/g'				|\
    sed 's/mttopen(\([0-9][0-9]*\))/mttopen(\1-1)/g'				|\
									 \
    # next sections tidy the code up a bit, but are not necessary
    sed 's/1\-1\([\,\)]\)/0\1/g'       			       		|\
    sed 's/2\-1\([\,\)]\)/1\1/g'       			       		|\
    sed 's/3\-1\([\,\)]\)/2\1/g'       			       		|\
    sed 's/4\-1\([\,\)]\)/3\1/g'       			       		|\
    sed 's/5\-1\([\,\)]\)/4\1/g'       			       		|\
    sed 's/6\-1\([\,\)]\)/5\1/g'       			       		|\
    sed 's/7\-1\([\,\)]\)/6\1/g'       			       		|\
    sed 's/8\-1\([\,\)]\)/7\1/g'       			       		|\
    sed 's/9\-1\([\,\)]\)/8\1/g'       			       		|\
									 \
    sed 's/10\-1\([\,\)]\)/09\1/g'     			       		|\
    sed 's/20\-1\([\,\)]\)/19\1/g'	       		       		|\
    sed 's/30\-1\([\,\)]\)/29\1/g'		       	       		|\
    sed 's/40\-1\([\,\)]\)/39\1/g'			       		|\
    sed 's/50\-1\([\,\)]\)/49\1/g'		               		|\
    sed 's/60\-1\([\,\)]\)/59\1/g'		      	       		|\
    sed 's/70\-1\([\,\)]\)/69\1/g'		       	       		|\
    sed 's/80\-1\([\,\)]\)/79\1/g'		       	       		|\
    sed 's/90\-1\([\,\)]\)/89\1/g'		       	       		|\
									 \
    sed 's/100\-1\([\,\)]\)/099\1/g'   			       		|\
    sed 's/200\-1\([\,\)]\)/199\1/g'	       		       		|\
    sed 's/300\-1\([\,\)]\)/299\1/g'		       	       		|\
    sed 's/400\-1\([\,\)]\)/399\1/g'			       		|\
    sed 's/500\-1\([\,\)]\)/499\1/g'			       		|\
    sed 's/600\-1\([\,\)]\)/599\1/g'			       		|\
    sed 's/700\-1\([\,\)]\)/699\1/g'			       	       	|\
    sed 's/800\-1\([\,\)]\)/799\1/g'			       		|\
    sed 's/900\-1\([\,\)]\)/899\1/g'			       		|\
    sed 's/\([(,]\)0\([^,)]\)/\1\2/g'
};

fix_pow ()
{
    # matches number^number where number is one or more digits and one or zero decimal points
    # converts to pow (number, number)
    sed 's/\([0-9]*\)\(\.\)\{0,1\}\([0-9]*\)\^\([0-9]*\)\(\.\)\{0,1\}\([0-9]*\)/pow \(\1\2\3,\4\5\6\)/g'
};


echo Creating ${OUT}

mtt_header ${SYS} ${REP} "oct"	>  ${TMP}
find_code ${TMP} head		>  ${OUT}
rep_declarations		>> ${OUT}
find_code ${IN} body   		|\
        decrement_indices	|\
	fix_comment_delimiter	|\
	fix_pow			|\
	strip_junk		|\
	${PARSER}      		>> ${OUT}
rep_footer			>> ${OUT}
find_code ${TMP} foot		>> ${OUT}
rm ${TMP}

#include <octave/oct.h>

#include <octave/toplev.h>
#include <octave/LSODE.h>
#include <octave/ov-struct.h>
#include <octave/oct-map.h>

#include "$_def.h"
#include "$_sympar.h"

octave_value_list
mtt_cse (ColumnVector x, ColumnVector u, double t, 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);
  f = feval ("$_cse", args, 2);
  return (f);
}

ColumnVector
mtt_cseo (ColumnVector x, ColumnVector u, double t, 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);
  ColumnVector f;
  f = feval ("$_cseo", args, 1)(0).vector_value ();
  return (f);
}

#define mtt_implicit(x,dx,AA,AAx,ddt,nx,open) call_mtt_implicit((x),(dx),(AA),(AAx),(ddt),(nx),(open))
ColumnVector
call_mtt_implicit (ColumnVector x,
		   ColumnVector dx,
		   Matrix AA,
		   ColumnVector AAx,
		   double ddt,
		   int nx,
		   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);
  f = feval ("mtt_implicit", args, 1);
  return f(0).vector_value ();
}


ColumnVector
mtt_input (ColumnVector x, ColumnVector y, const double t, 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);
  ColumnVector f;
  f = feval ("$_input", args, 1)(0).vector_value ();
  return (f);
}

ColumnVector
mtt_numpar (void)
{
  octave_value_list args;
  ColumnVector f;
  f = feval ("$_numpar", args, 1)(0).vector_value ();
  return (f);
}

Octave_map
mtt_simpar (void)
{
  octave_value_list args;
  Octave_map f;
  f["first"]		= feval ("$_simpar", args, 1)(0).map_value ()["first"];
  f["dt"]		= feval ("$_simpar", args, 1)(0).map_value ()["dt"];
  f["last"]		= feval ("$_simpar", args, 1)(0).map_value ()["last"];
  f["stepfactor"]     	= feval ("$_simpar", args, 1)(0).map_value ()["stepfactor"];
  f["wmin"]		= feval ("$_simpar", args, 1)(0).map_value ()["wmin"];
  f["wmax"]		= feval ("$_simpar", args, 1)(0).map_value ()["wmax"];
  f["wsteps"]		= feval ("$_simpar", args, 1)(0).map_value ()["wsteps"];
  f["input"]		= feval ("$_simpar", args, 1)(0).map_value ()["input"];
  return (f);
}

Matrix
mtt_smxa (ColumnVector x, ColumnVector u, double t, 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);
  Matrix f;
  f = feval ("$_smxa", args, 1)(0).matrix_value ();
  return (f);
}

ColumnVector
mtt_smxax (ColumnVector x, ColumnVector u, double t, 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);
  ColumnVector f;
  f = feval ("$_smxa", args, 1)(0).vector_value ();
  return (f);
}

ColumnVector
mtt_state (ColumnVector x)
{
  octave_value_list args;
  args (0) = octave_value (x);
  ColumnVector f;
  f = feval ("$_state", args, 1)(0).vector_value ();
  return (f);
}

ColumnVector
mtt_switchopen (ColumnVector x)
{
  octave_value_list args;
  args (0) = octave_value (x);
  ColumnVector f;
  f = feval ("$_switchopen", args, 1)(0).vector_value ();
  return (f);
}

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)
    : 0x0;
}


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

  ColumnVector	x;
  ColumnVector	par;
  Octave_map	simpar;

  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"];
      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 ();
      break;
    case 2:
      simpar["first"]		= mtt_simpar ()["first"];
      simpar["dt"]		= mtt_simpar ()["dt"];
      simpar["last"]		= mtt_simpar ()["last"];
      simpar["stepfactor"]     	= mtt_simpar ()["stepfactor"];
      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 ();
      break;
    case 1:
      simpar["first"]		= mtt_simpar ()["first"];
      simpar["dt"]		= mtt_simpar ()["dt"];
      simpar["last"]		= mtt_simpar ()["last"];
      simpar["stepfactor"]     	= mtt_simpar ()["stepfactor"];
      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 ();
      break;
    case 0:
      simpar["first"]		= mtt_simpar ()["first"];
      simpar["dt"]		= mtt_simpar ()["dt"];
      simpar["last"]		= mtt_simpar ()["last"];
      simpar["stepfactor"]     	= mtt_simpar ()["stepfactor"];
      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);
      break;
    default:
      usage("$_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	= (int)round (simpar ["last"].double_value () / ddt);

  // cse translation
  // LSODE will need ODEFUNC

  for (register int j = 0, i = 1; i <= ilast; i++)
    {
      y	= mtt_cseo (x, u, t, par);
      u	= mtt_input (x, y, t, par);
      if (0 == j)
	{
	  mtt_write (t, x, y, MTTNX, MTTNY);
	}
      dx = mtt_cse (x, u, t, par)(0).vector_value ();
      
      AA = mtt_smxa (x, u, ddt, par);
      AAx = mtt_smxax (x, u, ddt, par);

      open_switches = mtt_switchopen (x);
      x = mtt_implicit (x, dx, AA, AAx, ddt, 1, open_switches);
      t += ddt;
      j++;
      j = (j == (int)simpar ["stepfactor"].double_value ()) ? j : 0;
    }

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

/* $Id$
 * $Log$
 * Revision 1.1  2000/10/31 04:29:50  geraint
 * Initial revision
 *
 */



#include <iostream>
#include <map>
#include <stack>
#include <string>



/*
 * Bracket.hh deals with nesting levels of parenthesis
 * just add Bracket pointer to string / stream
 */
#include "Bracket.hh"



/*
 * use lbrace, etc. in expressions to automate nesting calculations
 */
LeftBrace	*lbrace = new LeftBrace;
RightBrace	*rbrace = new RightBrace;

LeftParen	*lparen = new LeftParen;
RightParen	*rparen = new RightParen;



/*
 * use brace nesting depth to determine indentation
 */
string indent (void) { return lbrace->indentation (); }



/*
 * map contains keyword to look for and function to call when found
 */
map <string, void (*)(void)> keyword;




/*
 * stack records current nest type
 * so that "end" can be handled correctly
 */
enum nest_type
{
  null,
  if_statement,
  switch_statement,
  while_statement
};

stack <enum nest_type> current_nest;



/*
 * assuming cin has '(' as next character
 * find the matching ')' and return the string contained within
 */
string get_test (void)
{
  const int entry_nesting = lparen->get_nesting_depth ();
  char c;
  string buf = "";
  cin.setf (ios::skipws);
  cin >> c;
  if ('(' != c)
    {
      cerr << "Sorry, test must be enclosed in ( )" << endl;
      cerr << "current character is " << c << endl;
      /*
       * replace this with a call to get_expression maybe?
       */
      exit (-1);
    }
  else
    {
      buf += lparen;
    }
  cin.unsetf (ios::skipws);
  while (entry_nesting != lparen->get_nesting_depth ()) {
    cin >> c;
    switch (c)
      {	
      case '(':
	buf += lparen;
	break;
      case ')':
	buf += rparen;
	break;
      case EOF:
	cerr << "Oops! EOF reached" << endl;
	exit (-1);
      default:
	buf += c;
	break;
      }
  }
  return buf;
}



/*
 * functions to be called upon detection of keyword
 */
void got_if (void)
{
  current_nest.push (if_statement);
  cout << "if ";
  cout << get_test ();
  cout << lbrace;
}

void got_else (void)
{
  cout << rbrace;
  cout << "else";
  cout << lbrace;
}

void got_elseif (void)
{
  cout << rbrace;
  got_if ();
} 

void got_switch (void)
{
  current_nest.push (switch_statement);
  cout << "switch ";
  cout << get_test ();
  cout << lbrace;
  /*
   * open a second brace so that each "case" can be enclosed
   */
  cout << lbrace;
}

void got_case (void)
{
  cout << rbrace;
  cout << "case ";
  cout << get_test ();
  cout << ':';
  cout << lbrace;
}

void got_otherwise (void)
{
  cout << rbrace;
  cout << "default:";
  cout << lbrace;
}

void got_while (void)
{
  current_nest.push (while_statement);
  cout << "while ";
  cout << get_test ();
  cout << lbrace;
}

void got_end (void)
{
  enum nest_type n = current_nest.top ();
  switch (n)
    {
    case if_statement:
    case while_statement:
      cout << rbrace;      
      break;
    case switch_statement:
      cout << rbrace
	   << rbrace;
      break;
    default:
      cerr << "Oops! Unmatched end ... aborting" << endl;
      exit (-1);
      break;
    }
  current_nest.pop ();
}

void got_comment (void)
{
  /*
   * get remainder of line in case there are any commented keywords
   */
  char c;
  cout << " // ";
  cin >> c;
  while (c != '\n' && c != EOF)
    {
      cout << c;
      cin >> c;
    }
  cout << endl
       << indent () << ';' << endl
       << indent ();
}
  


/*
 * map contains keyword to look for and functions to call when found
 * functions must be of the form "void f (void) {...;}"
 * new structures should also be added to enum nest_type
 */
void set_keyword (void)
{
  // if
  keyword ["if"]		= got_if;
  keyword ["else"]		= got_else;
  keyword ["elseif"]		= got_elseif;
  keyword ["endif"]		= got_end;

  // switch
  keyword ["switch"]		= got_switch;
  keyword ["case"]		= got_case;
  keyword ["otherwise"] 	= got_otherwise;
  keyword ["endswitch"] 	= got_end;

  // while
  keyword ["while"]		= got_while;
  keyword ["endwhile"]		= got_end;

  // general
  keyword ["end"]		= got_end;

  /*
   * this won't work inside a string
   * or if there is no space around the delimiter
   */
  keyword ["#"]			= got_comment;
  keyword ["//"]		= got_comment;
  
/*
 * unimplemented keywords
 *
 keyword ["for"]		= got_for;
 keyword ["endfor"]		= got_end;
 keyword ["break"]		= got_break;
 keyword ["continue"]		= got_continue;
*/
}



/* read in one character at a time looking for tokens (delimited by whitespace)
 * if the token is a keyword, call the appropriate function
 */
int main (void)
{
  set_keyword ();
  string buf = "";
  char c;
  cin.unsetf (ios::skipws);
  while (cin >> c)
    {
      switch (c)
	{
	case EOF:
	  return 0;
	case ' ':
	case '\t':
	  if (keyword [buf])
	    {
	      keyword [buf]();
	    }
	  else
	    {
	      if (! lbrace->get_nesting_depth ())
		{
		  buf += c;
		}
	      cout << buf;
	    }
	  buf = "";
	  break;
	case '\n':
	  if (keyword [buf])
	    {
	      /*
	       * keyword found, call function
	       */
	      keyword [buf] ();
	    }
	  else
	    {
	      cout << buf
		/*
		 * keep newline in case this line has an EOL-comment
		 */
		   << endl
		/*
		 * EOL is end-of-statement in Octave, add ;
		 */
		   << indent () << ';' << endl
		   << indent ();
	    }
	  buf = "";
	  break;
	default:
	  buf += c;
	}
      if (lbrace->get_nesting_depth ())
	{
	  cout.setf (ios::skipws);
	}
      else
	{
	  cout.unsetf (ios::skipws);
	}
    }
  return 0;
}

#include <octave/oct.h>

#include <octave/toplev.h>
#include <octave/LSODE.h>
#include <octave/ov-struct.h>
#include <octave/oct-map.h>

#include "rc_def.h"
#include "rc_sympar.h"

octave_value_list
mtt_csex (ColumnVector x, ColumnVector u, double t, 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);
  f = feval ("rc_csex", args, 2);
  return (f);
}

ColumnVector
mtt_cseo (ColumnVector x, ColumnVector u, double t, 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);
  ColumnVector f;
  f = feval ("rc_cseo", args, 1)(0).vector_value ();
  return (f);
}

#define mtt_implicit(x,dx,AA,AAx,ddt,nx,open) call_mtt_implicit((x),(dx),(AA),(AAx),(ddt),(nx),(open))
ColumnVector
call_mtt_implicit (ColumnVector x,
		   ColumnVector dx,
		   Matrix AA,
		   ColumnVector AAx,
		   double ddt,
		   int nx,
		   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);
  f = feval ("mtt_implicit", args, 1);
  return f(0).vector_value ();
}


ColumnVector
mtt_input (ColumnVector x, ColumnVector y, const double t, 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);
  ColumnVector f;
  f = feval ("rc_input", args, 1)(0).vector_value ();
  return (f);
}

ColumnVector
mtt_numpar (void)
{
  octave_value_list args;
  ColumnVector f;
  f = feval ("rc_numpar", args, 1)(0).vector_value ();
  return (f);
}

Octave_map
mtt_simpar (void)
{
  octave_value_list args;
  Octave_map f;
  f["first"]		= feval ("rc_simpar", args, 1)(0).map_value ()["first"];
  f["dt"]		= feval ("rc_simpar", args, 1)(0).map_value ()["dt"];
  f["last"]		= feval ("rc_simpar", args, 1)(0).map_value ()["last"];
  f["stepfactor"]     	= feval ("rc_simpar", args, 1)(0).map_value ()["stepfactor"];
  f["wmin"]		= feval ("rc_simpar", args, 1)(0).map_value ()["wmin"];
  f["wmax"]		= feval ("rc_simpar", args, 1)(0).map_value ()["wmax"];
  f["wsteps"]		= feval ("rc_simpar", args, 1)(0).map_value ()["wsteps"];
  f["input"]		= feval ("rc_simpar", args, 1)(0).map_value ()["input"];
  return (f);
}

Matrix
mtt_smxa (ColumnVector x, ColumnVector u, double t, 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);
  Matrix f;
  f = feval ("rc_smxa", args, 1)(0).matrix_value ();
  return (f);
}

ColumnVector
mtt_smxax (ColumnVector x, ColumnVector u, double t, 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);
  ColumnVector f;
  f = feval ("rc_smxa", args, 1)(0).vector_value ();
  return (f);
}

ColumnVector
mtt_state (ColumnVector x)
{
  octave_value_list args;
  args (0) = octave_value (x);
  ColumnVector f;
  f = feval ("rc_state", args, 1)(0).vector_value ();
  return (f);
}

ColumnVector
mtt_logic (ColumnVector x, ColumnVector u, double t, 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);

  ColumnVector f;
  f = feval ("rc_logic", args, 1)(0).vector_value ();
  return (f);
}

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)
    : 0x0;
}


DEFUN_DLD (rc_ode2odes, args, ,
"Octave ode2odes representation of system $
rc_ode2odes (x, par, simpar)
")
{
  octave_value_list retval;

  ColumnVector	x;
  ColumnVector	par;
  Octave_map	simpar;

  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"];
      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 ();
      break;
    case 2:
      simpar["first"]		= mtt_simpar ()["first"];
      simpar["dt"]		= mtt_simpar ()["dt"];
      simpar["last"]		= mtt_simpar ()["last"];
      simpar["stepfactor"]     	= mtt_simpar ()["stepfactor"];
      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 ();
      break;
    case 1:
      simpar["first"]		= mtt_simpar ()["first"];
      simpar["dt"]		= mtt_simpar ()["dt"];
      simpar["last"]		= mtt_simpar ()["last"];
      simpar["stepfactor"]     	= mtt_simpar ()["stepfactor"];
      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 ();
      break;
    case 0:
      simpar["first"]		= mtt_simpar ()["first"];
      simpar["dt"]		= mtt_simpar ()["dt"];
      simpar["last"]		= mtt_simpar ()["last"];
      simpar["stepfactor"]     	= mtt_simpar ()["stepfactor"];
      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);
      break;
    default:
      usage("rc_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	= (int)round (simpar ["last"].double_value () / ddt);

  // cse translation
  // LSODE will need ODEFUNC

  for (register int j = 0, i = 1; i <= ilast; i++)
    {
      y	= mtt_cseo (x, u, t, par);
      u	= mtt_input (x, y, t, par);
      if (0 == j)
	{
	  // mtt_write (t, x, y, MTTNX, MTTNY);
	}
      dx = mtt_csex (x, u, t, par)(0).vector_value ();
      
      AA = mtt_smxa (x, u, ddt, par);
      AAx = mtt_smxax (x, u, ddt, par);

      open_switches = mtt_logic (x, u, t, par);
      x = mtt_implicit (x, dx, AA, AAx, ddt, 1, open_switches);
      t += ddt;
      j++;
      j = (j == (int)simpar ["stepfactor"].double_value ()) ? j : 0;
    }

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

#! /bin/sh
# $Id$
# $Log$
# Revision 1.5  2000/12/05 12:44:55  peterg
# Changed $() to ``
#
# Revision 1.4  2000/12/05 12:16:02  peterg
# Changed function name to name()
#
# Revision 1.3  2000/12/04 11:05:01  peterg
# Removed () -- geraint
#
# Revision 1.2  2000/11/07 17:29:27  peterg
# Changed echo
#
# Revision 1.1  2000/11/07 17:28:53  peterg
# Initial revision
#
# Revision 1.1  2000/10/31 04:32:02  geraint
# Initial revision
#

SYS=$1

if [ $# -gt 1 ]
then
  NUM_OF_TMP_VAR=$2;
  shift; shift;
else
  NUM_OF_TMP_VAR=500
  shift;
fi
TMP_VAR_NAMES="mtt_tmp mtt_o $*"

IN=${SYS}_sympar.txt
OUT=${SYS}_sympar.h

declare_sys_param ()
{
cat ${IN} | awk '{printf ("double %s;\t// %s\n", $1, $2)}'
}

declare_temp_vars ()
{
for name in ${TMP_VAR_NAMES}
do
    echo ""
    i=0
    while [ ${i} -le ${NUM_OF_TMP_VAR} ]
    do
	echo "double ${name}${i};"
	i=`expr ${i} + 1`
    done
done
}

echo Creating ${OUT}
declare_sys_param	>  ${OUT}
declare_temp_vars	>> ${OUT}