# Copyright (C) 2000 by Peter J. Gawthrop

###############################################################
## Version control history
###############################################################
## $Header$
## $Log$
## Revision 1.293.2.6  2001/03/03 06:50:38  geraint
## Added #SUMMARY lines for ode2odes.
##
## Revision 1.293.2.5  2001/03/03 00:27:14  geraint
## Fixed ar options to work with GNU ar. Allow mtt to create dependencies for mtt_%.cc when making mtt_%.oct.
##
## Revision 1.293.2.4  2001/03/02 00:45:21  geraint
## Separated Euler and Implicit methods in .cc code and dependencies.
##
## Revision 1.293.2.3  2001/03/01 05:05:53  geraint
## Minor revisions.
##
## Revision 1.293.2.2  2001/02/23 03:53:53  geraint
## Minor revisions. (ode2odes.exe works on GNU/Linux and ~GNU/Irix)
##
## Revision 1.293.2.1  2001/02/22 06:41:40  geraint
## Initial merge of ode2odes.exe into main mtt.
## standalone_rep.make deleted: rules moved to mtt, variables to mttrc.
##
## Revision 1.293  2001/02/06 04:39:08  geraint
## Added _standalone.exe representation.
##
## Revision 1.292  2001/02/05 17:27:40  gawthrop
## Make sure _def.r exists before creating _state.txt
##
## Revision 1.291  2000/12/27 14:50:40  peterg
## This is the first CVS version (4.9).
## Commented out code now deleted
##

			integration_method=euler;
                        mtt_switches="$mtt_switches euler";
			;;
		    implicit)
			integration_method=implicit;
                        mtt_switches="$mtt_switches implicit";
			;;
		    rk4)
			integration_method=rk4;
			mtt_switches="$mtt_switches rk4";
			;;
		    *)
			echo $1 is an unknown integration method - use euler, rk4 or implicit;
                        exit;;
		esac;;
	-s )
                mtt_switches="$mtt_switches $1";
		sensitivity=sensitivity ;;
	-ss )
                mtt_switches="$mtt_switches $1";

    echo '         -D  debug -- leave log files etc'
    echo '         -I  prints more information'
    echo '         -abg start at abg.m representation'
    echo '         -c  c-code generation'
    echo '         -d  <dir>  use directory <dir>'
    echo '         -dc Maximise derivative (not integral) causality'
    echo '         -dc Maximise derivative (not integral) causality'
    echo '         -i <implicit|euler|rk4>   Use implicit, euler or rk4 integration'
    echo '         -o ode is same as dae'
    echo '         -oct use oct files in place of m files where appropriate'
    echo '         -opt optimise code generation'
    echo '         -p  print environment variables'
    echo '         -partition partition hierachical system'
    echo '         -r  reset time stamp on representation'
    echo '         -s  try to generate sensitivity BG (experimental)'

fi


if [ -n "$sensitivity" ]; then
  sys=`echo $1 | cut -c 2-`
  sys_abg=$sys"_abg"
fi

## Check octave version
if [ `octave --version | awk -F\. '{print $2}'` ]; then
    define_octave_dev="-DOCTAVE_DEV"
else
    define_octave_dev=""
fi

################################
# This is the main mtt programme
################################


##echo Target is $target, Subsystem is $Subsystem, options are $mtt_switches
$MAKE -S $verbose  -f  -  $target << EOF

# Cancel implicit rules I don't want
%.dvi: %.tex

# MTT implicit rules
## .oct files

$1_%.oct:  $1_%.cc $1_def.h $1_sympar.h $1_cr.h
	echo Creating $1_\$*.oct; $MKOCTFILE  $1_\$*.cc


.PRECIOUS: mtt_%.oct
mtt_%.oct: mtt_%.cc
	echo Compiling \$<
	${MKOCTFILE} $define_octave_dev \$<

## .cc files
.PRECIOUS: %.cc # Don't let mtt delete them

$1_%.cc:  $1_%.m
	mtt_m2cc.sh  $1 \$* cat 

mtt_%.cc: ${MTT_LIB}/cc/mtt_%.cc
	cp ${MTT_LIB}/cc/mtt_\$*.cc mtt_\$*.cc

## .o files
.PRECIOUS: $1_%.o
$1_%.o: $1_%.cc $1_def.h $1_sympar.h $1_cr.h
	echo Compiling $1_\$*.cc
	${MTT_CXX} ${MTT_CXXFLAGS} ${MTT_CXXINCS} -c $< -DSTANDALONE

.PRECIOUS: mtt_%.o
mtt_%.o: mtt_%.cc
	echo Compiling mtt_\$*.cc
	${MTT_CXX} ${MTT_CXXFLAGS} ${MTT_CXXINCS} -c $< -DSTANDALONE $define_octave_dev

# Generates code to write reduce code
$1_ode_write.r: $1_def.m
	def2write_r  $1 $2

ifeq ("$sensitivity","sensitivity")
ifeq ("$level","0")

