File mttroot/mtt/bin/trans/mtt_make_sim artifact d2886556d7 part of check-in f1d97fb1f0

#! /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.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 | awk '{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 -ascii $1_input.dat ut 

  ## Create the state file
  xi = [[1:$Nx]' x0]; #'
  save -ascii   $1_state.dat xi

  ## Create the sympar file
  save -ascii   $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 -ascii   $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 -ascii   $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