Index: mttroot/mtt/bin/trans/make_ode2odes
==================================================================
--- mttroot/mtt/bin/trans/make_ode2odes
+++ mttroot/mtt/bin/trans/make_ode2odes
@@ -7,10 +7,13 @@
 ###############################################################
 ## Version control history
 ###############################################################
 ## $Id$
 ## $Log$
+## 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
@@ -129,24 +132,26 @@
 # Bourne shell script: make_ode2odes
 
 # Copyright (c) P.J.Gawthrop July 1998.
 
 # Tell user
-Sys=$1
+sys=$1
+lang=$2
+filename=${sys}_ode2odes.${lang}
 
-if [ -n "$2" ]; then
-  method=$2    
+if [ -n "$3" ]; then
+  method=$3    
 else
   method=implicit  
 fi
 
-echo  "Creating $1_ode2odes.m with $method integration method"
+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  
+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)"
@@ -154,31 +159,26 @@
     ode=ode
     odeo=odeo
     algorithm="mtt_euler(x,dx,ddt,$Nx,open_switches)"
 fi
 
-#cat << EOF > $1_ode2odes.m
-# Program $1_ode2odes
-#EOF
-
-# Do the globals
-#sympar2global_txt2m $1 >> $1_ode2odes.m
-lang_header  $1 ode2odes m  'x,par,simpar' '[Y,X,t]' > $1_ode2odes.m
-
-cat >> $1_ode2odes.m <<EOF
+function 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 = $1_simpar; 
+    simpar = ${sys}_simpar; 
     [simpar.dt] = mtt_simpar_update;
   endif
   if nargin<2
-    par = $1_numpar; 
+    par = ${sys}_numpar; 
     [par] = mtt_numpar_update(par);
   endif
   if nargin<1
-    [x]  =  $1_state(par); 
+    [x]  =  ${sys}_state(par); 
     [x] = mtt_state_update(x);
   endif
 
   ## Initialise
   t = 0.0;
@@ -191,29 +191,29 @@
     u(MTTi) = 0;
   endfor;
 
   mttj = 0;
     for it = 1:ilast		#Integration loop
-    [y] = $1_$odeo(x,u,t,par);  # Output 
-    [u] = $1_input(x,y,t,par);	# Input
+    [y] = ${sys}_$odeo(x,u,t,par);  # Output 
+    [u] = ${sys}_input(x,y,t,par);	# Input
     if mttj==0
      mtt_write(t,x,y,$Nx,$Ny);   # Write it out
     endif
-    [dx] = $1_$ode(x,u,t,par);	# State derivative
+    [dx] = ${sys}_$ode(x,u,t,par);	# State derivative
 EOF
 
 if [ "$method" = "implicit" ]; then
-cat<<EOF >> $1_ode2odes.m
+cat<< EOF >> $filename
 
-    [AA] = $1_smxa(x,u,ddt,par);	# (I-Adt) and (I-Adt)x
-    [AAx] = $1_smxax(x,u,ddt,par); # (I-Adt) and (I-Adt)x
+    [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 >> $1_ode2odes.m
-    [open_switches] = $1_logic(x,u,t,par); # Switch logic
+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
@@ -226,59 +226,328 @@
   X  = MTT_data(:,4);
 
 endfunction
 
 EOF
-
-exit
-
-### old stuff follows
-if [ "$method" = "euler" ]; then
-cat << EOF >> $1_ode2odes.m
-ddt = mttdt/mttstepfactor; # The small sample interval
-EOF
-fi
-
-
-cat << EOF >> $1_ode2odes.m
-for MTTit = 1:MTTilast		#Integration loop
-  [MTTy] = $1_$odeo(MTTx,MTTu,MTTt,MTTpar);     # Output 
-  [MTTu] = $1_input(MTTt,MTTx,MTTy);	# Input
-  mtt_write(MTTt,MTTx,MTTy,$Nx,$Ny);    # Write it out
-
-  if $Nx>0			# Dont if no states
-EOF
-
-if [ "$method" = "euler" ]; then
-cat << EOF >> $1_ode2odes.m
-#    if mttmethod==1		# Euler
-      for MTTjt = 1:mttstepfactor
-	[MTTdx] = $1_$ode(MTTx,MTTu,MTTt,MTTpar); # State derivative
-        [MTTopen] = $1_logic(MTTx,MTTu,MTTt,MTTpar); # Switch logic
-	[MTTx] = mtt_euler(MTTx,MTTdx,ddt,$Nx,MTTopen); # Euler update
-	MTTt = MTTt + ddt;
-      endfor;
-#    endif;
-EOF
-fi
-
-if [ "$method" = "implicit" ]; then
-cat << EOF >> $1_ode2odes.m
-#    if mttmethod==2		# Implicit 
-      [MTTdx] = $1_$ode(MTTx,MTTu,MTTt,MTTpar); # State derivative
-      [mttAA] = $1_smxa(MTTx,MTTu,mttdt,MTTpar);	# (I-Adt) and (I-Adt)x
-      [mttAAx] = $1_smxax(MTTx,MTTu,mttdt,MTTpar);	# (I-Adt) and (I-Adt)x
-      [MTTopen] = $1_logic(MTTx,MTTu,MTTt,MTTpar); # Switch logic
-      [MTTx] = $algorithm(MTTx,MTTdx,mttAA,mttAAx,mttdt,$Nx,MTTopen); # Implicit update
-      MTTt = MTTt + mttdt;
-#    endif;
-EOF
-fi
-
-cat << EOF >> $1_ode2odes.m
-    else # NX is 0 - no states
-      MTTt = MTTt + mttdt;
-  endif;			# $Nx>0
-
-endfor;				# Integration loop
-
-EOF
+} # make_m
+
+function make_cc() {
+
+cat <<EOF  > $filename
+#include <octave/oct.h>
+
+#include <octave/toplev.h>
+#include <octave/LSODE.h>
+#include <octave/ov-struct.h>
+#include <octave/oct-map.h>
+
+#include "${sys}_def.h"
+#include "${sys}_sympar.h"
+
+octave_value_list
+mtt_${ode} (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 ("${sys}_${ode}", 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 ("${sys}_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 ("${sys}_input", args, 1)(0).vector_value ();
+  return (f);
+}
+
+ColumnVector
+mtt_numpar (void)
+{
+  octave_value_list args;
+  ColumnVector f;
+  f = feval ("${sys}_numpar", args, 1)(0).vector_value ();
+  return (f);
+}
+
+Octave_map
+mtt_simpar (void)
+{
+  octave_value_list args;
+  Octave_map f;
+  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"];
+  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 ("${sys}_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 ("${sys}_smxax", 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 ("${sys}_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 ("${sys}_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 (${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;
+
+  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("${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	= (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_${ode} (x, u, t, par)(0).vector_value ();
+EOF
+
+if [ "$method" = "implicit" ]; then
+echo Hi $filename
+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 = 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);
+}
+
+EOF
+}
+
+case ${lang} in
+    m)
+        make_m	
+	;;
+    cc)
+	make_cc
+	;;
+    *)
+	echo Language ${lang} is not supported
+esac