#SUMMARY simpar	simulation information (m)
$1_simpar.m: $1_simpar.txt $1_def.r $1_sympar.txt
	mtt_txt2m $1 simpar


#SUMMARY numpar	numerical parameter declaration -- default (txt) 
$1_numpar.txt : $1_sympar.txt $1_def.r


ifeq ($target,$1_numpar.txt)
	mtt_update $1 numpar update
else
	mtt_update $1 numpar
endif

#SUMMARY state	numerical state declaration  -- default (txt)
$1_state.txt : $1_struc.txt $1_sympar.txt $1_def.r



ifeq ($target,$1_state.txt)
	mtt_update $1 state update
else
	mtt_update $1 state
endif

#SUMMARY input	numerical state declaration  -- default (txt)
$1_input.txt : $1_struc.txt $1_sympar.txt


ifeq ($target,$1_input.txt)
	mtt_update $1 input update
else
	mtt_update $1 input
endif

#SUMMARY logic	Dynamic switch logic (txt)
$1_logic.txt : $1_def.r $1_struc.txt $1_sympar.txt



ifeq ($target,$1_logic.txt)
	mtt_update $1 logic update
else
	mtt_update $1 logic
endif

# Dummy target
FORCE:

$1_ode2odes_common.m    : $1_input.m   $1_logic.m   $1_numpar.m   $1_simpar.m   $1_state.m
	@echo > /dev/null
$1_ode2odes_common.cc   : $1_input.cc  $1_logic.cc  $1_numpar.cc  $1_simpar.cc  $1_state.cc
	@echo > /dev/null
$1_ode2odes_common.o    : $1_input.o   $1_logic.o   $1_numpar.o   $1_simpar.o   $1_state.o
	@echo "Creating $1_ode2odes_common.o"
	ar -cr \$@ \$^
$1_ode2odes_common.oct  : $1_input.oct $1_logic.oct $1_numpar.oct $1_simpar.oct $1_state.oct
	@echo > /dev/null

$1_ode2odes_euler.m  $1_ode2odes_rk4.m    : $1_ode.m     $1_odeo.m
	@echo > /dev/null
$1_ode2odes_euler.cc $1_ode2odes_rk4.cc   : $1_ode.cc    $1_odeo.cc
	@echo > /dev/null
$1_ode2odes_euler.o  $1_ode2odes_rk4.o    : $1_ode.o     $1_odeo.o   mtt_euler.o
	@echo "Creating \$@"
	ar -cr \$@ \$^
$1_ode2odes_euler.oct   : $1_ode.oct   $1_odeo.oct mtt_euler.oct
	@echo > /dev/null

$1_ode2odes_implicit.m  : $1_cseo.m    $1_csex.m   $1_smxa.m   $1_smxax.m
	@echo > /dev/null
$1_ode2odes_implicit.cc : $1_cseo.cc   $1_csex.cc  $1_smxa.cc  $1_smxax.cc
	@echo > /dev/null
$1_ode2odes_implicit.o  : $1_cseo.o    $1_csex.o   $1_smxa.o   $1_smxax.o   mtt_implicit.o
	@echo "Creating $1_ode2odes_implicit.o"
	ar -cr \$@ \$^
$1_ode2odes_implicit.oct: $1_cseo.oct  $1_csex.oct $1_smxa.oct $1_smxax.oct mtt_implicit.oct
	@echo > /dev/null


#SUMMARY numpar	numerical parameter declaration (m) 
$1_numpar.m:  $1_numpar.txt $1_sympars.txt
	mtt_txt2m $1 numpar

#SUMMARY numpar	numerical parameter declaration (c) 
#SUMMARY numpar	numerical parameter declaration (view) 

	touch $1_aliased.txt
$1_struc.tex: $1_struc.txt
	struc_txt2tex $1
$1_struc.m: $1_struc.txt 
	struc_txt2m $1
$1_def.m: $1_def.r
	def_r2m $1; matlab_tidy $1_def.m;
$1_def.h: $1_def.m $1_sympar.txt
	def_m2h.sh $1

#SUMMARY rdae	raw differential-algebraic equations (r)
#SUMMARY dae	differential-algebraic equations (r)
#SUMMARY dae*	differential-algebraic equations (m)
#SUMMARY dae	differential-algebraic equations (tex)
#SUMMARY dae*	differential-algebraic equations (view)

# The main simulation programme
ifeq ($integration_method,implicit)
$1_ode2odes.m : $1_def.r $1_sympars.txt \
                $1_smxa.m $1_smxax.m\
                $1_simpar.m $1_numpar.m $1_state.m $1_input.m \
                $1_csex.m $1_cseo.m  $1_logic.m
ifeq ($using_oct,yes)












	touch $1_ode2odes.m # Create a dummy which wont' be used
	mtt $mtt_switches -q -u $1 ode2odes oct
else
	make_ode2odes $1 m $integration_method
endif
endif
ifneq ($integration_method,implicit)
$1_ode2odes.m : $1_def.r $1_sympars.txt\
		$1_simpar.m $1_numpar.m $1_state.m $1_input.m \
		$1_ode.m $1_odeo.m  $1_logic.m
