#! /bin/sh

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

###############################################################
## 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
##
## Revision 1.37  2000/08/01 12:25:06  peterg

# 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)"
else
    ode=ode
    odeo=odeo
    algorithm="mtt_euler(x,dx,ddt,$Nx,open_switches)"
fi

#cat << EOF > $1_ode2odes.m
# Program $1_ode2odes
#EOF

function make_m() {
# 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
#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;
  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:$Nu
    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
    endif

  endfor;			# Integration loop

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

endfunction

EOF
} # make_m

function make_cc() {
exit

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>
### old stuff follows
if [ "$method" = "euler" ]; then
cat << EOF >> $1_ode2odes.m
ddt = mttdt/mttstepfactor; # The small sample interval
EOF

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

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

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

	MTTt = MTTt + ddt;
      endfor;
#    endif;
  ColumnVector f;
  f = feval ("${sys}_logic", args, 1)(0).vector_value ();
  return (f);
}

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
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)
      [MTTx] = $algorithm(MTTx,MTTdx,mttAA,mttAAx,mttdt,$Nx,MTTopen); # Implicit update
      MTTt = MTTt + mttdt;
#    endif;
    {
    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 ();
EOF
fi

cat << EOF >> $1_ode2odes.m
    else # NX is 0 - no states
      MTTt = MTTt + mttdt;
  endif;			# $Nx>0
      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);
endfor;				# Integration loop

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