#! /bin/sh
######################################
##### Model Transformation Tools #####
######################################
# Bourne shell script: mtt_make_sim
# Copyright (C) 2000 by Peter J. Gawthrop
###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.8 2002/04/28 18:41:27 geraint
## Fixed [ 549658 ] awk should be gawk.
## Replaced calls to awk with call to gawk.
##
## Revision 1.7 2000/05/19 17:47:35 peterg
## Fixed bug in oct version but still needs proper check
##
## Revision 1.6 2000/05/17 17:20:49 peterg
## Fixed bugs with ny>1. Could be made faster by not generating y when
## y_sim >1
##
## Revision 1.5 2000/05/17 16:01:42 peterg
## Fixed bug for n_y>1
##
## Revision 1.4 2000/05/16 18:57:15 peterg
## Still debugging
##
## Revision 1.3 2000/05/16 11:59:34 peterg
## Stard new version with data files not argv.
##
## Revision 1.2 2000/05/11 19:32:29 peterg
## Put in c version + sensitivity computation
##
## Revision 1.1 2000/04/08 10:43:26 peterg
## Initial revision
##
###############################################################
# Tell user
Sys=$1
method=$2
computation=$3
if [ -z "$method" ]; then
method=implicit
fi
if [ -n "$computation" ]; then
blurb="for language $computation"
fi
echo "Creating $1_sim.m with $method integration method $blurb"
if [ $method = "implicit" ]; then
ode=cse
odeo=cseo
else
ode=ode
odeo=odeo
fi
# Find system constants
Nx=`mtt_getsize $Sys x` # States
Nu=`mtt_getsize $Sys u` # Inputs
Ny=`mtt_getsize $Sys y` # Inputs
Npar=`wc -l $Sys\_sympar.txt | gawk '{print $1}'`
# Header
lang_header -noglobals $1 sim m 'x0,u,t,par,sensitivities' '[y,ys]' > $1_sim.m
cat >> $1_sim.m <<EOF
## tick=time;
if nargin<5
sensitivities = [];
endif;
##Sizes
N = length(t);
[n_u,N_u] = size(u);
## Doing sensitivities (assumes sensitivity system is invoked)
n_sens = length(sensitivities);
doing_sensitivities = (n_sens>0);
if doing_sensitivities
n_y = $Ny/2;
doing_state=0;
ys = [];
else
n_y = $Ny;
doing_state=(nargout>1);
endif;
## Initialise
ui = zeros(n_u,1); # Initial control
[xi] = x0; # Read in initial state
## Timing parameters
first = t(1);
dt = t(2) - t(1);
n_t = length(t);
last = t(n_t);
EOF
if [ "$computation" = "c" ]; then
cat >> $1_sim.m <<EOF
T = "";
S = "";
S = [S;"Save\t"];
tim = time;
## Create the system input file
t1 = [0:N_u-1]*dt; # Create time vector from zero (to fit u);
ut = [t1' u'];
save -text $1_input.dat ut
## Create the state file
xi = [[1:$Nx]' x0]; #'
save -text $1_state.dat xi
## Create the sympar file
save -text $1_simpar.dat dt
## Create the numpar file
[n_par,m_par] = size(par);
if m_par==1
if n_par!=$Npar
error(sprintf("Number of parameters is %i, should be %i", n_par, $Npar));
else
ipar= [[1:$Npar]' par]; #'
endif
elseif m_par==2
ipar = par;
else
error(sprintf("Number of parameter columns is %i, should be 1 or 2", m_par));
endif;
if !doing_sensitivities
save -text $1_numpar.dat ipar
endif;
dtime = time-tim;
T = [T; num2str(dtime)];
## main simulation loop
n_sim = max(1,n_sens);
for i_sim=1:n_sim
i_1 = 2;
if doing_state
i_2 = i_1 + n_y - 1;
i_3 = i_1 + n_y + 1;
i_4 = i_1 + n_y + $Nx;
else
if doing_sensitivities
i_2 = 1 + 2*n_y;
ipars = [ipar; [sensitivities(i_sim) 1]];
save -text $1_numpar.dat ipars
else
i_2 = i_1 + n_y -1;
endif
endif;
if doing_state # Need to cut twice
command = sprintf("./$1_ode2odes.out< $1_input.dat | cut -f %i-%i,%i-%i | tail -%i;", i_1,i_2,i_3,i_4,n_t);
else # not doing_state
command = sprintf("./$1_ode2odes.out< $1_input.dat | cut -f %i-%i | tail -%i;", i_1,i_2,n_t);
endif
## Execute external programme
S = [S;sprintf("Run %i\t",i_sim)];
tim = time;
yy_str=system(command);
dtime = time-tim;
T = [T; num2str(dtime)];
## Convert data
S = [S;sprintf("Conv %i\t",i_sim)];
tim = time;
yy = str2num(yy_str)'; #'
dtime = time-tim;
T = [T; num2str(dtime)];
[N_yy,M_yy] = size(yy);
if i_sim==1
if doing_sensitivities
y = yy(1:n_y:2*n_y,:); # Output
ys = yy(2:n_y:2*n_y,:);
else
y = yy(1:n_y,:);
endif;
if doing_state
i_1 = n_y+1;
i_2 = i_1 + $Nx - 1;
ys = yy(i_1:i_2,:);
endif;
else # i_sim>1
if doing_sensitivities
ys = [ys; yy(2:n_y:2*n_y,:)];
endif;
endif
endfor;
##RealTime = [S T]
endfunction
EOF
else
cat >> $1_sim.m <<EOF
A = zeros($Nx,$Nx);
Ax = zeros($Nx,1);
dx = zeros($Nx,1);
## Step size
dt = t(2)-t(1);
iFirst = first/dt;
for i = 1:N
ti = t(i);
ui = u(:,i);
yi = $1_cseo(xi,ui,ti,par); # Output
y(:,i) = yi;
x(:,i) = xi;
dxi = $1_cse(xi,ui,ti,par); # State derivative
A = $1_smxa(xi,ui,dt,par); # (I-Adt)
#A = reshape(A,$Nx,$Nx);
Ax = $1_smxax(xi,ui,dt,par); # (I-Adt)x
#open = eval(sprintf("%s_switchopen(x);", system_name)); # Open switches
#x = mtt_implicit(x,dx,A,Ax,dt,$Nx,zeros(20,1)); # Implicit update
xi = A\(Ax + dxi*dt); # Implicit update
endfor;
endfunction
EOF
fi