ifeq ($using_oct,yes)









	echo "*** Warning: Shouldn't be here! Creating dummy $1_ode2odes.m"
	touch $1_ode2odes.m # Create a dummy which wont' be used
	mtt $mtt_switches -q -u $1 ode2odes oct
else
	make_ode2odes $1 m $integration_method
endif

endif

#SUMMARY ode2odes Simulation function (m)
#SUMMARY ode2odes Simulation function (cc)
#SUMMARY ode2odes Simulation function (oct)
#SUMMARY ode2odes Simulation function (exe)
$1_ode2odes.exe: $1_def.h $1_sympar.h\
		 $1_ode2odes.o $1_ode2odes_common.o $1_ode2odes_${integration_method}.o
	echo Creating $1_ode2odes.exe
	${MTT_CXX} ${MTT_CXXFLAGS} -o $1_ode2odes.exe\
	$1_ode2odes.o $1_ode2odes_common.o $1_ode2odes_${integration_method}.o\
	${MTT_LDFLAGS} ${MTT_CXXLIBS}

$1_ode2odes.o: $1_ode2odes.cc $1_ode2odes_common.o $1_ode2odes_${integration_method}.o
	echo Creating $1_ode2odes.o
	${MTT_CXX} ${MTT_CXXFLAGS} ${MTT_CXXINCS} -c $1_ode2odes.cc -DSTANDALONE

$1_ode2odes.oct: $1_ode2odes.cc $1_ode2odes_common.oct $1_ode2odes_${integration_method}.oct
	touch $1_ode2odes.m
	echo Creating $1_ode2odes.oct
	$MKOCTFILE $1_ode2odes.cc

$1_ode2odes.cc: $1_def.r $1_sympars.txt\
		$1_ode2odes_common.m $1_ode2odes_common.cc\
		$1_ode2odes_${integration_method}.m $1_ode2odes_${integration_method}.cc
	touch $1_ode2odes.m
	make_ode2odes $1 cc $integration_method

#Conversion of m to p to c
#SUMMARY ode	ordinary differential equations (c)
#SUMMARY ode	ordinary differential equations (p)
#SUMMARY state	state declaration (c) 
#SUMMARY state	state declaration (p) 
$1_simpar.p : $1_def.r $1_simpar.m

ifeq ($integration_method,implicit)
$1_ode2odes.p : $1_ode2odes.m $1_def.r $1_smxa.p $1_smxax.p\
		$1_simpar.p $1_numpar.p $1_state.p $1_input.p \
		$1_csex.p $1_cseo.p  $1_logic.p
	mtt_m2p $1_ode2odes.m $integration_method $stdin
endif
ifneq ($integration_method,implicit)
$1_ode2odes.p : $1_ode2odes.m $1_def.r\
		$1_simpar.p $1_numpar.p $1_state.p $1_input.p \
		$1_ode.p $1_odeo.p  $1_logic.p
	mtt_m2p $1_ode2odes.m $integration_method $stdin
endif

$1_ode2odes.c: $1_ode2odes.p 

	m2dat $1_sro
1_odes.c: $1_ode.c $1_def.r $1_sympar.r
	ode2odes_r2c $1 
$1_sms.c: $1_sm.c $1_def.r $1_sympar.r
	sm2sms_r2c $1

ifeq ($computation,octave)
ifeq ($using_oct,yes)
$1_odes.dat2: $1_ode2odes.oct
	time octave_ode2odes $1
else
$1_odes.dat2: $1_ode2odes.m
	time octave_ode2odes $1
endif
endif

ifeq ($computation,c)
$1_ode2odes.out: $1_ode2odes.c
	echo Creating $1_odes.out
	$CC $1_ode2odes.c -lp2c -lm -o $1_ode2odes.out

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.8  2000/12/27 13:11:42  peterg
## *** empty log message ***
##
## Revision 1.7  2000/08/04 15:48:30  peterg
## Added PDF
##
## Revision 1.6  2000/05/22 12:37:24  peterg
## *** empty log message ***
##
## Revision 1.5  2000/05/22 12:34:07  peterg

