#! /bin/sh
######################################
##### Model Transformation Tools #####
######################################
###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## 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 "${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(mtt_write)
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
ColumnVector
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);
#ifdef STANDALONE
f = F${sys}_${ode} (args, 1);
#else
f = feval ("${sys}_${ode}", args, 1);
#endif
return f(0).${vector_value} ();
}
ColumnVector
mtt_${odeo} (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);
#ifdef STANDALONE
f = F${sys}_${odeo} (args, 1);
#else
f = feval ("${sys}_${odeo}", args, 1);
#endif
return f(0).${vector_value} ();
}
// #define mtt_implicit(x,dx,AA,AAx,ddt,nx,open) call_mtt_implicit((x),(dx),(AA),(AAx),(ddt),(nx),(open))
ColumnVector
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);
#ifdef STANDALONE
f = Fmtt_implicit (args, 1);
#else
f = feval ("mtt_implicit", args, 1);
#endif
return f(0).${vector_value} ();
}
ColumnVector
mtt_euler (ColumnVector x,
ColumnVector dx,
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 (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} ();
}
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);
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} ();
}
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} ();
}
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);
}
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);
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 ();
}
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);
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} ();
}
ColumnVector
mtt_state (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} ();
}
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);
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} ();
}
void
mtt_write (double t, ColumnVector x, ColumnVector y){
octave_value_list args;
args (0) = octave_value (t);
args (1) = octave_value (x);
args (2) = octave_value (y);
#ifdef STANDALONE
Fmtt_write (args, 1);
#else
feval ("mtt_write", args, 1);
#endif
}
//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)
: ColumnVector ();
}
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 = static_cast<int> (round (simpar ["last"].double_value () / ddt)) + 1;
// cse translation
// LSODE will need ODEFUNC
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, MTTNX, MTTNY);
mtt_write (t, x, y);
}
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> (simpar ["stepfactor"].double_value ())) ? 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