mtt_check_var "$FIG" "FIG" "$print" "MTT"
mtt_check_var "$PSVIEW" "PSVIEW" "$print" "MTT"
mtt_check_var "$PDFVIEW" "PDFVIEW" "$print" "MTT"
mtt_check_var "$HTMLVIEW" "HTMLVIEW" "$print" "MTT"
mtt_check_var "$DVIVIEW" "DVIVIEW" "$print" "MTT"
mtt_check_var "$CC" "CC" "$print" "MTT"
mtt_check_var "$LATEX2HTML" "LATEX2HTML" "$print" "MTT"
mtt_check_var "$MKOCTFILE" "MKOCTFILE" "$print" "MTT"
mtt_check_var "$MTT_CXX" "MTT_CXX" "$print" "MTT"
mtt_check_var "$MTT_CXXFLAGS" "MTT_CXXFLAGS" "$print" "MTT"
mtt_check_var "$MTT_CXXLIBS" "MTT_CXXLIBS" "$print" "MTT"
mtt_check_var "$MTT_CXXINCS" "MTT_CXXINCS" "$print" "MTT"
mtt_check_var "$MTT_LDFLAGS" "MTT_LDFLAGS" "$print" "MTT"

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.17.2.3  2001/03/01 05:05:53  geraint
## Minor revisions.
##
## Revision 1.17.2.2  2001/02/23 03:53:53  geraint
## Minor revisions. (ode2odes.exe works on GNU/Linux and ~GNU/Irix)
##
## Revision 1.17.2.1  2001/02/22 06:41:40  geraint
## Initial merge of ode2odes.exe into main mtt.
## standalone_rep.make deleted: rules moved to mtt, variables to mttrc.
##
## Revision 1.17  2000/12/27 16:46:13  peterg
## Stripped the mtt- from paths
##
## Revision 1.16  2000/12/27 15:16:44  peterg
## If then else format
##
## Revision 1.15  2000/12/27 14:57:43  peterg
## Now takes the base path as an argument
##
## Revision 1.14  2000/12/27 13:11:43  peterg

  export LATEX2HTML="latex2html -contents_in_navigation -index_in_navigation -address http://mtt.sourceforge.net"
  
  # Ascend stuff
  export ASCENDLIBRARY=$MTTPATH/ascend/lib
  
  # Oct file generation - use version with no optimisation.
  export MKOCTFILE=$MTT_LIB/octave/mkoctfile

  # ode2odes.exe stuff

    # local system

    PLAT="i686-pc-linux-gnu"
    PREFIX="/usr/local"
    GCCVERS="2.95.2"
    SRCOCTAVE="/cvs/octave"

#    PLAT="mips-sgi-irix6.5"
#    PREFIX="/usr/people/bevangp/GNU"
#    GCCVERS="2.95.2"
#    SRCOCTAVE="${PREFIX}/../build/octave-2.1.33"

    # include paths

    IOCTAVE="-I${PREFIX}/include/octave"

    # library paths

    LOCTAVE="-L${PREFIX}/lib/octave -loctave -lcruft -loctinterp"
    LKPATHSEA="-L${SRCOCTAVE}/kpathsea -lkpathsea"
    LREADLINE=" -L${SRCOCTAVE}/readline -lreadline"
    LSYSTEM="-ldl -lm -lncurses"
    LF2C="-L${PREFIX}/lib/gcc-lib/${PLAT}/${GCCVERS} -lg2c"

    # compiler options

    DEBUG="-g"
    OPTIM="-O3"
    FLAGS="-fno-rtti -fno-exceptions -fno-implicit-templates"

    # exported variables

    export MTT_CXX="g++"
    export MTT_CXXFLAGS="${DEBUG} ${OPTIM} ${FLAGS}"
    export MTT_CXXLIBS="${LOCTAVE} ${LKPATHSEA} ${LREADLINE} ${LF2C} ${LSYSTEM}"
    export MTT_CXXINCS="-I. ${IOCTAVE}"
    export MTT_LDFLAGS=" "
fi

# Copyright (c) P.J.Gawthrop 1991, 1992, 1994.

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.14  2000/12/28 12:24:35  peterg
## *** empty log message ***
##
## Revision 1.13  2000/10/11 08:52:46  peterg
## Creates csex (cse with dxe only) rep.
##
## Revision 1.12  2000/10/10 21:00:58  peterg
## New code genration
##
## Revision 1.11  1998/11/26 09:18:55  peterg

## Sorted out bug when MTTNz=0
##
## Revision 1.1  1996/08/15 16:47:02  peter
## Initial revision
##
###############################################################






#Explicit solution option
solve=0
while [ -n "`echo $1 | grep '^-'`" ]; do
  case $1 in
        -A )
                solve=1;;
        *)
                echo "$1 is an invalid argument - ignoring" ;;
  esac
  shift
done

# Create the reduce output code
def2write_r $1 cse      
def2write_r $1 csex  # Version without E matrix
def2write_r $1 cseo

if [ "$solve" = "1" ]; then
    echo "Creating $1_cse.r (with explicit solution of algebraic equations)"
else
    echo "Creating $1_cse.r"
fi

touch $1_cseo.r1
touch $1_cseo.r2
cat $1_cse.r1 $1_cse.r2  > $1_cse.r
cat $1_csex.r1 $1_csex.r2  > $1_csex.r
cat $1_cseo.r1 $1_cseo.r2  > $1_cseo.r

if [ "$solve" = "1" ]; then
    echo "Setting MTTNyz=0 in $1_def.r and updating other $1_def files"
    awk '{
     if ($1=="MTTNyz") 
       print "MTTNyz := 0;" 
     else print $0
    }' $1_def.r > mtt_junk
  # Make sure it preserves the time stamp!!
  # and remove dependent reps
   touch -r $1_def.r mtt_junk

   mv mtt_junk $1_def.r
   for file in `ls $1_def.*`; do
      if [ $file != $1_def.r ]; then
          ext=`echo $file|awk -F\. '{print $2}'`
          mtt -q $1 def $ext
      fi
   done
fi


# Now invoke the standard error handling.
mtt_error_r dae2cse_r.log

% %%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Field separator
if nargin<3
  FS = ';';
end;

arg = '';
if strcmp(args, '')==0
  L = length(args);
  args_count = 0;
  for i=1:n
    arg_count = 0;
    arg = '';
    if args_count == L
      break;
    end;  
    while args_count < L
      args_count = args_count+1;
      arg_count = arg_count+1;
      ch = str2ch(args,args_count);
      if ch==FS
	break;
      end;
      arg = [arg ch];
    end;
  end;
end;

#! /bin/sh

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

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.51.2.6  2001/03/16 03:56:13  geraint
## Removed psignal/siginfo.h - problematic and unnecessary.
##
## Revision 1.51.2.5  2001/03/12 23:16:37  geraint
## Minor improvements to signal handling (.exe).
##
## Revision 1.51.2.4  2001/03/12 03:59:30  geraint
## SIGINT  (C-c C-c) now causes simulation data to be dumped to MTT.core.
## SIGQUIT (C-c C-\) as for SIGINT, then raises default SIGQUIT.
## SIGFPE  as for SIGINT, then raises default SIGABRT.
##
## Revision 1.51.2.3  2001/03/07 04:06:55  geraint
## Irix: catch SIGFPE and write data before aborting (.exe).
## GNU/Linux: nada.
##
## Revision 1.51.2.2  2001/03/02 00:45:21  geraint
## Separated Euler and Implicit methods in .cc code and dependencies.
##
## Revision 1.51.2.1  2001/03/01 05:05:53  geraint
## Minor revisions.
##
## Revision 1.51  2001/02/19 06:33:19  geraint
## Removed operation form loop.
##
## Revision 1.50  2001/02/18 09:18:49  geraint
## Removed temporary Matrices from mtt_implicit.cc
##
## Revision 1.49  2001/02/14 06:06:34  geraint
## Removed octave_value_list wrappers from standalone.exe - speed improvements
##
## Revision 1.48  2001/02/11 07:08:59  geraint

  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
EOF

if [ "$method" = "rk4" ]; then
cat << EOF >> $filename
    [k1] = ddt * ${sys}_${ode}(x,u,t,par);
    [k2] = ddt * ${sys}_${ode}(x+k1/2,u,t+ddt/2,par);
    [k3] = ddt * ${sys}_${ode}(x+k2/2,u,t+ddt/2,par);
    [k4] = ddt * ${sys}_${ode}(x+k3,u,t+ddt,par);
    [dx] = [k1 + 2.0 * [k2 + k3] + k4] / (6.0 * ddt);
EOF
else
cat << EOF >> $filename
    [dx] = ${sys}_$ode(x,u,t,par);	# State derivative
EOF
fi

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

    vector_value=vector_value
    feval_header=toplev.h
fi

cat <<EOF  > $filename
#include <octave/oct.h>
#include <octave/ov-struct.h>
#include <octave/load-save.h>
#include <octave/lo-mappers.h>
#include <octave/variables.h>

#ifndef STANDALONE
#include <octave/${feval_header}>
#endif

#include "${sys}_def.h"
#include "${sys}_sympar.h"

#ifdef STANDALONE
#include <csignal>
#include <fstream>

extern ColumnVector F${sys}_input (
	ColumnVector &x,
	ColumnVector &y,
	const double &t,

	ColumnVector &par);

extern ColumnVector F${sys}_logic (
	ColumnVector &x,
	ColumnVector &u,
	const double &t,
	ColumnVector &par);

extern ColumnVector F${sys}_numpar (
	void);

extern Octave_map F${sys}_simpar (
     	void);

extern ColumnVector F${sys}_state (
	ColumnVector &x);

extern ColumnVector F${sys}_${ode} (
	ColumnVector &x,
	ColumnVector &u,
	const double &t,
	ColumnVector &par);

extern ColumnVector F${sys}_${odeo} (
	ColumnVector &x,
	ColumnVector &u,
	const double &t,
	ColumnVector &par);

EOF
if [ "$method" != "implicit" ]; then
cat <<EOF >> $filename
extern ColumnVector Fmtt_euler (
	ColumnVector &x,
	const ColumnVector &dx,
	const double &ddt,
	const int &nx,
	const ColumnVector &open_switches);
 
EOF
else
cat <<EOF >> $filename
extern ColumnVector Fmtt_implicit (
	ColumnVector &x,
	ColumnVector &dx,
	Matrix &AA,
	ColumnVector &AAx,

	const double &ddt,
	const int &nx,
	const ColumnVector &open_switches);

extern Matrix F${sys}_smxa (
	ColumnVector &x,
	ColumnVector &u,
	const double &t,
	ColumnVector &par);

extern ColumnVector F${sys}_smxax (
	ColumnVector &x,
	ColumnVector &u,
	const double &t,
	ColumnVector &par);
 
EOF
fi
cat <<EOF >> $filename
#endif // STANDALONE


inline ColumnVector
mtt_input (ColumnVector &x,
	   ColumnVector &y,
	   const double &t,
	   ColumnVector &par)
{
#ifdef STANDALONE
  return F${sys}_input (x, y, t, par);
#else
  static octave_value_list args, f;
  args (0) = octave_value (x);
  args (1) = octave_value (y);
  args (2) = octave_value (t);
  args (3) = octave_value (par);
  f = feval ("${sys}_input", args, 1);
  return f(0).${vector_value} ();
#endif
}

inline ColumnVector
mtt_logic (ColumnVector &x,
	   ColumnVector &u,
	   const double &t,
	   ColumnVector &par)
{
#ifdef STANDALONE
  return F${sys}_logic (x, u, t, par);
#else
  static 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}_logic", args, 1);
  return f(0).${vector_value} ();
#endif
}

inline ColumnVector
mtt_numpar (void)
{
#ifdef STANDALONE
  return F${sys}_numpar ();
#else
  static octave_value_list args, f;
  f = feval ("${sys}_numpar", args, 1);
  return f(0).${vector_value} ();
#endif
}

inline Octave_map
mtt_simpar (void)
{
#ifdef STANDALONE
  return F${sys}_simpar ();
#else
  static octave_value_list args;
  static 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);
#endif
}

inline ColumnVector
mtt_state (ColumnVector &x)
{
#ifdef STANDALONE
  return F${sys}_state (x);
#else
  static octave_value_list args, f;
  args (0) = octave_value (x);
  f = feval ("${sys}_state", args, 1);
  return f(0).${vector_value} ();
#endif
}

inline ColumnVector
mtt_${ode} (ColumnVector &x,
	    ColumnVector &u,
	    const double &t,
	    ColumnVector &par)
{

  args (2) = octave_value (t);
  args (3) = octave_value (par);
  f = feval ("${sys}_${odeo}", args, 1);
  return f(0).${vector_value} ();
#endif
}

EOF
if [ "$method" != "implicit" ];then
cat <<EOF >> $filename
inline ColumnVector
mtt_euler (ColumnVector &x,
	   const ColumnVector &dx,
	   const double &ddt,
	   const int &nx,
	   const ColumnVector &open_switches)
{
#ifdef STANDALONE
  return Fmtt_euler (x, dx, ddt, nx, open_switches);
#else
  static 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);
  f = feval ("mtt_euler", args, 1);
  return f(0).${vector_value} ();
#endif
}

EOF
else
cat <<EOF >> $filename
inline ColumnVector
mtt_implicit (ColumnVector &x,
	      ColumnVector &dx,
	      Matrix &AA,
	      ColumnVector &AAx,
	      const double &ddt,
	      const int &nx,

  args (5) = octave_value ((double)nx);
  args (6) = octave_value (open_switches);
  f = feval ("mtt_implicit", args, 1);
  return f(0).${vector_value} ();
#endif
}









































































inline Matrix
mtt_smxa (ColumnVector &x,
	  ColumnVector &u,
	  const double &t,
	  ColumnVector &par)
{
#ifdef STANDALONE

  args (2) = octave_value (t);
  args (3) = octave_value (par);
  f = feval ("${sys}_smxax", args, 1);
  return f(0).${vector_value} ();
#endif
}













EOF
















fi

cat <<EOF >> $filename

inline void
mtt_write (const double &t,
	   ColumnVector &x,
	   ColumnVector &y,
	   const int &nrows,
	   const bool dump_data = false,
	   ostream &file = cout)
{
  static Matrix data;
  static int row;

  if (dump_data)
    {
       Matrix written_data = data.extract (0, 0, row-1, data.cols ()-1);
       save_ascii_data_for_plotting (file, written_data, "MTT_data");
       return;
    }

  const int nx = x.length (), ny = y.length ();
  register int col = 0;

  if (0 == row)
    data = Matrix (nrows, 1+ny+1+nx, 0.0);

  data.elem (row, col) = t;
  for (register int i = 0; i < ny; i++)
    data.elem (row, ++col) = y.elem (i);
  data.elem (row, ++col) = t;
  for (register int i = 0; i < nx; i++)
    data.elem (row, ++col) = x.elem (i);

  row++;

  if (nrows == row)
    {
#ifdef STANDALONE
      save_ascii_data_for_plotting (file, data, "MTT_data");
//    cout << data << endl;
#else // ! STANDALONE
  set_global_value ("MTT_data", data);
#endif
    }
}

#ifdef STANDALONE
void dump_data (ostream &file)
{
  ColumnVector null (0.0);
  mtt_write (0.0, null, null, 0, true, file);
}

void set_signal_handlers (void);

void handle_signal (int signum)
{
  // handle some signals to ensure data is written.
  cerr << "# Writing data to MTT.core (signal " << signum << ")" << endl;
  ofstream corefile ("MTT.core");
  dump_data (corefile);
  switch (signum)
    {
      case SIGFPE:
        // Intel chips do not raise SIGFPE for DIVZERO :-(
        raise (SIGABRT);
        break;
      case SIGINT:
	break;
      case SIGQUIT:
	signal (SIGQUIT, SIG_DFL);
	raise (SIGQUIT);
	break;
      default:
        cerr << "# Warning: make_ode2odes needs updating!" << endl;
	signal (signum, SIG_DFL);
	raise (signum);
        break;
    }
  corefile.close ();
  set_signal_handlers ();
}

void set_signal_handlers (void)
{
  signal (SIGFPE,	handle_signal);
  signal (SIGINT,	handle_signal);
  signal (SIGQUIT,	handle_signal);
}

int main (void) {
  set_signal_handlers ();
#else
DEFUN_DLD (${sys}_ode2odes, args, ,
"Octave ode2odes representation of system with $method integration method
Usage: ${sys}_ode2odes (x, par, simpar)
")
{
  static octave_value_list retval;
#endif // STANDALONE
  static ColumnVector	x;
  static ColumnVector	par;

  static ColumnVector	open_switches (MTTNX);

  register double t	= 0.0;

  const double	ddt	= dt / stepfactor;
  const int	ilast	= static_cast<int> (round ((last - first) / ddt)) + 1;
  const int	nrows	= static_cast<int> (round ((last - first) / dt)) + 1;

  for (register int i = 0; i < MTTNU; i++)
    {
      u (i) = 0.0;
    }

  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, nrows);
	}
EOF
if [ "$method" = "rk4" ]; then
cat << EOF >> $filename
      {
        static ColumnVector
          k1 (MTTNX,0.0),
          k2 (MTTNX,0.0),
          k3 (MTTNX,0.0),
          k4 (MTTNX,0.0);

        const double
          t1 = t + ddt/2.0,
          t2 = t + ddt;

        ColumnVector
          x1 (x),
          x2 (x),
          x3 (x);

        k1 = ddt * mtt_${ode} (x , u, t , par); x1 += k1 * 0.5;
        k2 = ddt * mtt_${ode} (x1, u, t1, par); x2 += k2 * 0.5;
        k3 = ddt * mtt_${ode} (x2, u, t1, par); x3 += k3;
        k4 = ddt * mtt_${ode} (x3, u, t2, par);
        dx = (k1 + 2.0 * (k2 + k3) + k4) / (6.0 * ddt);
      }
EOF
else
cat << EOF >> $filename
      dx = mtt_${ode} (x, u, t, par);
EOF
fi

if [ "$method" = "implicit" ]; then
cat <<EOF >> $filename

      AA = mtt_smxa (x, u, ddt, par);
      AAx = mtt_smxax (x, u, ddt, par);
EOF

#! /bin/sh

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

# Bourne shell script: struc_txt2m

# Structure file - txt to m conversion
# P.J.Gawthrop April 1997
# Copyright (c) P.J.Gawthrop, 1997.

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
###############################################################



# Inform user
echo "Creating $1_struc.m"

rm -f mtt_error

#Write some file headers
echo "function [input_name,output_name,state_name] = $1_struc" > $1_struc.m
echo "%% Structure file ($1_struc.m)" >> $1_struc.m
echo "%% Generated by MTT at `date`" >> $1_struc.m

# This is the main transformation using awk
awk '
function header(what){
  print what "_name = [";
}

function footer(what){
  print "];";
}

BEGIN{
Which="";
apostrophe = "\047";
}
{
  if ($1!=Which) {
    if (Which!="") footer(Which);
    header($1);
  }
  print  apostrophe $4 "_" $3 apostrophe;
  Which = $1;
}
END{
  footer(Which)
}
'  SYSTEM=$1  < $1_struc.txt >> $1_struc.m 2>mtt_error.txt

# Now invoke the standard error handling.
mtt_error mtt_error.txt

	mkoctfile $<
endif



all:	$(sys)_$(rep).cc

parser.exe: parse_m2cc.cc
	g++ $< $(CFLAGS) -o ${MTTPATH}/trans/parser.exe -I${MTTPATH}/../cc/include

ifeq ($(rep),ode2odes)
$(sys)_ode2odes.cc: ode2odes.cc
	cat ode2odes.cc | sed 's/\$$/$(sys)/' > $(sys)_ode2odes.cc
else
$(sys)_$(rep).cc: $(sys)_def.h $(sys)_sympar.h $(sys)_def.r $(sys)_$(rep).m
	mtt_m2cc.sh ${sys} ${rep} parser

#ifndef HAVE_USEFUL_FUNCTIONS_HH
#define HAVE_USEFUL_FUNCTIONS_HH


#ifdef __CPLUSPLUS
template <class class_t>
#else
#define inline			// strip from pre-processed file
#define class_t double
#endif // __CPLUSPLUS


// == Template functions ==

static inline class_t
max (const class_t &x1, const class_t &x2)
{
  return ((x1 >= x2) ? x1 : (x1 < x2) ? x2 : 0);
}

static inline class_t
min (const class_t &x1, const class_t &x2)
{
  return ((x1 <= x2) ? x1 : (x1 > x2) ? x2 : 0);
}

static inline class_t
sign (const class_t &x)
{
  return ((x > 0) ? +1 : (x < 0) ? -1 : 0);
}


// == Octave functions ==

static inline Matrix
ones (const int r = 1, const int c = 1)
{
  Matrix m (r, c, 1.0);
  return m;
}

static inline ColumnVector
nozeros (const ColumnVector v0, const double tol = 0.0)
{
  ColumnVector v (v0.length ());
  register int i, j;
  for (i = j = 0; i < v.length (); i++)
    if (tol < abs (v0 (i)))
      {
	v (j) = v0 (i);
	j++;
      }
  if (0 == j)
    {
      return *new ColumnVector ();
    }
  else
    {
      return (v.extract (0, --j));
    }
}

static inline ColumnVector
zeros (const int r)
{
  ColumnVector v (r, 0.0);
  return v;
}

static inline Matrix
zeros (const int r, const int c)
{
  Matrix m (r, c, 0.0);
  return m;
}

#endif // HAVE_USEFUL_FUNCTIONS_HH

/* $Id$
 * $Log$
 * Revision 1.1.2.1  2001/03/09 02:59:26  geraint
 * got_comment: (char)c no longer compared to (int)EOF.
 *
 * Revision 1.1  2000/12/28 09:46:05  peterg
 * put under RCS
 *
 * Revision 1.1  2000/10/31 04:29:50  geraint
 * Initial revision
 *
 */

{
  /*
   * get remainder of line in case there are any commented keywords
   */
  char c;
  cout << " // ";
  cin >> c;
  do {

    cout << c;
  } while (c != '\n' && cin >> c);

  cout << endl
       << indent () << ';' << endl
       << indent ();
}
  

		/*
		 * EOL is end-of-statement in Octave, add ;
		 */
		   << indent () << ';' << endl
		   << indent ();
	    }
	  buf = "";
	  break;
	case '\\':
	  cin >> c;
	  if ('\n' == c)
	    {
	      buf += '\n';
	    }
	  else
	    {
	      buf += '\\';
	      buf += c;
	    }
	  break;
	default:
	  buf += c;
	}
      if (lbrace->get_nesting_depth ())
	{
	  cout.setf (ios::skipws);

#! /bin/sh
# $Id$
# $Log$
# Revision 1.2  2001/02/05 13:03:19  geraint
# Restrict scope of variables to file (static).
# Warn GCC that variables may be unused.
#
# Revision 1.7  2001/01/09 15:43:50  geraint
# Warn gcc that variables may be unused.
#
# Revision 1.6  2001/01/08 05:47:56  geraint
# Restrict scope of variables to file (static)
#
# Revision 1.5  2000/12/05 12:44:55  peterg

TMP_VAR_NAMES="mtt_tmp mtt_o $*"

IN=${SYS}_sympar.txt
OUT=${SYS}_sympar.h

declare_sys_param ()
{
cat ${IN} | awk '{printf ("static double %s MTT_UNUSED;\t// %s\n", tolower($1), $2)}'
}

declare_temp_vars ()
{
for name in ${TMP_VAR_NAMES}
do
    echo ""

{
#else // !STANDALONE
DEFUN_DLD (mtt_implicit, args, ,
	   "implicit integration method")
{
#ifdef OCTAVE_DEV
  static ColumnVector  	x	= args(0).column_vector_value ();
  static ColumnVector	dx	= args(1).column_vector_value ();
  static Matrix		AA	= args(2).matrix_value ();
  static ColumnVector	AAx	= args(3).column_vector_value ();
  const  double		t	= args(4).double_value ();
  const  int		Nx	= (int) (args(5).double_value ());
  const  ColumnVector	openx	= args(6).column_vector_value ();
#else // !OCTAVE_DEV
  static ColumnVector	x	= args(0).vector_value ();
  static ColumnVector	dx	= args(1).vector_value ();
  static Matrix		AA	= args(2).matrix_value ();
  static ColumnVector	AAx	= args(3).vector_value ();
  const  double		t	= args(4).double_value ();
  const  int		Nx	= (int) (args(5).double_value ());
  const  ColumnVector	openx	= args(6).vector_value ();
#endif // OCTAVE_DEV
#endif // STANDALONE

  register int row, col;

  for (row = 0; row < Nx; row++)
    {
      if (openx (row) > 0.5)
	{
	  AAx (row) = 0.0;
          dx (row) = 0.0;
	  for (col = 0; col < Nx; col++)
	    {
	      AA (row,col) = 0.0;
	      AA (col,row) = 0.0;
	    }
	}
    }

  x = static_cast<ColumnVector> (xleftdiv (AA, static_cast<Matrix>(AAx + dx * t)));
  
  for (row = 0; row < Nx; row++)
    {
      if (openx (row) > 0.5)
	{
	  x (row) = 0.0;
	}
    }

#ifdef STANDALONE
  return x;