#! /bin/sh

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

# Bourne shell script: mtt

# P.J.Gawthrop Oct 1989, Dec 1989, 1990, 1991, Oct 1993, Dec 1993, Jan 1994, Dec 95
# Copyright (c) P.J.Gawthrop, 1989, 1990, 1991, 1993, 1994, 1995, 1996.

###############################################################
## Version control history
###############################################################
## $Header$
## $Log$
## Revision 1.20  1996/08/20 08:24:53  peter
## Now with version control.
##
## Revision 1.19  1996/08/19 15:33:22  peter
## Included new state matric (sm) rep.
## Generate step response (sr) from sm.
##
## Revision 1.18  1996/08/19 13:30:54  peter
## New CLEAN target.
##
## Revision 1.17  1996/08/19 06:48:14  peter
## Explicit figure editing.
##
## Revision 1.16  1996/08/18 19:58:36  peter
## Report generation included.
##
## Revision 1.15  1996/08/18 12:01:01  peter
## Unified format of time responses.
##
## Revision 1.14  1996/08/16 14:52:24  peter
## Added two alternative view options - dview and pview.
##
## Revision 1.13  1996/08/16 08:24:31  peter
## Correctly chooses how to make a ps file using $REPTYPE.
##
## Revision 1.12  1996/08/15 16:47:47  peter
## DAE solution (experimental) included.
##
## Revision 1.11  1996/08/15 11:55:30  peter
## Checks for changed argument.
## Handles frequency response.
##
## Revision 1.10  1996/08/15 07:44:24  peter
## Now handles generic transformations using %.
##
## Revision 1.9  1996/08/14 09:16:08  peter
## Step response now goes through  m, dat and gdat conversions
##
## Revision 1.8  1996/08/12 20:19:52  peter
## Arguments now passed via $1_args.m file
##
## Revision 1.7  1996/08/11 19:49:27  peter
## Parameter passing (via $4) added.
##
## Revision 1.6  1996/08/11 10:40:56  peter
## Added new numpar stuff.
## Added step response.
##
## Revision 1.5  1996/08/10 14:13:47  peter
## Added impulse and frequency response stuff.
##
## Revision 1.4  1996/08/10 09:19:48  peter
## Put in help, info and warranty stuff.
##
## Revision 1.3  1996/08/09 15:13:13  peter
## Generic viewing of any representation via ps files implemented.
##
## Revision 1.2  1996/08/05 19:50:55  peter
## Put in fig.fig target.
##
## Revision 1.1  1996/08/04 17:29:42  peter
## Initial revision
##
###############################################################
echo
echo 'MTT (Model Transformation Tools) version 2.0'
echo 'This is free software with ABSOLUTELY NO WARRANTY.'
echo 'Type `mtt warranty'\' 'for details.'
echo

if [ "$1" = "" ]; then
  echo 'Usage: mtt help'
  echo '       mtt info'
  echo '       mtt info topic'
  echo '       mtt manual'
  echo '       mtt warranty'
  echo '       mtt clean'
  echo '       mtt system_name clean'
  echo '       mtt system_name representation vc'
  echo '       mtt system_name representation language'
  echo '       mtt system_name representation language parameters'
  exit
fi


if [ "$1" = "info" ]  && [ "$3" = "" ]; then
  echo Invoking info
  info -f $MTTPATH/../doc/mtt.info $2
  exit
fi

if [ "$1" = "hinfo" ] && [$2 = ""]; then
   echo Invoking $HTMLVIEW
   (cd $MTTPATH/../doc; $HTMLVIEW ./mtt.html)&
   exit
fi

if [ "$1" = "manual" ] && [$2 = ""]; then
  echo Invoking $DVIVIEW
  $DVIVIEW  $MTTPATH/../doc/mtt.dvi &
  exit
fi

if [ "$1" = "warranty" ] && [$2 = ""]; then
  cat << EOF
Copyright (C) 1989, 1990, 1991, 1992, 1993, 1994, 1995,1996 Peter J. Gawthrop
This is free software with ABSOLUTELY NO WARRANTY.

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program. If not, write to the Free Software
Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.

EOF
  exit
fi

# Clean up 
if [ "$1" = "clean" ] && [ "$2" = "" ]; then
  echo 'Removing all generated files for all systems'
  rm -f *.log  mtt_info.txt warning.txt
  rm -f *_abg.m *_args.* 
  rm -f *_rbg.* *_cmp.* *_fig.fig *_*cbg.* *_ese.* *_def.*
  rm -f *_dae*.* *_cse.* *_ode*.* *_obs.* *_rfe.* *_ss.*
  rm -f *_dm.* *_csm.* *_sm.* *_tf.* *_sr*.* *_*fr.*
  rm -f *_rep.*
  exit
fi

# Clean up named system
if [ "$2" = "clean" ] && [ "$3" = "" ]; then
  echo 'Removing all generated files for system ' $1
  rm -f *.log  mtt_info.txt warning.txt
  rm -f $1_abg.m $1_args.* 
  rm -f $1_rbg.* $1_cmp.* $1_fig.fig $1_*cbg.* $1_ese.* $1_def.*
  rm -f $1_dae*.* $1_cse.* $1_ode*.* $1_obs.* $1_rfe.* $1_ss.*
  rm -f $1_dm.* $1_csm.* $1_sm.* $1_tf.* $1_sr*.* $1_*fr.*
  rm -f $1_rep.*
  exit
fi

# Invoke explicit requests for modification
if [ "$2" = "abg" ] && [ "$3" = "fig" ]; then
  echo Editing $1_$2.$3
  $FIG $1_$2.$3 &
  exit
fi


# Save up the argument list in a file; but only if argument has changed
ARGS=$4;
cat > $1_args.new <<EOF
$ARGS
EOF

DIFF=$(diff -bq $1_args.m $1_args.new 2>/dev/null)

if [ -n "$DIFF" ]; then
  mv $1_args.new $1_args.m
fi

# Classify the representation - needed for conversion route to postscript
case $2 in
	abg )
		REPTYPE='bg' ;;
	cbg )
		REPTYPE='bg' ;;
	odes)
		REPTYPE='data' ;;
	odeso)
		REPTYPE='data' ;;
	odess)
		REPTYPE='data' ;;
	odesso)
		REPTYPE='data' ;;
	daes)
		REPTYPE='data' ;;
	daeso)
		REPTYPE='data' ;;
	sr)
		REPTYPE='data' ;;
	sro)
		REPTYPE='data' ;;
	lmfr)
		REPTYPE='data' ;;
	lpfr)
		REPTYPE='data' ;;
	nyfr)
		REPTYPE='data' ;;
	nifr)
		REPTYPE='data' ;;
	*)
		REPTYPE='tex' ;;
esac


# Version control
VC=''

if [ "$3" = "vc" ]; then
  case $2 in
	abg )
		VC='ok'
		VCext='fig' ;;
	lbl )
		VC='ok'
		VCext='txt' ;;
	 cr)
		VC='ok'
		VCext='r' ;;
	sympar )
		VC='ok'
		VCext='r' ;;
	simp )
		VC='ok'
		VCext='r' ;;
	numpar )
		VC='ok'
		VCext='m' ;;
	input )
		VC='ok'
		VCext='m' ;;
	sspar )
		VC='ok'
		VCext='r' ;;
	*)
		echo Version control is not appropriate for representation $2
		exit  ;;
  esac
fi


# This is the main mtt programme
 
make -s  -f -  $1_$2.$3 << EOF

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

# Create an arg file if it doesn't exist
$1_args.m:
	touch $1_args.m

#Create acausal bond graph figure using xfig
$1_abg.fig:
	$FIG $1_abg.fig 

#Create empty files (with titles) if not already there
$1_lbl.txt:
	echo Creating $1_lbl.txt
	( \
	echo "%% Label file for system $1 ($1_lbl.txt)"; \
	cat $MTTPATH/m/rcs_header.txt; \
	echo "%% Each line should be of one of the following forms:"; \
	echo "%	a comment (ie starting with %)"; \
	echo "%	Component-name	CR_name	arg1,arg2,..argn"; \
	echo "%	blank"; \
        )> $1_lbl.txt
$1_cr.r:
	makecr $1
$1_sympar.r:
	echo Creating $1_sympar.r
	( \
	echo "%% Symbolic parameter file for system $1 ($1_sympar.r)"; \
	cat $MTTPATH/m/rcs_header.txt; \
	echo "%% Change the following line to reflect the number of variables"; \
	echo "MTTNVar := 0;"; \
	echo "IF MTTNVar>0 THEN MATRIX MTTVar(MTTNVar,1);"; \
 	echo "%% Add an assignment statement here for each variable"; \
 	echo "%% For example: MTTVar(1,1) := r"; \
        )> $1_sympar.r
$1_simp.r:
	echo Creating $1_simp.r
	( \
	echo "%% Reduce comands to simplify output for system $1 ($1_simp.r)"; \
	cat $MTTPATH/m/rcs_header.txt; \
        )> $1_simp.r

$1_numpar.m: $1_sympar.r
	mknumpar $1; matlab_tidy $1_numpar.m;
	mv $1_numpar.m junk; cat $MTTPATH/m/rcs_header.txt junk>$1_numpar.m

$1_input.m:
	echo Creating $1_input.m
	( \
	echo "function u = $1_input(t)"; \
	echo "%% Input file for system $1 ($1_input.m)"; \
	echo "%% This file provides the system input for simulation:";\
	cat $MTTPATH/m/rcs_header.txt; \
	echo; \
	echo "% Set up system definitions"; \
	echo "[nx,ny,nu,nz,nyz] = $1_def;"; \
	echo "% Change the following lines as you wish"; \
	echo "u = zeros(nu,1);"; \
	echo "u(1) = 1;"; \
        )> $1_input.m

$1_sspar.r:
	echo Creating $1_sspar.r
	(\
	echo "%% Reduce steady-state parameter file for $1 ($1_sspar.r)" \
	cat $MTTPATH/m/rcs_header.txt; \
	)> $1_sspar.r

$1_rep.txt:
	echo Creating $1_rep.txt
	( \
	echo "%% Outline report file for system $1 ($1_rep.txt)"; \
	cat $MTTPATH/m/rcs_header.txt; \
	echo "abg ps"; \
	echo "lbl txt"; \
	echo "cbg ps"; \
	echo "dae tex"; \
	)> $1_rep.txt


#Raw bond graph: fig file to mfile
$1_rbg.m: $1_abg.fig $1_lbl.txt
	rbg_fig2m $1
$1_cmp.m: $1_rbg.m
$1_fig.fig: $1_rbg.m

#Bond graph conversions using matlab
#Raw bond graph to acausal bond graph: mfile
$1_abg.m: $1_rbg.m
	rbg2abg_m $1	

#Causal bond graph: mfile to fig conversion
$1_cbg.fig: $1_cbg.m $1_rbg.m $1_fig.fig
	cbg_m2fig $1


#Acausal bond graph to causal bond graph: mfile
$1_cbg.m: $1_abg.m
	abg2cbg_m $1

#Elementary system equations + definitions
$1_ese.r: $1_cbg.m
	   cbg2ese_m2r $1
$1_def.r: $1_ese.r
	touch $1_def.r
$1_def.m: $1_def.r
	def_r2m $1; matlab_tidy $1_def.m;


#Differential-algebraic equations
$1_dae.r:	$1_ese.r $1_cr.r $1_sympar.r
	ese2dae_r $1; tidy $1_dae.r
$1_dae.m: $1_def.r $1_dae.r  $1_sympar.r
	dae_r2m $1; matlab_tidy $1_dae.m; matlab_tidy $1_daeo.m
$1_dae.c: $1_def.r $1_dae.r  $1_sympar.r
	dae_r2c $1; c_tidy $1_dae.c
$1_dae.tex:	$1_dae.r $1_simp.r
	dae_r2tex $1; latex_tidy $1_dae.tex

#Constrained-state equations
$1_cse.r: $1_dae.r $1_def.r $1_sympar.r
	dae2cse_r $1; tidy $1_cse.r
$1_cse.m: $1_def.r $1_cse.r $1_sympar.r
	cse_r2m $1; matlab_tidy $1_cse.m
$1_cse.tex: $1_cse.r  $1_sympar.r
	cse_r2tex $1; latex_tidy $1_cse.tex

#Ordinary differential equations
$1_ode.r: $1_cse.r $1_def.r $1_sympar.r
	cse2ode_r $1; tidy $1_ode.r
$1_ode.m: $1_def.r $1_ode.r $1_sympar.r
	ode_r2m $1; matlab_tidy $1_ode.m; matlab_tidy $1_odeo.m
$1_ode.tex: $1_ode.r  $1_sympar.r $1_simp.r
	ode_r2tex $1; latex_tidy $1_ode.tex

#Observer functions for GPC
$1_obs.r: $1_ode.r $1_def.r $1_sympar.r
	ode2obs_r $1; tidy $1_obs.r
$1_obs.m: $1_def.r $1_obs.r $1_sympar.r
	obs_r2m $1; matlab_tidy $1_obs.m
$1_obs.tex: $1_obs.r  $1_sympar.r $1_simp.r
	obs_r2tex $1; latex_tidy $1_obs.tex


#Robot-form equations
$1_rfe.r: $1_cse.r $1_csm.r $1_def.r $1_sympar.r
	cse2rfe_r $1; tidy $1_rfe.r
$1_rfe.tex: $1_rfe.r  $1_sympar.r
	rfe_r2tex $1; latex_tidy $1_rfe.tex

#Steady-states
$1_ss.r: $1_def.r $1_sspar.r $1_sympar.r $1_cr.r $1_dae.r
	sspar2ss_r $1; tidy $1_ss.r
$1_ss.m: $1_def.r $1_sympar.r $1_ss.r 
	ss_r2m $1; matlab_tidy $1_ss.m
$1_ss.tex: $1_def.r $1_ss.r  $1_sympar.r $1_simp.r
	ss_r2tex $1; latex_tidy $1_ss.tex 


#Linearised system: descriptor matrices
$1_dm.r: $1_dae.r $1_def.r $1_cr.r  $1_ss.r $1_sympar.r
	dae2dm_r $1; tidy $1_dm.r
$1_dm.m: $1_def.r $1_dm.r $1_sympar.r
		dm_r2m $1; matlab_tidy $1_dm.m
$1_dm.tex: $1_dm.r $1_sympar.r $1_simp.r 
	dm_r2tex $1; latex_tidy $1_dm.tex


#Linearised system: constrained-state matrices
$1_csm.r: $1_cse.r $1_def.r $1_cr.r $1_ss.r $1_sympar.r
	cse2csm_r $1; tidy $1_csm.r
$1_csm.m: $1_def.r $1_csm.r
		csm_r2m $1;  matlab_tidy $1_csm.m
$1_csm.tex: $1_csm.r  $1_sympar.r $1_simp.r
	csm_r2tex $1; latex_tidy $1_csm.tex

#Linearised system: state matrices
$1_sm.r: $1_csm.r $1_def.r $1_cr.r $1_ss.r $1_sympar.r
	csm2sm_r $1; tidy $1_sm.r
$1_sm.m: $1_def.r $1_sm.r
		sm_r2m $1;  matlab_tidy $1_sm.m
$1_sm.tex: $1_sm.r  $1_sympar.r $1_simp.r
	sm_r2tex $1; latex_tidy $1_sm.tex

#Linearised system: transfer function
$1_tf.r: $1_dm.r $1_sympar.r 
	dm2tf_r $1; tidy $1_tf.r
$1_tf.m: $1_def.r $1_tf.r
	tf_r2m $1; matlab_tidy $1_tf.m
$1_tf.tex: $1_def.r $1_tf.r $1_sympar.r $1_simp.r
	tf_r2tex $1; latex_tidy $1_tf.tex


#Time responses
$1_sr.m: $1_sm.m $1_numpar.m $1_args.m $1_def.m
	sm2sr_m $1 '$ARGS'
$1_sro.m: $1_sr.m
	touch $1_sro.m
$1_odes.m: $1_ode.m $1_numpar.m $1_args.m $1_def.m $1_input.m
	ode2odes_m $1 '$ARGS'
$1_odeso.m: $1_odes.m
	touch $1_odeso.m
$1_daes.m: $1_dae.m $1_numpar.m $1_args.m $1_def.m $1_input.m
	dae2daes_m $1 '$ARGS'
$1_daeso.m: $1_daes.m
	touch $1_daeso.m

#Numerical steady states
$1_odess.m: $1_ode.m $1_numpar.m $1_args.m $1_def.m $1_input.m
	ode2odess_m $1 '$ARGS'

$1_odesso.m: $1_odess.m
	touch $1_odesso.m


#Frequency responses
$1_fr.m: $1_dm.m $1_numpar.m $1_args.m $1_def.m
	dm2fr_m $1 '$ARGS'
$1_lmfr.m: $1_fr.m
	touch $1_lmfr.m
$1_lpfr.m: $1_fr.m
	touch $1_lpfr.m
$1_nyfr.m: $1_fr.m
	touch $1_nyfr.m
$1_nifr.m: $1_fr.m
	touch $1_nifr.m

#Partially-known system indentification structure matrices
$1_pkim.r: $1_tf.r $1_sympar.r 
	tf2pkim_r $1; tidy $1_pki.r
$1_pkim.tex: $1_def.r $1_pkim.r $1_sympar.r
	pkim_r2tex $1; latex_tidy $1_pkim.tex

#Generic conversion of data files from m to dat format
$1_$2.dat: $1_$2.m
	m2dat $1_$2

#Generic conversion of data files from dat to (gplot) gdat format
$1_$2.gdat: $1_$2.dat
	dat2gdat $1_$2

#Generic conversion of Latex to latex document
$1_$2.doc: $1_$2.tex
	tex2doc $1 $2

#Generic conversion of LaTeX doc to dvi
$1_$2.dvi: $1_$2.doc
	echo Creating $1_$2.dvi
	latex $1_$2.doc > /dev/null

#Create PostScript version of fig files
ifeq ($REPTYPE,bg)
$1_$2.ps : $1_$2.fig
	echo Creating $1_$2.ps
	fig2dev -Lps $1_$2.fig >$1_$2.ps
endif

ifeq ($REPTYPE,data)
#Create PostScript version of dat file
$1_$2.ps: $1_$2.gdat
	gdat2ps $1_$2 '$ARGS'
endif

ifeq ($REPTYPE,tex)
#Create PostScript version of dvi file
$1_$2.ps: $1_$2.dvi
	echo Creating $1_$2.ps
	dvips -o $1_$2.ps $1_$2
endif

#Default view
ifeq ($REPTYPE,bg)
$1_$2.view : $1_$2.pview
endif

ifeq ($REPTYPE,data)
$1_$2.view : $1_$2.dview
endif

ifeq ($REPTYPE,tex)
$1_$2.view : $1_$2.tview
endif

#View a ps file
$1_$2.pview: $1_$2.ps
	echo Creating view of $1_$2
	$PSVIEW $1_$2.ps&

#View a gdat file
$1_$2.dview: $1_$2.gdat
	echo Creating view of $1_$2
	gdat2view $1_$2 '$ARGS'&
	
#View a tex file
$1_$2.tview: $1_$2.dvi
	echo Creating view of $1_$2
	$DVIVIEW $1_$2.dvi&
	

#Report generation
$1_rep.make: $1_rep.txt
	rep_txt2make $1
$1_rep.tex: $1_rep.make $(cat $1_rep.make 2>/dev/null)
	rep_txt2tex $1;
#$1_rep.make: $1_rep.tex
#	touch $1_rep.make
#$1_rep.rep: $1_rep.make
#	sh<$1_rep.make; cp $1_rep.tex $1_rep.doc

$1_rep.html: $1_rep.dvi
	echo Creating $1_rep.html 
	latex2html $1_rep.tex

# Version control system - uses RCS
RCS:
	echo Creating RCS
	mkdir RCS
$1_$2.vc: RCS $1_$2.$VCext 
	echo Version control on $1_$2.$VCext
	ci -I $1_$2.$VCext </dev/tty  ; co -l $1_$2.$VCext 

EOF

#! /bin/sh

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

# Bourne shell script: abg2cbg_m
#
# Acausal bond graph to causal bond graph: mfile format

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.7  1996/08/25 09:22:55  peter
## # Print errors, if any.
## mtt_error error.txt
##
## Revision 1.6  1996/08/24 14:36:01  peter
## Error handling included.
##
## Revision 1.5  1996/08/19 15:38:31  peter
## Removed bug work round.
##
## Revision 1.4  1996/08/16 14:28:45  peter
## Some debugging lines removed.
##
## Revision 1.3  1996/08/08 18:06:55  peter
## Unified naming scheme.
##
## Revision 1.2  1996/08/05 11:25:00  peter
## Removed bonds argument from abg2cbg.
## Removed call to $1_abg.
##
## Revision 1.1  1996/08/04 17:45:11  peter
## Initial revision
##
###############################################################

# The causal BG is described by four matrices
# (the same structure as the acausal BG except that Causality is 1 ot -1):
#   junctions: describes the junctions -
#     Column 1   Index of component (see components)
#     Column 2.. Indices of bonds on junction (see jbonds)
#                  rows may be padded with zeros
#
#    cjbonds:    describes the true bonds -
#     Column 1     Index of component (see components)
#     Column 2     Direction
#                    1 away from junction
#                   -1 towards junction
#     Column 3    Effort causality 
#                     1 stroke at arrow end
#                    -1 stroke at other end
#                     0 no stroke (ie acausal)
#     Column 4    Flow causality 
#                     1 stroke at arrow end
#                    -1 stroke at other end
#                     0 no stroke (ie acausal)
#     Columns 5:6  x,y coordinates of non-arrow end (in Fig)
#     Columns 7:8  x,y coordinates of arrow end (in Fig)
#
#   mbonds:    describes the modulation bonds -
#     Column 1     Index of component at non-arrow end
#     Column 2     Index of component at arrow end
#     Columns 4:5  x,y coordinates of non-arrow end (in Fig)
#     Columns 6:7  x,y coordinates of arrow end (in Fig)
#
#   components: describes the comonents -
#     Column 1     Type of component
#                    0  0-junction
#                    1  1-junction
#                    2  R component
#                    3  C component
#                    4  I component
#                    11 TF component
#                    12 GY component
#                    13 A  component
#     Column 2     Flag for explicit junction
#                    4  Explicit component (4 means text type)
#                    0  Implicit component
#     Column 3-14    fig file information (see fig documentation)
#
# Lists of under- and over-causal junctions are displayed.
#
# P.J.Gawthrop May 1996
# Copyright (c) P.J.Gawthrop, 1996.


# P.J.Gawthrop May 1996
# Copyright (c) P.J.Gawthrop, 1996.

infofile='mtt_info.txt';

# Remove the old log file
rm -f abg2cbg_m.log
rm -f $1_cbg.m
rm -f mtt_info.txt


#Inform user
echo Creating $1_cbg.m

# Use matrix manipulation to accomplish the transformation
$MATRIX  >abg2cbg_m.log  2>mtt_error.txt << EOF

  infofile= '$infofile';

  %Convert from acausal to causal bond graph in m-file form.
  system_name = '$1';

  % THIS IS A NASTY BUG WORKROUND FOR OCTAVE
  % abg2cbg('RC','RC' ,'', [], infofile);
  
  port_bonds = [];
  [cbonds,status] = abg2cbg(system_name, '', '', port_bonds,infofile);

EOF

cat mtt_info.txt

cp $1_$1_cbg.m $1_cbg.m 2>> mtt_error.txt 

# Print errors, if any.
mtt_error mtt_error.txt

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

# gawk script: lbl2sympar.awk
# Label file to symbolic parameters conversion
# P.J.Gawthrop August 1996
# Copyright (c) P.J.Gawthrop, 1996.

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.1  1996/08/24 13:34:48  peter
## Initial revision
##
###############################################################



function exact_match(name1, name2) {
  return ((match(name1,name2)>0)&&(length(name1)==length(name2)))
    }

function matches(name, names) {
  n_matches = split(names,match_name);
  matched = 0;
  for (i_matches = 1; i_matches <= n_matches; i_matches++) {
    if ( exact_match(name,match_name[i_matches]) ) {
      matched = 1;
      break;
    }
  }
  return matched;
    }


BEGIN {
comment = "%";
arg_delimiter = ",";
not_an_arg = "effort flow internal external zero";
numeric = "[0-9]";
symbol_count = 0;
symbols = "";
}
{
  if ( (match($1,comment)==0) && (NF>=3) ) {
    args = $3;
    n_args = split(args,arg,arg_delimiter);
    for (i = 1; i <= n_args; i++) {
      first_char = substr(arg[i],1,1);
      if ( (matches(not_an_arg,arg[i])==0) \
	   && (match(first_char,numeric)==0) \
	   && (length(arg[i])>0) \
	   && (matches(symbols,arg[i]) ==0) ) {
	symbol_count++;
	symbols = sprintf("%s %s", symbols, arg[i]);
	  }
    }
  }
}
END {
# print the _sympar file
  printf("%% Symbolic parameter file - generated by MTT\n\n");
  printf("MTTNVar := %1.0f;\n", symbol_count);

  if (symbol_count>0) {
    printf("MATRIX MTTVar(MTTNVar,1);\n");
    split(symbols,symbol);
    for (i = 1; i <= symbol_count; i++) {
      printf("MTTVar(%1.0f,1) := %s;\n", i, symbol[i]);
    }
  }
  printf("END;\n\n");  
}

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

# gawk script: rbg_fig2m.awk
# Raw  bond-graph conversion from fig to matlab
# P.J.Gawthrop June 1996
# Copyright (c) P.J.Gawthrop, 1996.

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.11  1996/08/19 10:48:57  peter
## Added `-' to the component regexp.
##
## Revision 1.10  1996/08/19 09:03:13  peter
## Parses repetative components: ie suffixed by *n.
##
## Revision 1.9  1996/08/09 08:23:11  peter
## Fixed bug: ports not recognised.
##
## Revision 1.8  1996/08/05 20:12:43  peter
## Now writes a _fig.fig file.
##
## Revision 1.7  1996/08/05 18:44:56  peter
## Now writes out a _cbg file without ----- symbol.
##
## Revision 1.6  1996/08/05 12:17:37  peter
## n_ports now appear in the _abg file instead.
##
## Revision 1.5  1996/08/05 12:01:28  peter
## The _cmp function now returns the number of ports.
##
## Revision 1.4  1996/08/05 10:14:46  peter
## Made ports appear, in order, at top of component lists
##
## Revision 1.3  1996/08/04 20:32:28  peter
## Stopped complaint about missing lbl entry for port components
##
## Revision 1.2  1996/08/04 20:05:25  peter
## Included port components - eg SS:[1]
##
## Revision 1.1  1996/08/04 20:01:58  peter
## Initial revision
##
###############################################################


##############################################################
# This (g)awk script reads a fig file in fig 3.1 format.
# It interprets the picture as: bonds, arrows and components
# as follows:
#
# Bonds are firm (not dashed etc) polylines with 2 line segments - 
#  fig represents this by  a firstline record where
#    field 1 = 2 (always 2)
#    field 2 = 1 (polyline)
#    field 3 = 0 (style is a firm line)
#    field 14 = 0 (no forward arrow)
#    field 15 = 0 (backward arrow) 
#  a data field starting with a tab followed by 3 (x,y) cordinates
#
# Arrows are polylines with 1 line segment and an arrow 
#  fig represents this by  a firstline record where
#    field 1 = 2 (always 2)
#    field 2 = 1 (polyline)
#    field 3 = 0 (style is a firm line)
#    field 14 = 1 for a forward arrow  
#    field 15 = 1 for a backward arrow 
#  an additional data files
#  a data field starting with a tab followed by 2(x,y) cordinates
#
#
# Strokes are polylines with 1 line segment and and no arrow 
#  fig represents this by  a firstline record where
#    field 1 = 2 (always 2)
#    field 2 = 1 (polyline)
#    field 3 = 0 (style is a firm line)
#    field 14 = 0 (no forward arrow)
#    field 15 = 0 (backward arrow) 
#  a data field starting with a tab followed by 2(x,y) cordinates
#
# Components appear in two files -- the fig file and the lbl file
# these two files are concatenated with the lbl file first
#  	The lbl file represents components by 3 fields
#		field  1 is the name
#		field  2 is the CR name
#		field  3 is the CR arguments
#  	The fig file represents components by 14 fields
#             field 1 = 4
#		fields 12 and 13 are the coordinates
#		field 14 is the type:name string terminated by \001
# To prevent text being confused with components, components consist
# of alphanumeric characters and : and _ only.
# The lbl file is used to sort the components.
##############################################################

function exact_match(name1, name2) {
  return ((match(name1,name2)>0)&&(length(name1)==length(name2)))
    }

function process_lbl() {
# This puts the components in the lable file at the top of the list
# and saves up the corresponding CR and arguments
# note that there may be more than one component per label
  if ((match($1,"%")==0)&&(NF>0))
    { 
      i_label++;
      name = $1;
      CR   = $2;
      args = $3;
      label[i_label,1] = name; 
      label[i_label,2] = CR;
      label[i_label,3] = args
	}
}

function fig_info() {
# Grabs the fig-file information for a component
  return(sprintf("%s %s %s %s %s %s %s %s %s %s %s ", \
		 $1, $2, $3, $4, $5, $6, $7, \
		 $8, $9, $10, $11))
	 }

function process_text() {
# The text string is field 14 onwards
  str = $14; 
  for (i=15; i<=NF; i++) {
    str = sprintf("%s %s", str, $i)
      }
# It is terminated by /001 - so delete this termination
  str = substr(str,1,length(str)-4);

# A component string contains only alphanumeric  _ and :
  isa_plain_component = match(str, component_regexp)==0;

# A port is an integer within [] and no alpha characters
  isa_port = (match(str, port_regexp)>0)&&(match(str, nonport_regexp)==0);

# A port component is SS followed by : followed by a port string
  isa_port_component = 0;
  if (match(str, delimiter)) {
    split(str,a,delimiter);
    isa_port_component = (exact_match(a[1], "SS"))&&
      (match(a[2], port_regexp)>0)
      }
    
# A component is a plain or a port component
  isa_component = isa_plain_component||isa_port_component;

# Coordinates in fields 12 & 13
  x_coord = $12;
  y_coord = $13;

# Do the ports
  if (isa_port) {
    i_port++;
    port_index = substr(str,2,length(str)-2);
    ports[i_port] = sprintf("%s %s %s", x_coord, y_coord, port_index);
  }

# Do the port components
  if (isa_port_component) {
    i_port_component++;
    # Port number is the bit between the []
    port_number  = substr(a[2],2,length(a[2])-2);
    x_port[port_number] = x_coord;
    y_port[port_number] = y_coord;
    info_port[port_number] = fig_info();
  }

# Do the plain components
  if (isa_plain_component) {
    i_text++;

# Get repetitions (if any)
    if (match(str,repetition_regexp) > 0) {
      split(str,b,repetition_delimiter);
      repetitions = b[2];
      str = b[1];
	}
    else {
      repetitions = "1";
    };

    named_component = (match(str,delimiter) > 0);
    if (named_component) {
      split(str,a,delimiter);
      type = a[1];
      name = a[2];

# Check  if name is in label file and if used already
      found = 0; name_used = 0;
      for (i=1; i<=i_label; i++) {
	lname = label[i,1];
	if ( exact_match(name,lname) ) {
	  found = 1;
	  if (name in used) {
	    name_used = 1;
	    CR = label[i,2];
	    args = label[i,3];
	  }
	  else {
	    used[name] = 1
	      }
	  break
	    }
      }

	if (!found) {
	  if (isa_plain_component) {
	    printf(warning_f, name)
	  }
	  i_label++;
	  CR = default_cr;
	  args = "";
	  label[i_label,1] = name; 
	  label[i_label,2] = CR;
	  label[i_label,3] = args
	    }

# Give it a new entry if already used
      if (name_used) {
	i_label++;
	i_name++;
	name = sprintf("%1.0f", i_name);
	label[i_label,1] = name; 
	label[i_label,2] = CR;
	label[i_label,3] = args
	  }
    }

# Unnamed component
    if (named_component==0) {
      i_name++;
      name = sprintf("%1.0f", i_name);
      type = str;
      i_label++;
      label[i_label,1] = name;
      label[i_label,2] = default_cr;
      label[i_label,3] = default_args
	}

# Save in associative arrays by name	
    comp_type[name] = type;
    x[name] = x_coord;
    y[name] = y_coord;
    info[name] =  fig_info();
    reps[name] = repetitions;
  }

}

function process_bond() {

  arg_count++;
  if ( (arg_count-arrow)==1 ) 
    {

#Save up bond coords
      if (NF == (2*bond_coords+1) ) {
	i_bond++;
	bonds[i_bond] = sprintf("%s %s %s %s %s %s", \
				$2, $3, $4, $5, $6, $7);
      }
      
#Save up arrow coords
      if ( (arrow)&&(NF==(2*arrow_coords+1)) ) {
	i_arrow++;
	arrows[i_arrow] = sprintf("%s %s %s %s",  $2, $3, $4, $5);
      }
      
#Save up stroke coords
      if ( (!arrow)&&(NF==(2*stroke_coords+1)) ) {
	i_stroke++;
	strokes[i_stroke] = sprintf("%s %s %s %s",  $2, $3, $4, $5);
      }
    }
}

function write_fig() {
# Create _fig.fig file from _abg file - not components
  if ( (isa_fig_file)&&((object!=text)||(isa_component==0))) {
    if (exact_match($1,data_symbol)) {
      field_1 = out_data_symbol
	}
    else {
      field_1 = $1
	}

    printf field_1   >> fig_file
      for (i=2; i<=NF; i++)
	printf(" %s", $i)  >> fig_file;
    printf("\n") >> fig_file
      }
}

function process_fig() {
# Test for the fig format first line and data line
  data_line = (match($1,data_symbol)>0);
  first_line = (data_line==0)&&(NF>min_line_length);

#Process firstline
  if (first_line) {
    object = $1;
    sub_type = $2;
    style = $3;
    f_arrow = ($14==1)&&(object=polyline);
    b_arrow = ($15==1)&&(object=polyline);
    arrow = f_arrow||b_arrow;
    arg_count = 0;
  }

#Process text
  if (object==text) {
    process_text()
      }

# Process bond
  if ( \
(data_line)&& \
       (object==polyline)&& \
       (sub_type==sub_polyline)&& \
       (style==firm_style) \
       ) {
    process_bond()
      }   

  write_fig()

    }

BEGIN {
  sys_name = ARGV[1];
  delete ARGV[1];
  b_file = sprintf("%s_rbg.m", sys_name);
  c_file = sprintf("%s_cmp.m", sys_name);
  fig_file = sprintf("%s_fig.fig", sys_name);
  warning_f = "WARNING %s \t in fig file but not lbl file  - using\n";
  warning_l = "WARNING %s \t in lbl file but not fig file  - ignoring\n";
  warning_p = "WARNING system ports are not consecutively numbered\n";

  data_symbol = "----";
  out_data_symbol = "\t";
  default_cr = "";
  default_args = "";
  delimiter = ":";
  repetition_delimiter = "*";
  repetition_regexp = "\*";
  q = "\047";
  terminator = "\\001";
  component_regexp = "[^0-9a-zA-Z_:\*-]";
  port_regexp = "\[[0-9]*\]";
  nonport_regexp = "[a-zA-Z]";
    
  isa_fig_file = 0;
  min_line_length = 10;
  object = 0;
  polyline = 2;
  sub_polyline=1; 
  firm_style = 0;
  text = 4;
  bond_coords = 3;
  stroke_coords = 2;
  arrow_coords = 2;

  i_bond = 0;
  i_port = 0;
  i_stroke = 0;
  i_arrow = 0;
  i_label = 0;
  i_text = 0;
  i_name = 0;
  i_port_component = 0;

}
{
# Start of .fig file?
  if ( (NF>0) && (match("#FIG", $1) > 0) ) {
    isa_fig_file=1;
  }

  if (isa_fig_file==0) {
    process_lbl()    
      }
  else {
    process_fig()
      }
}

END {
#Print out the matlab functions
  printf("function [rbonds, rstrokes,rcomponents,rports,n_ports] = %s_rbg\n", sys_name) > b_file;
  printf("%% [rbonds,rstrokes,rcomponents,rports,n_ports] = %s_rbg\n", sys_name) > b_file;
  printf("%% Generated by MTT\n\n") > b_file;

  printf("function [comp_type, name, cr, arg, repetitions] = %s_cmp(i)\n",\
		 sys_name) > c_file;
  printf("%% [comp_type, name, cr, arg, repetitions] = %s_cmp\n", sys_name) > c_file;
  printf("%% Generated by MTT\n\n") > c_file;

  printf("rbonds = [\n") >> b_file;
  for (i = 1; i <= i_bond; i++)
    print  bonds[i] >> b_file;
  for (i = 1; i <= i_arrow; i++)
    print  arrows[i], "-1 -1" >> b_file;
  printf("];\n") >> b_file;
  
  printf("rstrokes = [\n") >> b_file;
  for (i = 1; i <= i_stroke; i++)
    print  strokes[i] >> b_file;
  printf("];\n") >> b_file;

  printf("rcomponents = [") >> b_file;
  j = 0;


# Do the port components, in order, first
  for (i = 1; i <= i_port_component; i++) {
    port_type = "SS";
    name = sprintf("[%1.0f]", i);
    cr   = i;
    arg  = "";

    if (length(x_port[i])==0)
      printf(warning_p);
    else {
      j++;
      print x_port[i], y_port[i], info_port[i] >> b_file;
      printf("if i==%1.0f\n", j) 	>> c_file;
      printf("\tcomp_type = %s%s%s;\n", q, port_type, q) >> c_file;
      printf("\tname = %s%s%s;\n", q, name, q) >> c_file;
      printf("\tcr = %s%s%s;\n", q, cr, q) >> c_file;
      printf("\targ = %s%s%s;\n", q, arg, q) >> c_file;
      print "end" >> c_file
	}
    
  }

# Now do the ordinary components (in no particular order)
  for (i = 1; i <= i_label; i++) {
    name = label[i,1];
    cr   = label[i,2];
    arg  = label[i,3];
    
    if (length(x[name])==0)
      printf(warning_l, name);
    else {
      j++;
      print x[name], y[name], info[name] >> b_file;
      printf("if i==%1.0f\n", j) 	>> c_file;
      printf("\tcomp_type = %s%s%s;\n", q, comp_type[name], q) >> c_file;
      printf("\tname = %s%s%s;\n", q, name, q) >> c_file;
      printf("\tcr = %s%s%s;\n", q, cr, q) >> c_file;
      printf("\targ = %s%s%s;\n", q, arg, q) >> c_file;
      printf("\trepetitions = %s;\n",  reps[name]) >> c_file;
      print "end" >> c_file
	}
  }
  printf("];\n") >> b_file;

# Print the ports list
  printf("rports = [\n") >> b_file;
  for (i = 1; i <= i_port; i++)
    print  ports[i] >> b_file;
  printf("];\n\n") >> b_file;

# Print the number of ports
  printf("n_ports = %1.0f;\n", i_port_component) >> b_file;

  

}

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

# gawk script: rep_txt2make
# Converts the text file describing a report to makefile dependencies
# P.J.Gawthrop August 1996
# Copyright (c) P.J.Gawthrop, 1996.

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


BEGIN {
  split(ARGV[1],a,"_");
  system_name = a[1];
    }
{
  if (NF==2) {
    if( match("tex txt r m ps",$2)>0) {
      printf("%s_%s.%s ", system_name, $1, $2)
    } 
  }
}
END {

}

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

# gawk script: rep_txt2tex
# Converts the text file describing a report to a report.
# P.J.Gawthrop August 1996
# Copyright (c) P.J.Gawthrop, 1996.

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


BEGIN {
  split(ARGV[1],a,"_");
  system_name = a[1];
  args = ARGV[2];
}
{
  if (NF==2) {
# tex files
    if( match("tex",$2)>0) {
      print "\\section{Representation " $1 ", language " $2 "}";
      printf("  \\input{%s_%s.%s}\n", system_name, $1, $2);
    }
# text files
    if( match("txt r m",$2)>0) {
      print "\\section{Representation " $1 ", language " $2 "}";
      print "  \\begin{verbatim}";
      command = sprintf("cat %s_%s.%s", system_name, $1, $2);
      system(command);
      print "  \\end{verbatim}";  
    }
# ps files
    if( match("ps",$2)>0) {
      print "\\section{Representation " $1 ", language " $2 "}";
      printf("This representation is given as Figure \\ref{fig:%s}.\n", $1);
      print "  \\begin{figure}";
      printf("    \\epsfig{file=%s_%s.%s,width=\\linewidth}\n", \
	     system_name,  $1, $2);
      printf("    \\caption{System %s, representation %s}\n", system_name, $1);
      printf("    \\label{fig:%s}\n", system_name, $1);
      print "  \\end{figure}";
    }
  }
}
END {

}

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

# gawk script: tex2doc
# Encapsulates a tex file in a document.
# P.J.Gawthrop August 1996
# Copyright (c) P.J.Gawthrop, 1996.

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


BEGIN {
  split(ARGV[1],a,"_");
  system_name = a[1];
  split(a[2],b,".");
  representation = b[1];
  makefile = sprintf("%s_rep.make", system_name, representation);
  print "\\documentclass[12pt,a4paper]{article}";
  printf("\\title{System %s, representation %s}\n", \
	 system_name, representation) 
  print "\\author{Generated by MTT}";
  print "  \\usepackage{epsfig}";
  print ;
  print "\\begin{document}";
  print "  \\maketitle";
  printf("  \\input{%s_%s}\n", system_name, representation);
  print "\\end{document}"
}

#! /bin/sh

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

# Bourne shell script: cbg2ese_m
#
#
# P.J.Gawthrop June 1996
# Copyright (c) P.J.Gawthrop, 1996.

###############################################################
## Version control history
###############################################################
## $Id$
##
## $Log$
# Revision 1.5  1996/08/24  14:38:11  peter
# Error handling included.
#
## Revision 1.4  1996/08/19 09:02:12  peter
## Errors go to user - not log.
##
## Revision 1.3  1996/08/18 20:05:55  peter
## Included zero outputs.
##
## Revision 1.2  1996/08/08 18:08:37  peter
## Sorted out file naming scheme.
##
## Revision 1.1  1996/08/08 15:57:15  peter
## Initial revision
##
###############################################################

infofile='mtt_info.txt';
eqnfile="$1_ese.r";
deffile="$1_def.r";

# Remove the old log file
rm -f cbg2ese_m.log
rm -f $1_ese.r
rm -f $1_def.r
rm -f $infofile

touch $infofile

#Inform user
echo Creating $eqnfile
echo Creating $1_def.r

# Use matrix manipulation to accomplish the transformation
$MATRIX > cbg2ese_m2r.log 2>mtt_error << EOF

infofile = '$infofile'
system_name = '$1'
deffile = '$1_def.r';

structure = cbg2ese(system_name,'','',infofile)
makedef(structure,deffile);
  
EOF

cat $infofile

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

#! /bin/sh

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

# Bourne shell script: cbg_m2fig
#
# Causal bond graph:  mfile format to fig file format
# The resultant fig file is the original _abg.fig with
# additional causal strokes superimposed.
#
# P.J.Gawthrop May 1996
# Copyright (c) P.J.Gawthrop, 1996.

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.3  1996/08/24 14:40:31  peter
## Error handling included.
##
## Revision 1.2  1996/08/08 18:07:31  peter
## Sorted out naming scheme for files
##
## Revision 1.1  1996/08/05 20:13:56  peter
## Initial revision
##
###############################################################


# Remove the old log file
rm -f cbg_m2fig.log
rm -f mtt_info.txt

# The following sets up the fig header file for the system
# and each sub system.
(ls $1_*_cbg.m \
| sed 's/\(.*_\)\(.*_\)cbg.m/cp \2fig.fig \1\2cbg.fig/' \
| sh ) 2>mtt_error.txt


#Inform user
echo Creating $1_cbg.fig

# Use  matrix manipulation to accomplish the transformation
$MATRIX > cbg_m2fig.log  2>>mtt_error.txt << EOF

  system_name = '$1';
  thick = 3; length = 250; red = 4; blue = 1; green = 12; font = 18;
  cbg2fig(system_name, ...
          '', '', ...
          length, thick, blue, ...
          font, green, red);
EOF

cp $1_$1_cbg.fig $1_cbg.fig 2>> mtt_error.txt

# Print errors, if any.
mtt_error mtt_error.txt

#! /bin/sh

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

# Bourne shell script: cse2csm_r
# Constrained-state equation to linear constrained-state matrices conversion
# P.J.Gawthrop  6th September 1991, May 1994
# Copyright (c) P.J.Gawthrop, 1991, 1994.

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

# Inform user
echo Creating $1_csm.r

# Remove the old log file
rm -f cse2csm_r.log

# Use reduce to accomplish the transformation
reduce >cse2csm_r.log << EOF

in "$1_def.r";
in "$1_cse.r";
in "$1_cr.r";
in "$1_sympar.r";

OFF Echo;
OFF Nat;

% Find MTTA : the A matrix
matrix MTTA(MTTNx,MTTNx);

  FOR j := 1:MTTNx DO
    BEGIN
    xj := MTTX(j,1);
    FOR i := 1:MTTNx DO
      MTTA(i,j) := df(MTTEdx(i,1), xj, 1);
    END;

% Find MTTB : the B matrix
matrix MTTB(MTTNx,MTTNu);

  FOR j := 1:MTTNu DO
    BEGIN
    uj := MTTU(j,1);
    FOR i := 1:MTTNx DO
      MTTB(i,j) := df(MTTEdx(i,1), uj, 1);
    END;

% Find MTTC : the C matrix
matrix MTTC(MTTNy,MTTNx);

  FOR i := 1:MTTNy DO
    FOR j := 1:MTTNx DO
    BEGIN
      xj := MTTX(j,1);
      MTTC(i,j) := df(MTTY(i,1), xj, 1);
    END;

% Find MTTD : the D matrix
matrix MTTD(MTTNy,MTTNu);

  FOR i := 1:MTTNy DO
    FOR j := 1:MTTNu DO
    BEGIN
      xj := MTTU(j,1);
      MTTD(i,j) := df(MTTY(i,1), xj, 1);
    END;

%Substitute the ss values
in "$1_ss.r";

%Create the output file
OUT "$1_csm.r";

%Write out the matrices.
IF MTTNx>0 THEN
BEGIN
  write "matrix MTTE(", MTTNx, ",", MTTNx, ");";
  FOR i := 1:MTTNx DO
   FOR j := 1:MTTNx DO IF MTTE(i,j) NEQ 0 THEN
      write "MTTE(", i, ",", j, ") := ", MTTE(i,j);

  write "matrix MTTA(", MTTNx, ",", MTTNx, ");";
  FOR i := 1:MTTNx DO
   FOR j := 1:MTTNx DO IF MTTA(i,j) NEQ 0 THEN
      write "MTTA(", i, ",", j, ") := ", MTTA(i,j);
END;

IF MTTNx>0 THEN
IF MTTNu>0 THEN
BEGIN
  write "matrix MTTB(", MTTNx, ",", MTTNu, ");";
  FOR i := 1:MTTNx DO
   FOR j := 1:MTTNu DO IF MTTB(i,j) NEQ 0 THEN
      write "MTTB(", i, ",", j, ") := ", MTTB(i,j);
END;

%Write it out
IF MTTNy>0 THEN
IF MTTNx>0 THEN
BEGIN
  write "matrix MTTC(", MTTNy, ",", MTTNx, ");";
  FOR i := 1:MTTNy DO
    FOR j := 1:MTTNx DO IF MTTC(i,j) NEQ 0 THEN
      write "MTTC(", i, ",", j, ") := ", MTTC(i,j);
END;

IF MTTNy>0 THEN IF MTTNu>0 THEN 
BEGIN
  write "matrix MTTD(", MTTNy, ",", MTTNu, ");";
  FOR i := 1:MTTNy DO
    FOR j := 1:MTTNu DO IF MTTD(i,j) NEQ 0 THEN
      write "MTTD(", i, ",", j, ") := ", MTTD(i,j);
END;

write "END;";

SHUT "$1_csm.r";

EOF

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

#! /bin/sh

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

# Bourne shell script: dae2daesol_m

# Transforms descriptor matrix rep to step response

# Copyright (c) P.J.Gawthrop, 1996.

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.3  1996/08/18 12:00:19  peter
## Unified format of responses.
##
## Revision 1.2  1996/08/16 13:17:57  peter
## Changed default args to include nx+nz
## Fixed bug with vector outputs.
##
## Revision 1.1  1996/08/15 16:24:55  peter
## Initial revision
##
###############################################################


echo Creating $1_daes.m
echo Creating $1_daeso.m
rm -f dae2daes_m.log
rm -f $1_daes.m
rm -f $1_daeso.m

if [ "$2" = "" ]; 
then
  PARAMS='T=[0:0.1:10];x0=zeros(nx+nz,1);dx0=zeros(nx+nz,1);'
  echo Using default parameter $PARAMS
else
  PARAMS=$2;
fi

PARAMS="$PARAMS ;"


$MATRIX << EOF > dae2daes_m.log 2>mtt_error

  %Read the parameters
  $1_numpar;

  [nx,ny,nu,nz,nyz] = $1_def;
  t=0;  %Just in case its in the parameter list
  $PARAMS

  %Defaults
  if exist('T')==0
    T=[0:0.1:10]
  end;

  if exist('x0')==0
    x0 = zeros(nx+nz,1);
  end;

  if exist('dx0')==0
   dx0 = zeros(nx+nz,1);
  end;

  [n,m]=size(T);
  if m>n
    T=T';
  end;


if nx>0
  x = dassl('$1_dae', x0, dx0,  T);
  write_matrix([T,x], '$1_daes');
else
  x = zeros(size(T));
end;

if ny>0
  i=0;
  for tt=T'
    i=i+1;
    y(i,:) = $1_daeo(x(i,:),tt)';
  end;
  write_matrix([T,y], '$1_daeso');
end;

EOF

# Test for errors and print if any
err_length=$(wc -c <mtt_error)
if [ $err_length != "0" ]
then
  echo MTT has failed with the following errors '...'
  cat mtt_error
  exit 1
else
  exit 0
fi

#! /bin/sh

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

# Bourne shell script: dae2dm_r
# Differential-algebraic equation to linear descriptor matrices conversion
# P.J.Gawthrop  8th May 1991, 7 June 1991, April 1994.
# Copyright (c) P.J.Gawthrop, 1991, 1994.

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.1  1996/08/18 20:03:51  peter
## Initial revision
##
###############################################################

#Inform user
echo Creating $1_dm.r

# Remove the old log file
rm -f dae2dm_r.log

# Use reduce to accomplish the transformation
reduce >dae2dm_r.log << EOF

in "$1_def.r";
in "$1_dae.r";
in "$1_cr.r";
in "$1_sympar.r";

OFF Echo;
OFF Nat;


%Total number of descriptor variables
MTTN := MTTNx + 2*MTTNz + MTTNyz;

IF MTTN>0 THEN
BEGIN
% Find MTTE : the E matrix
matrix MTTE(MTTN,MTTN);

% First row partition
IF MTTNx>0 THEN 
    FOR j := 1:MTTNx DO
      MTTE(j,j) := 1;

% Second row partition
  IF MTTNz>0 THEN 
     FOR j := 1:MTTNz DO
      MTTE(j+MTTNx, j+MTTNx) := 1;

% Find MTTA : the A matrix
matrix MTTA(MTTN,MTTN);

% First row partition: x
IF MTTNx>0 THEN 
  BEGIN
  FOR j := 1:MTTNx DO
    BEGIN
    xj := MTTX(j,1);
    FOR i := 1:MTTNx DO
      MTTA(i,j) := df(MTTdX(i,1), xj, 1);
    END;

  IF MTTNz>0 THEN 
     FOR j := 1:MTTNz DO
       BEGIN
       dzj := MTTdz(j,1);
       FOR i := 1:MTTNx DO
         MTTA(i, j + MTTNx + MTTNz) := df(MTTdX(i,1), dzj, 1);
       END;

  IF MTTNyz>0 THEN 
     FOR j := 1:MTTNyz DO
       BEGIN
       uij := MTTUi(j,1);
       FOR i := 1:MTTNx DO
         MTTA(i, j + MTTNx + 2*MTTNz) := df(MTTdX(i,1), uij, 1);
       END;
  END;

% Second row partition: z
IF MTTNz>0 THEN 
  FOR j := 1:MTTNz DO
    MTTA(j + MTTNx, j + MTTNx + MTTNz) := 1;

% Third row partition
IF MTTNz>0 THEN 
  BEGIN
  IF MTTNx>0 THEN 
  FOR j := 1:MTTNx DO
    BEGIN
    xj := MTTX(j,1);
    FOR i := 1:MTTNz DO
      MTTA(i + MTTNx + MTTNz,j) := df(MTTZ(i,1), xj, 1);
    END;

   FOR j := 1:MTTNz DO
	MTTA(j + MTTNx + MTTNz, j + MTTNx) := -1;

   FOR j := 1:MTTNz DO
     BEGIN
     dzj := MTTdz(j,1);
     FOR i := 1:MTTNz DO
       MTTA(i + MTTNx + MTTNz, j + MTTNx + MTTNz) := df(MTTZ(i,1), dzj, 1);
     END;

  IF MTTNyz>0 THEN 
     FOR j := 1:MTTNyz DO
       BEGIN
       uij := MTTUi(j,1);
       FOR i := 1:MTTNz DO
         MTTA(i+MTTNx+MTTNz, j+MTTNx+2*MTTNz) := df(MTTZ(i,1), uij, 1);
       END;
  END;


% Fourth row partition
IF MTTNyz>0 THEN 
  BEGIN
  IF MTTNx>0 THEN 
  FOR j := 1:MTTNx DO
    BEGIN
    xj := MTTX(j,1);
    FOR i := 1:MTTNyz DO
      MTTA(i + MTTNx + 2*MTTNz,j) := df(MTTYz(i,1), xj, 1);
    END;

  IF MTTNz>0 THEN 
  FOR j := 1:MTTNz DO
    BEGIN
     dzj := MTTdz(j,1);
    FOR i := 1:MTTNyz DO
      MTTA(i + MTTNx + 2*MTTNz, j + MTTNx + MTTNz) := df(MTTYz(i,1), dzj, 1);
    END;

  IF MTTNyz>0 THEN 
     FOR j := 1:MTTNyz DO
       BEGIN
       uij := MTTUi(j,1);
       FOR i := 1:MTTNyz DO
         MTTA(i+MTTNx+2*MTTNz, j+MTTNx+2*MTTNz) := df(MTTYz(i,1), uij, 1);
       END;
  END;

% Find MTTB : the B matrix
matrix MTTB(MTTN,MTTNu);

IF MTTNu>0 THEN
BEGIN

% First row partition
IF MTTNx>0 THEN 
  BEGIN
  FOR j := 1:MTTNu DO
    BEGIN
    uj := MTTU(j,1);
    FOR i := 1:MTTNx DO
      MTTB(i,j) := df(MTTdX(i,1), uj, 1);
    END;
  END;

% Third row partition
IF MTTNz>0 THEN 
  BEGIN
  FOR j := 1:MTTNu DO
    BEGIN
    uj := MTTU(j,1);
    FOR i := 1:MTTNz DO
      MTTB(i + MTTNx + MTTNz ,j) := df(MTTZ(i,1), uj, 1);
    END;
  END;

% Fourth row partition
IF MTTNyz>0 THEN 
  BEGIN
  FOR j := 1:MTTNu DO
    BEGIN
    uj := MTTU(j,1);
    FOR i := 1:MTTNyz DO
      MTTB(i + MTTNx + 2*MTTNz,j) := df(MTTYz(i,1), uj, 1);
    END;
  END;

END;

% Find MTTC : the C matrix
matrix MTTC(MTTNy,MTTN);
IF MTTNy>0 THEN
BEGIN

% First column partition
IF MTTNx>0 THEN 
  BEGIN
  FOR i := 1:MTTNy DO
    FOR j := 1:MTTNx DO
    BEGIN
      xj := MTTX(j,1);
      MTTC(i,j) := df(MTTY(i,1), xj, 1);
    END;
  END;

% Third column partition
IF MTTNz>0 THEN 
  BEGIN
  FOR i := 1:MTTNy DO
    FOR j := 1:MTTNz DO
    BEGIN
      dzj := MTTdZ(j,1);
      MTTC(i, j + MTTNx + MTTNz) := df(MTTY(i,1), dzj, 1);
    END;
  END;

% Fourth column partition
IF MTTNyz>0 THEN 
  BEGIN
  FOR i := 1:MTTNy DO
    FOR j := 1:MTTNyz DO
    BEGIN
      uij := MTTUi(j,1);
      MTTC(i, j + MTTNx + 2*MTTNz) := df(MTTY(i,1), uij, 1);
    END;
  END;

END;
END; %of MTTN>0

% Find MTTD : the D matrix
matrix MTTD(MTTNy,MTTNu);
IF MTTNy>0 THEN
IF MTTNu>0 THEN
BEGIN
  FOR i := 1:MTTNy DO
    FOR j := 1:MTTNu DO
    BEGIN
      uj := MTTU(j,1);
      MTTD(i,j) := df(MTTY(i,1), uj, 1);
    END;
END;

%Substitute the ss values
in "$1_ss.r";

%Create the output file
OUT "$1_dm.r";

%Write out the matrices.
IF MTTN>0 THEN
BEGIN
  write "matrix MTTE(", MTTN, ",", MTTN, ");";
  FOR i := 1:MTTN DO
   FOR j := 1:MTTN DO IF MTTE(i,j) NEQ 0 THEN
      write "MTTE(", i, ",", j, ") := ", MTTE(i,j);

  write "matrix MTTA(", MTTN, ",", MTTN, ");";
  FOR i := 1:MTTN DO
   FOR j := 1:MTTN DO IF MTTA(i,j) NEQ 0 THEN
      write "MTTA(", i, ",", j, ") := ", MTTA(i,j);
END;

IF MTTN>0 THEN
IF MTTNu>0 THEN
BEGIN
  write "matrix MTTB(", MTTN, ",", MTTNu, ");";
  FOR i := 1:MTTN DO
   FOR j := 1:MTTNu DO IF MTTB(i,j) NEQ 0 THEN
      write "MTTB(", i, ",", j, ") := ", MTTB(i,j);
END;

%Write it out
IF MTTNy>0 THEN
IF MTTN>0 THEN
BEGIN
  write "matrix MTTC(", MTTNy, ",", MTTN, ");";
  FOR i := 1:MTTNy DO
    FOR j := 1:MTTN DO IF MTTC(i,j) NEQ 0 THEN
      write "MTTC(", i, ",", j, ") := ", MTTC(i,j);
END;

IF MTTNy>0 THEN IF MTTNu>0 THEN 
BEGIN
  write "matrix MTTD(", MTTNy, ",", MTTNu, ");";
  FOR i := 1:MTTNy DO
    FOR j := 1:MTTNu DO IF MTTD(i,j) NEQ 0 THEN
      write "MTTD(", i, ",", j, ") := ", MTTD(i,j);
END;

write "END;";

SHUT "$1_dm.r";
quit;

EOF

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

#! /bin/sh

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

# Bourne shell script: dae_r2m
# Reduce DAE to simulab  DAE  
# P.J.Gawthrop 14 June 1991, 12 Jan 1994, April 1994, Jan 95.
# Copyright (c) P.J.Gawthrop 1991, 1994, 1995, 1996

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.4  1996/08/18 20:02:21  peter
## Include zero outputs.
##
## Revision 1.3  1996/08/16 13:10:53  peter
## Fixed bug Ny changed to Nx in redidual loop.
##
## Revision 1.2  1996/08/16 08:57:34  peter
## Reads $1_input.
##
## Revision 1.1  1996/08/15 16:46:06  peter
## Initial revision
##
###############################################################


# WARNING this is experimental!!
# It asumes that the output y is NOT dependent on derivatives of
#  z or u.


#Inform user
echo Creating $1_dae.m
echo Creating $1_daeo.m


# Remove the old log file
rm -f dae_r2m.log

# Use reduce to accomplish the transformation
reduce >dae_r2m.log << EOF

%Read the reduce definitions file
in "$1_def.r";

%Read the reduce  DAE   file
in "$1_dae.r";

%Set up the number of argument variables to zero in case the user has forgotten
MTTNVar := 0;

%Read the parameter file
in "$1_sympar.r";


ON NERO;        % Suppress zero elements

%Define the common part of the functions.

PROCEDURE common;
BEGIN
  IF MTTNvar>0 THEN
  BEGIN
    write "% Set the parameters";
    write "global ...;;";
    FOR i := 1:MTTNvar DO
    BEGIN
       IF i<MTTNvar THEN write MTTVar(i,1), " ..."
                    ELSE write MTTVar(i,1), ";"
    END;
  END;
  
  write "% Read in the input";
  write "u = $1_input(t)";

  write "% Read in the arguments";
  write "$1_args";
  
  write "% Set up the State variables";
  FOR i := 1:MTTNx DO
  BEGIN
    write "MTTx", i, " = MTTx(", i, ");";
  END;

  write "% Set up the non-state variables";
  FOR i := 1:MTTNz DO
  BEGIN
    write "MTTz", i, " = MTTx(", MTTNx+ i, ");";
  END;

   write "% Set up the internal inputs";
  FOR i := 1:MTTNz DO
  BEGIN
    write "MTTui", i, " = MTTx(", MTTNx+MTTz+i, ");";
  END;
  
  write "% Set up the Input variables";
  IF MTTNu>0 THEN
  FOR i := 1:MTTNu DO
  BEGIN
    write "MTTu", i, " = u(", i, ");";
  END;
END;
  
% Firstly do the residual = f(dx,x,t) function.

%MTTGx and u are the derivatives of z wrt x and u respectively

% Find MTTGx;
matrix MTTGx(MTTNz,MTTNx);
FOR j := 1:MTTNx DO
  BEGIN
  xj := MTTX(j,1);
  FOR i := 1:MTTNz DO
    MTTGx(i,j) := df(MTTZ(i,1), xj, 1);
  END;

% Find MTTGu;
matrix MTTGu(MTTNz,MTTNu);
FOR j := 1:MTTNu DO
  BEGIN
  uj := MTTu(j,1);
  FOR i := 1:MTTNz DO
    MTTGu(i,j) := df(MTTZ(i,1), uj, 1);
  END;


OUT "$1_dae.m";
write "function residual = $1_dae(MTTdX0,MTTx,t);";
write "% residuals = $1_dae(dx,x,t);";
write "%DAE in DASSL $1;;";
write "%File $1_dae.m;;";
write "%Generated by MTT;;";


common();
  write "% Set up the non-state derivatives";
  FOR i := 1:MTTNz DO
  BEGIN
    write "MTTdz", i, " = MTTdx0(", MTTNx+ i, ");";
  END;
  

%Fortran switches - one line expressions
OFF echo;
ON fort$
cardno!* := 1$
fortwidth!* := 100$
OFF period$


MTTdx := MTTdx;
MTTGx := MTTGx;

OFF fort;

write "for i=1:", MTTNx, ";;";
write "  residual(i) = MTTdX(i)-MTTdX0(i);";
write "end;";

IF MTTNz>0 THEN
BEGIN
  write "dxz = MTTGx*MTTdx";
  write "for i=1:", MTTNz, ";;";
  write "  residual(i+", MTTNx, ") = dxz(i)-MTTdX0(i+", MTTNx,");";
  write "end;";
END;



SHUT "$1_dae.m";

OFF fort;

% Now do the y = g(x,t) function.
OUT "$1_daeo.m";
  
write "function MTTy = $1_daeo(MTTx,t);";
write "% dX = $1_daeo(MTTx,t);";
write "%DAE in Simulab form for system $1;;";
write "%File $1_daeo.m;;";
write "%Generated by MTT;;";

common();

%Fortran switches - one line expressions
OFF echo;
ON fort$
cardno!* := 1$
fortwidth!* := 100$
OFF period$

MTTy := MTTy;

SHUT "$1_daeo.m";

#! /bin/sh

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

# Bourne shell script: dm2fr_m

# Transformation: descriptor matrix to frequency response

# Copyright (c) P.J.Gawthrop, 1996.

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.6  1996/08/16 14:26:07  peter
## Fixed multi-output bug.
## Made sure that u0 is ok.
##
## Revision 1.5  1996/08/15 16:23:02  peter
## Now uses W in place of w for consistancy with time responses.
##
## Revision 1.4  1996/08/15 11:52:42  peter
## Now creats a number of versions:
## fr 	complex frequency response
## lmfr 	log  magnitude of fr
## lpfr  	phase of fr
## nyfr  	real and imag parts - Nyquist style
## nifr 	Nichols style.
##
## Revision 1.3  1996/08/11 19:08:08  peter
## Parameter passing now enabled.
##
## Revision 1.2  1996/08/11 09:32:12  peter
## Now takes the numpar parameters correctly
##
## Revision 1.1  1996/08/10 14:11:11  peter
## Initial revision
##
###############################################################

echo Creating $1_fr.m
echo Creating $1_lmfr.m
echo Creating $1_lpfr.m
echo Creating $1_nyfr.m
echo Creating $1_nifr.m

if [ "$2" = "" ]; 
then
  PARAMS='W=logspace(-1,2,100); u0=zeros(nu,1); u0(1)=1;'
  echo Using default parameter $PARAMS
else
  PARAMS=$2;
fi

$MATRIX << EOF > dm2fr_m.log 2>mtt_error

  $1_numpar;

  [nx,ny,nu,nz,nyz] = $1_def;
  $PARAMS

  %Defaults
  if exist('W')==0
    W = logspace(-1,2,100)';
  end;

  if exist('u0')==0
    u0 = zeros(nu,1);
    u0(1) = 1;
  end;

  [n,m]=size(W);
  if m>n
    W=W';
  end;

  [A,B,C,D,E] = $1_dm;
  fr = dm2fr(A,B,C,D,E,W,u0);
  lw = log10(W);
  lmfr = log10(abs(fr));
  pfr = angle(fr)*180/pi;

% Complex frequency response
  write_matrix([W fr], '$1_fr');

% Log magnitude v log frequency
  write_matrix([lw lmfr], '$1_lmfr');

% Angle v log frequency
  write_matrix([lw pfr], '$1_lpfr');

% Nyquist style
  re = real(fr);
  im = imag(fr);
  nyq = [re(:,1) im(:,1)];
  for i = 2:ny
    nyq = [nyq re(:,i) im(:,i)]
  end;
  write_matrix(nyq, '$1_nyfr');

% Nichols style
  re = lmfr;
  im = pfr;
  nic = [im(:,1) re(:,1)];
  for i = 2:ny
    nic = [nic im(:,i) re(:,i)]
  end;
  write_matrix(nic, '$1_nifr');

EOF

# Test for errors and print if any
err_length=$(wc -c <mtt_error)
if [ $err_length != "0" ]
then
  echo MTT has failed with the following errors '...'
  cat mtt_error
  exit 1
else
  exit 0
fi

#! /bin/sh

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

# Bourne shell script: dat2ps

# Tramsforms descriptor matrix rep to frequency response dat.


# Copyright (c) P.J.Gawthrop, 1996.

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.2  1996/08/11 10:40:00  peter
## Now handles e and NaN numbers.
##
## Revision 1.1  1996/08/10 14:15:03  peter
## Initial revision
##
###############################################################

echo Creating $1_ir.dat

if [ "$2" = "" ]; 
then
  PARAMS='t=[0:0.1:10]'
  echo Using default parameter $PARAMS
else
  PARAMS=$2;
fi


$MATRIX << EOF | awk '$1 !~ /[A-MO-Zb-cf-z]/ && NF>0 {print} '  >$1_ir.dat

  $PARAMS;
  [n,m]=size(t);
  if m>n
    t=t';
  end;

  [A,B,C,D,E] = $1_dm($1_numpar);
  y = dm2ir(A,B,C,D,E,t);

  [Ny,Nu]=size(D);
  NN = Ny*Nu;

  xy = [t y(:,1)];
  for i=2:NN
    xy = [xy t y(:,i)];
  end;
xy

EOF

#! /bin/sh

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

# Bourne shell script: dm2tf_r
# Reduce  descriptor matrices to transfer function
# P.J.Gawthrop 8th May 1991, Dec 1993, April 1994.
# Copyright (c) P.J.Gawthrop, 1991, 1993, 1994.

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

#Inform user
echo Creating $1_tf.r

# Remove the old log file
rm -f dm2tf_r.log

# Use reduce to accomplish the transformation
reduce >dm2tf_r.log << EOF

%ON FLOAT;

IN "$1_def.r";
IN "$1_dm.r";

OFF Echo;
OFF Nat;

%create sE-A
%MTT_SEA := s*MTTE-MTTA;

%Find the denominator of the TF - det(sE-A);
%comden := det(MTT_SEA);

%Find the Adjoint transpose.
%matrix AdjT(MTTNx,MTTNx);
%FOR i := 1:MTTNx DO
%  BEGIN
%  FOR j := 1:MTTNx DO
%    AdjT(i,j) := cofactor(MTT_SEA,i,j);
%  END;

%Adj := TP(AdjT);

%Find the numerator matrix
%Num := MTTC*Adj*MTTB + MTTD*comden;


%Create the transfer function matrix
MTTTF := MTTD;
IF MTTNy>0 THEN MTTTF := MTTTF + (MTTC * ((s*MTTE-MTTA)^-1) * MTTB);

%MTTTF := Num/comden;


OUT "$1_tf.r";

%Declare the transfer function matrix
write "matrix MTTTF(", MTTNy, ",", MTTNu, ")$"$

%And write it.
%MTTTF := MTTTF;
FOR i := 1:MTTNy DO
  BEGIN
  FOR j := 1:MTTNu DO
    IF MTTTF(i,j) NEQ 0 THEN write "MTTTF(", i, ",", j, ") := ", MTTTF(i,j)$
  END;

write ";END;"$
SHUT "$1_tf.r";

EOF

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

#! /bin/sh

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

# Bourne shell script: dm_r2m
# Reduce descriptor state-space A,B,C,D and E matrices to Matlab
# P.J.Gawthrop 29th July 1990, 6th Feb 1991, 28 May 1991, Dec 1993, Jan 94
# Copyright (c) P.J.Gawthrop, 1989, 1990, 1991, 1993,1994.

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.3  1996/08/18 12:02:46  peter
## Fixed bug when some matrices are empty.
##
## Revision 1.2  1996/08/12 20:22:20  peter
## Paramaters passed via an internal call to _numpar.
## _args script can overide these parameters - sneaky.
##
## Revision 1.1  1996/08/12 19:09:13  peter
## Initial revision
##
###############################################################

#Inform user
echo Creating $1_dm.m

# Remove the old log file
rm -f dm_r2m.log

# Use reduce to accomplish the transformation
reduce >dm_r2m.log << EOF

ON BigFloat, NumVal;
PRECISION 16; %Compatible with Matlab

%Read in the definitions file
in "$1_def.r";

%Set up the number of argument variables to zero in case the user has forgotten
MTTNVar := 0;

%Read in the parameter file
in "$1_sympar.r";

%Read the reduce state-space A,B,C and D matrices file
in "$1_dm.r";


OUT "$1_dm.m";

%Headings - Matlab style
%(Note. The ;; are deleted by for2mat)

write "function [MTTA,MTTB,MTTC,MTTD,MTTE] = $1_dm;;";
write "% [MTTA,MTTB,MTTC,MTTD,MTTE] = $1_dm;;";
write "% Linearised descriptor matrices for system $1";
write "% File $1_dm.m";
write "% Generated by MTT";

IF MTTNvar>0 THEN
BEGIN
  write "% Parameters";
  write "global ...;;";
  FOR i := 1:MTTNvar DO
  BEGIN
     IF i<MTTNvar THEN write MTTVar(i,1), " ..."
                  ELSE write MTTVar(i,1), ";"
  END;
END;


%Fortran switches - one line expressions
OFF echo;
ON fort$
cardno!* := 1$
fortwidth!* := 100$
OFF period$

MTTN := MTTNx + 2*MTTNz + MTTNyz;
write "MTTA = zeros(", MTTN, ",", MTTN, ");";
write "MTTB = zeros(", MTTN, ",", MTTNu, ");";
write "MTTC = zeros(", MTTNy, ",", MTTN, ");";
write "MTTD = zeros(", MTTNy, ",", MTTNu, ");";
write "MTTE = zeros(", MTTN, ",", MTTN, ");";


ON NERO;	% Suppress zero elements.
MTTA := MTTA;
MTTB := MTTB;
MTTC := MTTC;
MTTD := MTTD;
MTTE := MTTE;

SHUT "$1_dm.m";

#! /bin/sh

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

# Bourne shell script: gdat2ps

# Converts a data file in gplot format to a ps file

# Copyright (c) P.J.Gawthrop, 1996.

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
# Revision 1.6  1996/08/16  14:51:21  peter
# Put in some cosmetic postscript options.
#
## Revision 1.5  1996/08/15 09:37:58  peter
## Parameters in title.
##
## Revision 1.4  1996/08/14 09:21:17  peter
## Gnu plot conversion now done externally in dat2gdat
##
## Revision 1.3  1996/08/10 14:12:48  peter
## Revised to new dat format: [x y x y x y] in place of [x y y y].
##
## Revision 1.2  1996/08/10 12:54:56  peter
## Reorganised data file for gnuplot input.
##
## Revision 1.1  1996/08/10 09:52:39  peter
## Initial revision
##
###############################################################

echo Creating $1.ps

gnuplot << EOF
  set terminal postscript eps color "Times-Roman" 24
  set size 1.5,1.5
  set output '$1.ps'
  set grid
  set title "`date` $2 "
  plot '$1.gdat' with lines
  exit
EOF

#! /bin/sh

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

# Bourne shell script: gdat2view

# Converts a data file in gplot format to  view 

# Copyright (c) P.J.Gawthrop, 1996.

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
# Revision 1.2  1996/08/15  09:36:00  peter
# Parameters in title.
#
## Revision 1.1  1996/08/15 08:50:20  peter
## Initial revision
##
###############################################################


# echo Creating view of $1

gnuplot << EOF
  set grid
#  set title "$1. Generated by MTT on $(date)"
  set title "`date` $2 "
  plot '$1.gdat' with lines
  pause 100000
EOF

#! /bin/sh

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

# Bourne shell script: latex_tidy
# Tidies up LaTex format files
# P.J.Gawthrop 12 June 1990, Dec 1994
# Copyright (c) P.J.Gawthrop, 1990, 1994

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.1  1996/08/18 20:03:12  peter
## Initial revision
##
###############################################################

cat $1 | \
tr "[A-Z]" "[a-z]" | \
sed -e "s/\// \\\over /g" \
    -e "s/mttex/\\\chi/g" \
    -e "s/mttdx/\\\dot x/g" \
    -e "s/mttdz/\\\dot z_/g" \
    -e "s/mttdu/\\\dot u_/g" \
    -e "s/dot_/ \\\dot /g" \
    -e "s/mtta/A/g" \
    -e "s/mttb/B/g" \
    -e "s/mttc/C/g" \
    -e "s/mttd/D/g" \
    -e "s/mtte/E/g" \
    -e "s/mttg/G/g" \
    -e "s/mtti/I/g" \
    -e "s/mttm/M/g" \
    -e "s/mtto/O/g" \
    -e "s/mttv/V/g" \
    -e "s/mttx0/X_0/g" \
    -e "s/mttx1/x_1/g" \
    -e "s/mttx2/x_2/g" \
    -e "s/mttx3/x_3/g" \
    -e "s/mttx4/x_4/g" \
    -e "s/mttx5/x_5/g" \
    -e "s/mttx6/x_6/g" \
    -e "s/mttui0/v_0/g" \
    -e "s/mttui1/v_1/g" \
    -e "s/mttui2/v_2/g" \
    -e "s/mttui3/v_3/g" \
    -e "s/mttui4/v_4/g" \
    -e "s/mttui5/v_5/g" \
    -e "s/mttui6/v_6/g" \
    -e "s/mttui7/v_7/g" \
    -e "s/mttui8/v_8/g" \
    -e "s/mttui9/v_9/g" \
    -e "s/mttu/u_/g" \
    -e "s/mtty0/y_0/g" \
    -e "s/mttu0/u_0/g" \
    -e "s/mttx/x/g" \
    -e "s/mttz/z/g" \
    -e "s/mttu/u/g" \
    -e "s/mttyy/Y/g" \
    -e "s/mtty/y/g" \
    -e "s/mttw/w/g" \
    -e "s/mtttf/G/g" \
    -e "s/mttpar//g" \
    -e "s/\*\*/\^/g" \
    -e "s/\*/ /g" \
    -e "s/!/\\\/g" \
    -e "s/\\$//g" \
    -e "s/(/{(/g" \
    -e "s/)/)}/g" \
    -e "s/\left/\left (/g" \
    -e "s/\right/\right )/g" \
    -e "s/sqrt/\\\sqrt/g" \
    -e "s/theta/\\\theta/g" \
    -e "s/tau/\\\tau/g" \
    -e "s/alpha/\\\alpha/g" \
    -e "s/beta/\\\beta/g" \
    -e "s/epsilon/\\\epsilon/g" \
    -e "s/lambda/\\\lambda/g" \
    -e "s/mu/\\\mu/g" \
    -e "s/phi/\\\phi/g" \
    -e "s/rho/\\\rho/g" \
    -e "s/sigma/\\\sigma/g" \
    -e "s/omega/\\\omega/g" \
    -e "s/zeta/\\\zeta/g" \
    -e "s/sin/\\\sin/g" \
    -e "s/cos/\\\cos/g" \
    -e "s/23/{23}/g" \
    -e "s/+-/ -/g" |\
cat -s  | \
sed -e "s/ $//g" | \
tr  -s '\012' '\012' \
>junk

mv junk $1

#! /bin/sh

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

# Bourne shell script: lbl2sympar_txt2r

# Label file to symbolic parameters conversion
# P.J.Gawthrop August 1996
# Copyright (c) P.J.Gawthrop, 1996.

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


# Inform user
echo "Creating $1_sympar.r"

rm -f mtt_error

# This is the main transformation using gawk
gawk -f $MTTPATH/lbl2sympar.awk $1_lbl.txt > $1_sympar.r 2>mtt_error

err_length=$(wc -c <mtt_error)

# Test for errors and print if any
if [ $err_length != "0" ]
then
  echo lbl2sympar_txt2 $1 failed '...'
  cat mtt_error
  exit 1
else
  exit 0
fi

function [port_bonds, status] = abg2cbg(system_name, ...
    system_type, full_name, ...
    port_bonds,infofile)
% [bonds,status] = abg2cbg(system_name, system_type, full_name, port_bonds, infofile)
% 
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%     %%%%% Model Transformation Tools %%%%%
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% Matlab function  abg2cbg.m
% Acausal bond graph to causal bond graph: mfile format

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %%Revision 1.7  1996/08/16  12:58:58  peter
% %% Now does preferred causality of I and C.
% %% Revision 1.6  1996/08/09 08:27:29  peter
% %% Added a few deguging lines
% %%
% %% Revision 1.5  1996/08/08 18:06:18  peter
% %% Unified file naming scheme
% %%
% %% Revision 1.4  1996/08/08 08:30:06  peter
% %% The cbg filename contains the system name - this makes things easier
% %% when setting up the m to fig translation and m to ese translation
% %%
% %% Revision 1.3  1996/08/05 18:53:21  peter
% %% Fixed bug passing causality from subsystems.
% %%
% %% Revision 1.2  1996/08/05 15:41:41  peter
% %% Now recursively does causality on subsystems.
% %%
% %% Revision 1.1  1996/08/04 17:55:55  peter
% %% Initial revision
% %%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


pc = '%';
if nargin<1
  system_name = 'no_name';
end;

if nargin<4
  port_bonds = [];
end;

if nargin<5
  infofile = 'stdout';
end;

% Create the (full) system name
if length(full_name)==0
  full_name = system_name;
  system_type = system_name;
else
  full_name = [full_name, '_', system_name];
end;

fun_name = [system_type, '_abg']

% If no such function - then there is nothing to be done.
if exist(fun_name)~=2
  mtt_info(['m-file ', fun_name, ' does not exist'], infofile);
  bonds = [];
  status = [];
  return
end;

% Evaluate the system function to get the bonds
eval(['[bonds,components,n_ports]=', fun_name, ';']);

% Find number of bonds
[n_bonds,columns] = size(bonds);
if (columns ~= 2)&(n_bonds>0)
  error('Incorrect bonds matrix: must have 2 columns');
end;

 
% Find number of components
[n_components,columns] = size(components);
if n_components==0 % there is nothing to be done
  return
end;

% Find number of port bonds
[n_port_bonds,columns] = size(port_bonds);

% Check compatibility - if ok copy port bonds to the internal bonds list.
if n_port_bonds~=n_ports
  mtt_info(sprintf('%1.0f port bonds incompatible with %1.0f ports', ...
      n_port_bonds, n_ports), infofile);
else % Copy the port bonds
  for i = 1:n_ports
    bonds(i,:) = port_bonds(i,:);
  end;
end;


% Set initial status
status = -ones(n_components,1);
total = 2*n_bonds;

done = sum(sum(abs(bonds)))/total*100;

% Outer while loop sets preferred causality
ci_index=1;
while ci_index>0
  old_done = inf;
  % Inner loop propogates causality
  while done~=old_done
    disp(sprintf('Causality is %3.0f%s complete.', done, pc));
    old_done = done;
  
    for i = 1:n_components
      comp = nozeros(components(i,:));
      bond_list = abs(comp);
      direction = sign(comp)'*[1 1];
      % Convert from arrow orientated to component orientated causality
      comp_bonds = bonds(bond_list,:).*direction;
      eval([ '[comp_type,name,cr,arg] = ', system_type, '_cmp(i);' ]);
      % change name of 0 and 1 components -- matlab doesn't like numbers here
      if strcmp(comp_type,'0')
	comp_type = 'zero';
      end;
      if strcmp(comp_type,'1')
	comp_type = 'one';
      end;
      
      % Component causality procedure name
      cause_name = [comp_type, '_cause'];
      
      % Invoke  the appropriate causality procedure
      if exist(cause_name)~=2 % Try a compound component
	% disp('------------PUSH-----------------');
	[comp_bonds,s] = abg2cbg(name, comp_type, full_name, comp_bonds, ...
	    infofile);
	status(i)=max(abs(s));
	% disp('------------POP-----------------');
      else % its a simple component
	% disp(['---', name, ' (', cause_name, ') ---']);
	% comp_bonds
	eval([ '[comp_bonds,status(i)] = ', cause_name, '(comp_bonds);' ]);
	% comp_bonds
      end;
      
      % Update the full bonds list
      % and convert from component orientated to arrow orientated causality
      bonds(bond_list,:) = comp_bonds.*direction; 
    end;
    
    
    done = sum(sum(abs(bonds)))/total*100;
    %  mtt_info(sprintf('Causality is %3.0f%s complete.', done, pc), infofile);
    
  end;
  
  % Set causality of C and I which is not set
  [ci_index,prefered] = getdynamic(status,system_type);
  if ci_index>0
    bond_index = components(ci_index,1) % its a one port
    bonds(bond_index,1) = prefered;
    bonds(bond_index,2) = prefered;
  end;
  
end;

% Print final causality
final_done =  (sum(status==zeros(n_components,1))/n_components)*100;;
mtt_info(sprintf('Final causality is %3.0f%s complete.', final_done, pc), infofile);

% List overcausal bonds
[over_causal_bonds,n] = getindex(status,1);
if n>0
  for i=over_causal_bonds'
    eval([ '[comp_type,name] = ', system_type, '_cmp(i);' ]);
    mtt_info(sprintf('Component %s (%s) is overcausal', name, comp_type), ...
      infofile);
  end;
end;

% List undercausal bonds
[under_causal_bonds,n] = getindex(status,-1);
if n>0
  for i=under_causal_bonds'
    eval([ '[comp_type,name] = ', system_type, '_cmp(i);' ]);
    mtt_info(sprintf('Component %s (%s) is undercausal', name, comp_type), ...
      infofile);
  end;
end;

file_name = [full_name, '_', system_type]
write_cbg(file_name,bonds,status);

% Return the port bonds
for i = 1:n_ports
  port_bonds(i,:) = bonds(i,:);
end;

disp('----------------------');

function [index,distance] = adjbond(point,arrow_end,other_end);
% adjbond:  Determines the bond closest to the point 
%  [index,distance] = adjbond(point,arrow_end,other_end);


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


[n,m] = size(other_end);
one = ones(n,1);
arrow_distance = length2d(one*point - arrow_end);
min_arrow_distance = min(arrow_distance);

other_distance = length2d(one*point - other_end);
min_other_distance = min(other_distance);

min_distance = min([arrow_distance; other_distance]);

adjacent = [arrow_distance, other_distance] == min_distance*[one one];
[index,n] = getindex(adjacent,1);
distance = min_distance;

function indices = adjcomp(arrow_end,other_end,components);
% adjcomp: Determines the two components at each end of the bond



% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


[n,m] = size(components);
one = ones(n,1);

arrow_distance = length2d(one*arrow_end - components(:,1:2));
min_arrow_distance = min(arrow_distance);

other_distance = length2d(one*other_end - components(:,1:2));
min_other_distance = min(other_distance);

arrow_adjacent = arrow_distance==min_arrow_distance*one;
other_adjacent = other_distance==min_other_distance*one;
[index,n] = getindex([arrow_adjacent,other_adjacent],1);

if index(1,2)==1
  indices = index(1:2,1)';
else
  indices = index(2:-1:1,1)';
end;

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$



% %% Revision 1.1  1996/08/27  12:50:43  peterg
% %% Initial revision
% %%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

arg = '';

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;

function structure = cbg2ese(system_name, system_type, full_name, ...
    repetition,...
    structure, infofile)
% 
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%     %%%%% Model Transformation Tools %%%%%
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% Matlab function  cbg2ese.m
% Acausal bond graph to causal bond graph: mfile format
% Structure matrix [states,nonstates,inputs,outputs,zero_outputs]

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %% Revision 1.5  1996/08/24  15:02:23  peter
% %% Writes `END;' to keep reduce happy.
% %%
% %% Revision 1.4  1996/08/19 09:03:41  peter
% %% Handles repeating components.
% %%
% %% Revision 1.3  1996/08/18 20:08:02  peter
% %% Included additional structure: structure(5) = zero_outputs.
% %%
% %% Revision 1.2  1996/08/08 18:08:11  peter
% %% Sorted out file naming sceme
% %%
% %% Revision 1.1  1996/08/08 15:53:23  peter
% %% Initial revision
% %%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

system_name, system_type, full_name,

pc = '%';
if nargin<3
  eqnfile = 'stdout';
end;
if nargin<4
  infofile = 'stdout';
end;

% Create the (full) system name
if length(full_name)==0
  full_name = system_name;
  system_type = system_name;
else
  full_name = [full_name, '_', system_name];
end;

full_name_type = [full_name, '_', system_type];
cbg_name = [full_name_type, '_cbg'];
abg_name = [system_type, '_abg'];
cmp_name = [system_type, '_cmp'];
% Return if cbg file doesn't exist
if exist(cbg_name)~=2
  return
end;

% Setup file
ese_name = sprintf('%s_%1.0f_ese.r', full_name, repetition);
filenum = fopen(ese_name, 'w');

% Evaluate the system function to get the bonds
  eval(['[junk,components]=', abg_name, ';']);
  eval(['bonds=', cbg_name, ';']);
  % Find number of bonds
  [n_bonds,columns] = size(bonds);
  if (columns ~= 2)&(n_bonds>0)
    error('Incorrect bonds matrix: must have 2 columns');
  end;

  % Find number of components
  [n_components,columns] = size(components);
  n_components = n_components
  for i = 1:n_components
    comp = nozeros(components(i,:));
    bond_list = abs(comp);
    direction = sign(comp)'*[1 1];
    % Convert from arrow orientated to component orientated causality
    comp_bonds = bonds(bond_list,:).*direction;
    eval([ '[comp_type,comp_name,cr,args,repetitions] = ', cmp_name, '(i);' ]);
    % change name of 0 and 1 components -- matlab doesn't like numbers here
    if strcmp(comp_type,'0')
      comp_type = 'zero';
    end;
    if strcmp(comp_type,'1')
      comp_type = 'one';
    end;
    
    comp_type = comp_type
    ports = length(bond_list)
    repetitions = repetitions
    
    if repetitions>1
      port_pairs = ports/2;
      if round(port_pairs)~=port_pairs;
	mtt_info(['Repeated component ', comp_name, ...
	      ' has an odd number of ports - ignoring repetitions']);
	repetitions = 1;
      end;
    end;
    
    if repetitions>1
      odd_bonds = bond_list(1:2:ports-1);
      even_bonds = bond_list(2:2:ports);
      next_bond = max(max(abs(components)))+1;
    end;
    
    for k = 1:repetitions
      
      if repetitions>1
	if k==1
	  bond_list(1:2:ports-1) = odd_bonds;
	else
	  bond_list(1:2:ports-1) = bond_list(2:2:ports);
	end;
	
	if k==repetitions
	  bond_list(2:2:ports) = even_bonds;
	else
	  new_bonds = [next_bond:next_bond+port_pairs-1];
	  next_bond = next_bond+port_pairs;
	  bond_list(2:2:ports) = new_bonds;
	end;
      end;
	
      % Invoke the appropriate equation-generating procedure
      name_r = sprintf('%s_%1.0f', full_name, repetition);
      eqn_name = [comp_type, '_eqn']
      if exist(eqn_name)~=2 % Try a compound component
	disp('---PUSH---');
	structure = cbg2ese(comp_name, comp_type, full_name, k, ...
	    structure,  infofile);
	% Link up the bonds
	name_k = sprintf('%s_%1.0f', full_name, k);
	name_comp_name = sprintf('%s_%s_%1.0f', full_name, comp_name, k);
	for port_number=1:length(bond_list)

	  % Effort
	  if comp_bonds(port_number,1)==1 % Source
	     fprintf(filenum, '%s_MTT_inport_%1.0f := %s;\n', ...
		 name_comp_name, port_number, varname(name_r, ...
		 bond_list(port_number),1));
	   else % sensor
	     	 fprintf(filenum, '%s := %s_MTT_outport_%1.0f;\n', ...
		 varname(name_r, ...
		 bond_list(port_number),1), name_comp_name, port_number);
	  end;
	  % flow
	  if comp_bonds(port_number,2)==-1 % Source
	     fprintf(filenum, '%s_MTT_inport_%1.0f := %s;\n', ...
		 name_comp_name, port_number, varname(name_r, ...
		 bond_list(port_number),-1));
	   else % sensor
	     	 fprintf(filenum, '%s := %s_MTT_outport_%1.0f;\n', ...
		 varname(name_r, ...
		 bond_list(port_number),-1), name_comp_name, port_number);
	   end;	
	 end;
	
	disp('---POP---');
      else % its a simple component
	eval(['structure = ', ...
	      eqn_name, ...
	      '(name_r,bond_list,comp_bonds,direction,cr,args,structure,filenum);' ]);
      end;
    end;
  end;

  fclose(filenum);

function cbg2fig(system_name, ...
                 system_type, full_name, ...
                 stroke_length, stroke_thickness, stroke_colour, ...
                 comp_font, comp_colour_u, comp_colour_o)
 
% $$$ cbg2fig(system_name, ...
% $$$                  system_type, full_name, ...
% $$$                  stroke_length, stroke_thickness, stroke_colour, ...
% $$$                  comp_font, comp_colour_u, comp_colour_o)
 
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%     %%%%% Model Transformation Tools %%%%%
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Matlab function  cbg_m2fig
% Causal bond graph:  mfile format to fig file format
% The resultant fig file is the original _abg.fig with
% additional causal strokes superimposed.
%
% P.J.Gawthrop May 1996
% Copyright (c) P.J.Gawthrop, 1996.

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %% Revision 1.2  1996/08/05 20:15:39  peter
% %% Prepared for recursive version.
% %%
% %% Revision 1.1  1996/08/05 18:12:25  peter
% %% Initial revision
% %%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


if nargin<4
  stroke_length = 20;
end;

if nargin<5
  stroke_thickness = 2;
end;

if nargin<6
  stroke_colour = 1; %Blue
end;

if nargin<7
  comp_font = 18; %Helvetica bold
end;

if nargin<8
  comp_colour_u = 12; %Green
end;

if nargin<9
  comp_colour_o = 4; %Red
end;

% Create the (full) system name
if length(full_name)==0
  full_name = system_name;
  system_type = system_name;
else
  full_name = [full_name, '_', system_name];
end;

full_name_type = [full_name, '_', system_type];
fig_name = [full_name_type, '_cbg.fig'];

% Return if cbg file doesn't exist
if exist(fig_name)~=2
  return
end;

% Setup file - append to the fig file
filenum = fopen(fig_name, 'a');

% Get the raw and the processed bonds
eval(['[rbonds,rstrokes,rcomponents] = ', system_type, '_rbg;']);
eval(['[bonds] = ', system_type, '_abg;']);

% Get the causal bonds
eval(['[cbonds,status]=', full_name_type, '_cbg;']);

% Check sizes
[N_components,Columns] = size(rcomponents);
if (Columns ~= 13)
  error('Incorrect rcomponents matrix: must have 13 columns');
end;
M_components = Columns;

% Rotation matrix
rot = [0 -1; 1 0];


% Determine coordinates of the arrow end of the bond and the other end
% and other geometry
other_end_1 = rbonds(:,1:2);
arrow_end = rbonds(:,3:4);
other_end_2 = rbonds(:,5:6);

distance_1 = length2d(other_end_1 - arrow_end);
distance_2 = length2d(other_end_2 - arrow_end);
which_end = (distance_1>distance_2)*[1 1];
one = ones(size(which_end));
other_end = which_end.*other_end_1 + (one-which_end).*other_end_2;
arrow_barb = which_end.*other_end_2 + (one-which_end).*other_end_1;
arrow_vector =  arrow_barb-arrow_end;
unit_arrow_vector = arrow_vector./(length2d(arrow_vector)*[1 1]);
bond_vector = (arrow_end - other_end);
unit_bond_vector = bond_vector./(length2d(bond_vector)*[1 1]);
unit_stroke_vector = (rot*unit_bond_vector')';
  
% Get indices of bonds with changed causality
changed_e = bonds(:,1)~=cbonds(:,1);
changed_f = bonds(:,2)~=cbonds(:,2);
changed = changed_e|changed_f;
index_e  = getindex(changed_e,1)';
index_f  = getindex(changed_f,1)';
index  = getindex(changed,1)';

% Print the new strokes in fig format
if index(1,1)>0
  for i = index_e % Do the effort stroke - opp. side to arrow

    if cbonds(i,1)==1 % Stroke at arrow end
      stroke_end_1 = arrow_end(i,:);
    else
      stroke_end_1 = other_end(i,:);
    end;
    
    sig = sign(unit_arrow_vector(i,:)*unit_stroke_vector(i,:)');
    stroke_end_2 = stroke_end_1 - stroke_length*sig*unit_stroke_vector(i,:);
    

    %print the fig3 format firstline spec.
    polyline = 2; 
    firstline = fig3(polyline,stroke_thickness,stroke_colour);
    fprintf(filenum, '%s\n', firstline);

    fprintf(filenum, '	%4.0f %4.0f %4.0f %4.0f \n', ...
	stroke_end_1(1), stroke_end_1(2), ...
	stroke_end_2(1), stroke_end_2(2) );
  end;

  for i = index_f % Do the flow stroke - same side as arrow

    if cbonds(i,2)==1 % Stroke at arrow end
      stroke_end_1 = arrow_end(i,:);
    else
      stroke_end_1 = other_end(i,:);
    end;
    
    sig = sign(unit_arrow_vector(i,:)*unit_stroke_vector(i,:)');
    stroke_end_2 = stroke_end_1 + stroke_length*sig*unit_stroke_vector(i,:);
    

    %print the fig3 format firstline spec.
    polyline = 2; 
    firstline = fig3(polyline,stroke_thickness,stroke_colour);
    fprintf(filenum, '%s\n', firstline);

    fprintf(filenum, '	%4.0f %4.0f %4.0f %4.0f \n', ...
	stroke_end_1(1), stroke_end_1(2), ...
	stroke_end_2(1), stroke_end_2(2) );
  end;
end;

% Print all the components - coloured acording to causality.
for i = 1:N_components
  fig_params = rcomponents(i,3:M_components);
  coords = rcomponents(i,1:2);
  
  if status(i)==-1  %Then under causal
    fig_params(3) = comp_colour_u;
    fig_params(6) = comp_font;
  end;

  if status(i)==1  %Then over causal
    fig_params(3) = comp_colour_o;
    fig_params(6) = comp_font;
  end;


  %Now print the component in fig format
  eval(['[comp_type,comp_name] = ', system_type, '_cmp(i);']);

  Terminator = '\\001';   
  for j = 1:length(fig_params)
    fprintf(filenum, '%1.0f ', fig_params(j));
  end;
  fprintf(filenum, '%1.0f %1.0f ', coords(1), coords(2)); 
  fprintf(filenum, '%s:%s%s\n', comp_type, comp_name, Terminator);
  
  % If it's a subsystem, do the fig file for that as well
  cbg2fig(comp_name, ...
                 comp_type, full_name, ...
                 stroke_length, stroke_thickness, stroke_colour, ...
                 comp_font, comp_colour_u, comp_colour_o);
end;

% Close the file
fclose(filenum);
return

function name = cr2name(RHS_number,LHS_cause,RHS_cause,cr,args);
% cr2name. Constructs a string for the cr of a component.

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %% Revision 1.1  1996/08/22 13:14:39  peter
% %% Initial revision
% %%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


if length(args)==0
  a_comma = '';
else
  a_comma = ',';
end;

if length(cr)==0
  cause_name = '';
  lp = '';
  rp = '';
  c_comma = '';
else
  cause_name = cause2name(LHS_cause);
  lp = '(';
  rp = ')';
  c_comma = ',';
end;

name = sprintf('%s%s%s%s%s%s%s%s', ...
    cr,lp, args, a_comma, cause_name, c_comma, ...
    varname(RHS_number,RHS_cause),rp);
    
    

function [A,B,C,D] = csm2sm(A,B,C,D,E);
% [A,B,C,D] = csm2sm(A,B,C,D,E);
% Converts constrained-state matrices to state matrices.

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


A = E\A;
B = E\B;

function  fr = dm2fr(A,B,C,D,E,W,u0)
% fr = dm2fr(A,B,C,D,E,W,u0)
% Descriptor matrix to frequency response.
% A,B,C,D,E - descriptor matrices
% W vector of frequency points
% u0 input gain vector: u = u0*unit phasor

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %% Revision 1.5  1996/08/16 14:26:37  peter
% %% Check and fix size of u0.
% %%
% %% Revision 1.4  1996/08/15 12:50:51  peter
% %% Put in a conj to undo effect of transpose.
% %%
% %% Revision 1.3  1996/08/15 11:53:44  peter
% %% Now has u0 input vector
% %%
% %% Revision 1.2  1996/08/15 10:24:28  peter
% %% Includes u0 argument.
% %%
% %% Revision 1.1  1996/08/10 14:11:28  peter
% %% Initial revision
% %%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[Ny,Nu] = size(D);
[Ny,Nx] = size(C);
N = length(W);

if nargin<7
  U0 = zeros(Nu,1);
  U0(1) = 1;
else
  for i=1:Nu
    U0(i) = u0(i);
  end;
end;

u0 = U0;

[n,m]=size(W);
if m>n
  W=W';
end;

[n,m]=size(u0);
if m>n
  u0=u0';
end;


fr = zeros(N,Ny);
i = 0;
for w = W'
  i = i+1;
  FR = C*( (E*j*w - A) \ B*u0 ) + D*u0;
  fr(i,:) = conj(FR');
end;

function ir = dm2ir(A,B,C,D,E,T);
% ir = dm2ir(A,B,C,D,E,T);
% Descriptor matrix to impulse response.
% NB At the moment - this assumes that E is unity .....
% A,B,C,D,E - descriptor matrices
% T vector of time points

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %% Revision 1.3  1996/08/10 14:16:04  peter
% %% Now has T vector input - it was tmax & tstep
% %%
% %% Revision 1.2  1996/08/10 11:49:39  peter
% %% Fixed multi-input/output problem
% %%
% %% Revision 1.1  1996/08/10 10:26:00  peter
% %% Initial revision
% %%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


[N,M] = size(T);
if M>N
  T = T';
  N = M;
end;

[Ny,Nu] = size(D);
NN=Ny*Nu;

ir = zeros(N,NN);
i = 0;
for t = T'
  i=i+1;
  ir(i,:) = reshape(C*expm(A*t)*B, 1,NN);
end;

function [Y,X] = dm2sr(A,B,C,D,E,T,u0,x0);
% [Y,X] = dm2sr(A,B,C,D,E,T,u0,x0);
% Descriptor matrix to impulse response.
% NB At the moment - this assumes that E is unity .....
% A,B,C,D,E - descriptor matrices
% T vector of time points
% u0 input gain vector: u = u0*unit step.

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %% Revision 1.4  1996/08/15 08:34:08  peter
% %% Added step gain (u0) and initial condition (x0)
% %%
% %% Revision 1.3  1996/08/11 19:33:24  peter
% %% Replaced exp by expm - whoops!
% %%
% %% Revision 1.2  1996/08/11 10:37:40  peter
% %% Corrected mistake in step-response calculation.
% %%
% %% Revision 1.1  1996/08/11 09:42:40  peter
% %% Initial revision
% %%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[Ny,Nu] = size(D);
[Ny,Nx] = size(C);

if nargin<7
  u0 = zeros(Nu,1);
  u0(1) = 1;
end;

if nargin<8
  x0 = zeros(Nx,1);
end;

[N,M] = size(T);
if M>N
  T = T';
  N = M;
end;

one = eye(Nx);

Y = zeros(N,Ny);
X = zeros(N,Nx);
i = 0;
for t = T'
  i=i+1;
  if Nx>0
    expAt = expm(A*t);
    x = ( A\(expAt-one) )*B*u0 + expAt*x0;
    X(i,:) = x';
    if Ny>0
      y = C*x + D*u0;
      Y(i,:) = y';
    end;
  elseif Ny>0
    y = D*u0;
    Y(i,:) = y';
  end;
end;

function [index,prefered] = getdynamic(status,system_type);
% Get the index of a dynamic components which is not set.

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %% Revision 1.1  1996/08/16 12:50:41  peter
% %% Initial revision
% %%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

index=0; prefered=0;
n_components = length(status);

for i = 1:n_components
  if status(i)==-1 % Undercausal
    eval([ '[comp_type,name,cr,arg] = ', system_type, '_cmp(i);' ]); 
    if strcmp(comp_type,'C')
      index=i;
      prefered=-1;
      break;
    end;
    if strcmp(comp_type,'I')
      index=i;
      prefered=1;
      break;
    end;
  end;
end;

function [index, n, otherindex, m] = getindex(array,value);
% Finds the n indices of the elements of array equal to value
% otherindex contains indeces of the the m other  elements.


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[N,M] = size(array);
n=0; m=0;
for i=1:N
  for j = 1:M
    if array(i,j)==value
      n=n+1;
      index(n,:) = [i j];
    else
      m=m+1;
      otherindex(m,:) = [i j];
    end;
  end;
end;


if (M==1)&(n>0)
  index = index(:,1);
end;

if (M==1)&(m>0)
  otherindex = otherindex(:,1);
end;

% Octave doesn't like empty matrices
if n==0
  index=0;
end;

if m==0
  otherindex=0;
end;

function makedef(structure,deffile);

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %% Revision 1.2  1996/08/18 20:05:20  peter
% %% Put unded version control
% %%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



filenum = fopen(deffile,'w');

states = structure(1);
nonstates=structure(2);
inputs=structure(3);
outputs=structure(4);
zero_outputs = structure(5);

pc = '%';
% Declare reduce constants;
fprintf(filenum, 'MTTNx := %1.0f;\n', states);
fprintf(filenum, 'MTTNz := %1.0f;\n', nonstates);
fprintf(filenum, 'MTTNu := %1.0f;\n', inputs);
fprintf(filenum, 'MTTNy := %1.0f;\n', outputs);
fprintf(filenum, 'MTTNyz := %1.0f;\n', zero_outputs);

% Declare reduce matrices
fprintf(filenum, '%s Declare reduce matrices\n', pc);
if states>0
  fprintf(filenum, 'matrix MTTx(%1.0f,1);\n', states);
  fprintf(filenum, 'matrix MTTdx(%1.0f,1);\n', states);
end;
if nonstates>0
  fprintf(filenum, 'matrix MTTz(%1.0f,1);\n', nonstates);
  fprintf(filenum, 'matrix MTTdz(%1.0f,1);\n', nonstates);
end;
if inputs>0
  fprintf(filenum, 'matrix MTTu(%1.0f,1);\n', inputs);
end;
if outputs>0
  fprintf(filenum, 'matrix MTTy(%1.0f,1);\n', outputs);
end;
if zero_outputs>0
  fprintf(filenum, 'matrix MTTyz(%1.0f,1);\n', zero_outputs);
  fprintf(filenum, 'matrix MTTui(%1.0f,1);\n', zero_outputs);
end;

% Make an Nx x Nx unit matrix
if states>0
  fprintf(filenum, 'matrix MTTI(%1.0f,%1.0f);\n', states,states);
  for i = 1:states
    fprintf(filenum, 'MTTI(%1.0f,%1.0f) := 1;\n', i, i);
  end
end;

% Set the y, yz, u, x and dx matrices
fprintf(filenum, '%s Set the y, yz, u and x matrices\n', pc);
for i=1:outputs
  fprintf(filenum, 'MTTy(%1.0f,1) := MTTy%1.0f;\n', i, i);
end;
for i=1:zero_outputs
  fprintf(filenum, 'MTTyz(%1.0f,1) := MTTyz%1.0f;\n', i, i);
  fprintf(filenum, 'MTTui(%1.0f,1) := MTTui%1.0f;\n', i, i);
end;
for i=1:inputs
  fprintf(filenum, 'MTTu(%1.0f,1) := MTTu%1.0f;\n', i, i);
end;
for i=1:states
  fprintf(filenum, 'MTTx(%1.0f,1) := MTTx%1.0f;\n', i, i);
end;
for i=1:nonstates
  fprintf(filenum, 'MTTdz(%1.0f,1) := MTTdz%1.0f;\n', i, i);
end;

fprintf(filenum, 'END;');
fclose(filenum);
  

function mtt_info(info, infofile);
% mtt_info(info, infofile);

% function [bonds, status] = abg2cbg(bonds,components,system_name,filename)
% 
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%     %%%%% Model Transformation Tools %%%%%
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



if nargin<2
  infofile = 'mtt_info.txt';
end;

fprintf(infofile, 'INFORMATION: %s\n', info);
% fflush(infofile);

function [bonds,components] = rbg2abg(rbonds,rstrokes,rcomponents,rports,infofile)

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %% Revision 1.5  1996/08/24 19:21:26  peter
% %% More specific error messages.
% %%
% %% Revision 1.4  1996/08/24 18:00:33  peter
% %% Fixed bug with finding ports.
% %%
% %% Revision 1.3  1996/08/09 08:26:35  peter
% %% Cosmetic tidy up.
% %%
% %% Revision 1.2  1996/08/04 18:37:57  peter
% %% Fixed  no causal strokes bug.
% %%
% %% Revision 1.1  1996/08/04 18:30:14  peter
% %% Initial revision
% %%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


if nargin<5
  infofile='stdout';
else
  fnum = fopen(infofile, 'w');
end;

% Xfig scaling factor
scale = 1200.0;

% Rotation matrix
rot = [0 -1; 1 0];

% Find number of strokes
[n_strokes,columns] = size(rstrokes);
if (columns ~= 4)&(n_strokes>0)
  error('Incorrect rstrokes matrix: must have 4 columns');
end;

% Find number of bonds
[n_bonds,columns] = size(rbonds);
if (columns ~= 6)&(n_bonds>0)
  error('Incorrect rbonds matrix: must have 6 columns');
end;

% Find number of components
[n_components,columns] = size(rcomponents);

% Find number of ports
[n_ports,columns] = size(rports);

% Determine coordinates of the arrow end of the bond and the other end
other_end_1 = rbonds(:,1:2);
arrow_end   = rbonds(:,3:4);
other_end_2 = rbonds(:,5:6);

distance_1   = length2d(other_end_1 - arrow_end);
distance_2   = length2d(other_end_2 - arrow_end);
which_end    = (distance_1>distance_2)*[1 1];
one          = ones(size(which_end));
other_end    = which_end.*other_end_1 + (one-which_end).*other_end_2;
arrow_vector = ( which_end.*other_end_2 + (one-which_end).*other_end_1 ) - ...
    arrow_end;
 
% Locate the bond end nearest to each port
for i = 1:n_ports
  near_bond = adjbond(rports(i,1:2),arrow_end,other_end);
  port_near_bond(i,:) = [near_bond, rports(i,3)];
end;

% Locate the components at the ends of each bond
for i = 1:n_bonds
  comp_near_bond(i,:) = adjcomp(arrow_end(i,:),other_end(i,:),rcomponents);
end;

% Produce a list of bonds on each component - sorted if explicit port numbers
for i = 1:n_components    
  [index,n] = getindex(comp_near_bond,i);
    % Error message in case we need it!
  port_error = sprintf(... 
      'Component at (%1.3f,%1.3f) has inconsistent port numbers', ...
      rcomponents(i,1)/scale, rcomponents(i,2)/scale);

  if index(1,1) ~= 0 % Then its a genuine component 
    one = ones(n,1);
    bond_list = index(:,1); %  bond at component
    bond_end = index(:,2);  % which end of bond at component?


    % Default sort of bonds (ie no change)
    sort_index = [1:n]'; 
    
    if n_ports>0
      % Are the component ports numbered? (either they all are or none are)
      k=0;
      for j = 1:n
	b = bond_list(j); e = bond_end(j);
	% Find all ports on this bond
	[port_indices,m] = getindex(port_near_bond(:,1),b);
	% Now find the one at this end  - if any
	port_index = 0;
	for l=1:m
	  if port_near_bond(port_indices(l),2)==e
	    port_index = port_indices(l);
	    break;
	  end;
	end;
	if port_index>0
	  % and put the corresponding number in the list
	  k=k+1;
	  port_number(k,1) = port_near_bond(port_index,3);
	end;
      end;

	% Must have a lable for each port or non at all
      if k==n
	[sorted,sort_index]=sort(port_number);
	if sum(sorted==[1:n]')~=n % The there is something wrong
	  mtt_info(port_error,fnum);
	  mtt_info(sprintf('\t it must have ports from 1:%1.0f\n', n), fnum);
	end;
      elseif k~=0
	mtt_info(port_error,fnum);
	mtt_info(sprintf('\t it must have 0 or %1.0f ports\n', n), fnum);
      end;
 
    end;
  end;

  % direction is 1 if arrow at component else -1
  direction = -sign(index(:,2)-1.5*one);
  signed_bond_list = bond_list.*direction;
  
  % Write out bond list sorted by port number (if any)
  for j = 1:length(sort_index)
    jj = sort_index(j);
    components(i,j) = signed_bond_list(jj);
  end;
end;

% Deduce causality from the strokes (if any).
causality = zeros(n_bonds,2);
if n_strokes>0
  % Find out location of centre and ends of stroke.
  stroke_end_1 = [rstrokes(:,1) rstrokes(:,2)];
  stroke_end_2 = [rstrokes(:,3) rstrokes(:,4)];
  
  stroke_centre = (stroke_end_1 + stroke_end_2)/2;
  stroke_vector = (stroke_end_1 - stroke_end_2);
  stroke_length = length2d(stroke_vector);

% Deduce bond causality from the strokes 
  for i = 1:n_strokes
    stroke = [stroke_centre(i,:) 
    stroke_end_1(i,:)
    stroke_end_2(i,:)];


    % Find the nearest bond end.
    [index,distance] = adjbond(stroke(1,:),arrow_end,other_end);
    if (distance>2*stroke_length(i))
      info = sprintf('Stroke at (%4.3f,%4.3f) is %4.3f away from the nearest bond\n', ...
	stroke(1,1)/scale, stroke(1,2)/scale, distance/scale);
    end;
  
    % Bond end coordinates
    j = index(1,1);
    which_end = index(1,2)==1;
    bond_end = arrow_end(j,:)*which_end + other_end(j,:)*(1-which_end);
  
    % Now decide which bit of the stroke is nearest
    stroke_index = adjbond(bond_end,stroke,zeros(size(stroke)));
  
    if stroke_index(1)==1 % uni-causal stroke
      causality(j,1:2) = (2*which_end-1)*[1 1];
    else % bicausal stroke
      % Find out whether stroke is on flow side of bond
      stroke_direction = stroke(1,:) - stroke(stroke_index(1),:);
      flow_side = stroke_direction*arrow_vector(j,:)'>0;
      causality(j,1+flow_side) = 2*which_end-1;
    end;
  end;
end;

bonds = causality;
fclose(fnum);

function [Y,X] = sm2sr(A,B,C,D,T,u0,x0);
% [Y,X] = sm2sr(A,B,C,D,T,u0,x0);
% Constrained-state matrix to impulse response.
% A,B,C,D,E - (constrained) state matrices
% T vector of time points
% u0 input gain vector: u = u0*unit step.

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


[Ny,Nu] = size(D);
[Ny,Nx] = size(C);

if nargin<7
  u0 = zeros(Nu,1);
  u0(1) = 1;
end;

if nargin<8
  x0 = zeros(Nx,1);
end;

[N,M] = size(T);
if M>N
  T = T';
  N = M;
end;



one = eye(Nx);

Y = zeros(N,Ny);
X = zeros(N,Nx);
i = 0;
for t = T'
  i=i+1;
  if Nx>0
    expAt = expm(A*t);
    x = ( A\(expAt-one) )*B*u0 + expAt*x0;
    X(i,:) = x';
    if Ny>0
      y = C*x + D*u0;
      Y(i,:) = y';
    end;
  elseif Ny>0
    y = D*u0;
    Y(i,:) = y';
  end;
end;

function ch = str2ch(str,n)
% str2ch - pulls the nth char from a string
% 
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%     %%%%% Model Transformation Tools %%%%%
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% Matlab function  str2ch.m
% ch = str2ch(str,n)
% Copyright (c) P.J. Gawthrop, 1996.


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

implicit_str_to_num_ok = "true";
str = abs(str);

if (n>0) &  (n<=length(str))
  ch = setstr(str(n));
else
  ch = '';
end;

function name = varname(name,index,causality);
% varname - Creates name of bond graph variable

% 
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%     %%%%% Model Transformation Tools %%%%%
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% Matlab function  varname.m
% Acausal bond graph to causal bond graph: mfile format


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


bond_name = [name,'_bond'];

name =sprintf('%s%1.0f_%s', bond_name, index, cause2name(causality));

function write_matrix(matrix,name);

% Writes the matrix function file

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %% Revision 1.3  1996/08/14 19:20:41  peter
% %% Fixed output naming bug.
% %%
% %% Revision 1.2  1996/08/14 08:36:52  peter
% %% Puts a tab between columns.
% %%
% %% Revision 1.1  1996/08/14 08:21:27  peter
% %% Initial revision
% %%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


filename = [name, '.m'];
filenum = fopen(filename,'w');

% Write the function m-file for the causal bond graph
pc = '%';
fprintf(filenum, 'function data = %s\n', name);
fprintf(filenum, '%s data = %s\n\n', pc, name);

fprintf(filename, 'data = [\n');

[N,M] = size(matrix);
for row = 1:N
  for col = 1:M
    re = real(matrix(row,col));
    im = imag(matrix(row,col));
    fprintf(filename, '%g', re);
    if im ~= 0
      fprintf(filename, '+ %g*i', im);
    end
    if col<M
      fprintf(filename, '\t');
    end
  end;
  fprintf(filename, '\n');
end;

fprintf(filename, '];\n');
fprintf(filename, '\n');

fclose(filenum);

#! /bin/sh

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

# Bourne shell script: m2dat

# Transforms data in .m file to just the data - ie strips the matlab stuff

# Copyright (c) P.J.Gawthrop, 1996.

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

echo Creating $1.dat

awk '$1 ~ /[0-9]/ && NF>0 {print} ' <$1.m  >$1.dat

#! /bin/sh

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

# Bourne shell script: mtt

# Usage: mtt system_name  representation_name
# E.g.: 
# 	mtt dc mtf
# converts dc to matlab transfer function form.
#
# P.J.Gawthrop Oct 1989, Dec 1989, 1990, 1991, Oct 1993, Dec 1993, Jan 1994.
# Copyright (c) P.J.Gawthrop, 1989, 1990, 1991, 1993, 1994.

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

#Inform user
echo Creating $1_cr.r

cat - > $1_cr.r << EOF
%Default CR file
%File $1_cr.r
%Generated by MTT

%Unity CR
OPERATOR Unity;
FOR ALL Causality,Input
LET Unity(Causality,Input) = Input;

%Linear CR: e = Gain*f (if DefaultCausality = flow) 
%           f = Gain*e (if DefaultCausality = effort)
OPERATOR Lin;
FOR ALL DefaultCausality, Gain, Causality, Input
SUCH THAT Causality = DefaultCausality
LET Lin(DefaultCausality, Gain, Causality, Input) = Gain*Input;

%Linear CR: e = (1/Gain)*f (if DefaultCausality = flow) 
%           f = (1/Gain)*e (if DefaultCausality = effort)
FOR ALL DefaultCausality, Gain, Causality, Input
SUCH THAT Causality NEQ DefaultCausality
LET Lin(DefaultCausality, Gain, Causality, Input) = (1/Gain)*Input;

%Default CR
% -- Unity CR
OPERATOR Default;
FOR ALL Causality,Input
LET Default(Causality,Input) = Unity(Causality,Input);

END;

EOF

#! /bin/sh

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

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

# Bourne shell script: matlab_tidy
# Tidies up Matlab format files
# P.J.Gawthrop 14 June 1990 13 Nov 1990, 19 Nov 93, April 1994, Dec 1994, 
# June 1995
# Copyright (c) P.J.Gawthrop, 1990, 1994, 1995.

#rm junk
cat $1 | \
tr "[A-Z]" "[a-z]" | \
sed -e "s/      //" \
    -e "s/      //" \
    -e "s/\$\$//" \
    -e "s/\.\.;;/../" \
    -e "s/\$/;/g" \
    -e "s/;;;//g" \
    -e "s/;;/;/g" \
    -e "s/\.;/./g" \
    -e "s/,;/,/g" \
    -e "s/\*\*/^/g" \
    -e "s/nocr /nocr/g" \
    -e "s/nocr\\nl/nocr/g" \
    -e "s/;;/;/g" \
    -e "s/:=/=/g" \
    -e "s/ e  /e/g" \
    -e "s/*;/* .../g" \
    -e "s/+;/+ .../g" \
    -e "s/     . //g" \
    -e "s/^;$//g" \
    -e "s/mttx0/X_0/g" \
    -e "s/mttu0/u_0/g" \
    -e "s/mtty0/y_0/g" \
    -e "s/mttdx/dX/g" \
    -e "s/mttx/X/g" \
    -e "s/mttu/u/g" \
    -e "s/mtty/y/g" \
    -e "s/mtta/A/g" \
    -e "s/mttb/B/g" \
    -e "s/mttc/C/g" \
    -e "s/mttd/D/g" \
    -e "s/mtte/E/g" \
    -e "s/mtttf/G/g" \
    -e "s/mttpar/MTTPAR/g"  \
    -e "s/\.\.\\$/../g"  \
    -e "s/\\$/;/g"  \
    -e "s/inv(mtte)\*/mtte\\\/g"\
    -e "s/\\$//" |\
tr -s '\012' '\012'\
>junk
mv junk $1

#! /bin/sh

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

# Bourne shell script: mtt_error
# Does the error handing for transformations and forces error exit.
# P.J.Gawthrop  August 1996
# Copyright (c) P.J.Gawthrop, 1996

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
# Revision 1.1  1996/08/25  09:20:00  peter
# Initial revision
#
###############################################################

error_file=$1;

# Test for errors and print if any
err_length=`wc -c <$error_file`
if [ $err_length != "0" ]
then
  echo MTT has failed with the following errors '...'
  cat $error_file
  exit 1
else
  exit 0
fi

#! /bin/sh

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

# Bourne shell script: ode2obs_r
# Odrinary differential equations to observer function equations
# P.J.Gawthrop 14 June 1991, 8 Aug 1991, 2 April 1992, 14 April 1994, 28 Dec 94,
#		12th July 1995, April 1996
# Copyright (c) P.J.Gawthrop 1991, 1992, 1994, 1995, 1996.

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.2  1996/08/25 08:38:14  peter
## Error handling added.
##
## Revision 1.1  1996/08/25 08:37:44  peter
## Initial revision
##
###############################################################

#Inform user
echo Creating $1_obs.r

# Remove the old log file
rm -f ode2obs_r.log

# Use reduce to accomplish the transformation
reduce >ode2obs_r.log << EOF


OFF Echo;
OFF Nat;
ON NERO;

in "$1_def.r";
MTTdxs := MTTdX;  %Save the symbolic form of dX

%Set default values - reset by sympar file.
MTTGPCNy := 2;
MTTGPCNu := 0;

in "$1_sympar.r";
in "$1_ode.r";


%Create the U vector of input derivatives.
MTTNuu := (MTTGPCNu+1)*MTTNu;
MTTNuu1 := MTTGPCNu*MTTNu;
Matrix MTTUU(MTTNuu,1);

FOR i := 1:MTTNu DO
  MTTUU(i,1) := MTTU(i,1);
END;

IF MTTGPCNu>0 THEN
BEGIN
  FOR i := 1:MTTNu DO
    MTTUU(i+MTTNu,1) := MTTdU(i,1);
  END;
END;

MTTU := MTTU;
%  MTTdU := MTTdU; ---- removed temporarily, needs def file change? 
MTTUU := MTTUU;

%Create the Y vector of output derivatives.
MTTNyy := (MTTGPCNy+1)*MTTNy;
Matrix MTTYY(MTTNyy,1);

FOR i := 1:MTTNy DO
  MTTYY(i,1) := MTTY(i,1);
END;

l := MTTNy;
FOR i := 1:MTTGPCNy DO
    FOR j := 1:MTTNy DO
    BEGIN
    l := l+1;
    MTTYY(l,1) := 0;
    FOR k := 1:MTTNx DO %Derivatives wrt x
    BEGIN
      xk := MTTX(k,1);
      MTTYY(l,1) := MTTYY(l,1) + df(MTTYY(l-MTTNy,1), xk, 1)*MTTdX(k,1);
    END;
    IF MTTNuu1>0 THEN
    FOR k := 1:MTTNuu1 DO %Derivatives wrt u
    BEGIN
      uk := MTTUU(k,1);
      MTTYY(l,1) := MTTYY(l,1) + df(MTTYY(l-MTTNy,1), uk, 1)*MTTUU(k+MTTNu,1);
    END;
  END;
END;

%Create O_x - derivative of YY wrt x
Matrix MTTO_x(MTTNyy,MTTNx);
FOR j := 1:MTTNx DO
  BEGIN
  xj := MTTX(j,1);
  FOR i := 1:MTTNyy DO
    BEGIN
    MTTO_x(i,j) := df(MTTYY(i,1), xj);
    END;
  END;

%Create O_u - derivative of YY wrt u (Assumes GPC Nu = 0)
MTTNNu := (MTTGPCNu+1)*MTTNu;
Matrix MTTO_u(MTTNyy,MTTNNu);
FOR j := 1:MTTNNu DO
  BEGIN
  uj := MTTu(j,1);	
  FOR i := 1:MTTNyy DO
    BEGIN
    MTTO_u(i,j) := df(MTTYY(i,1), uj);
    END;
  END;


%Create O_uu - derivative of O_u wrt u (Assumes GPC Nu = 0)
%This is a multi-dimensional matrix kth elements stacked sideways.
Matrix MTTO_uu(MTTNyy,MTTNNu*MTTNNu);
FOR k := 1:MTTNNu DO
  BEGIN
  uk := MTTu(k,1);	
  FOR j := 1:MTTNNu DO
    BEGIN
    FOR i := 1:MTTNyy DO
      BEGIN
      jk := j+(k-1)*MTTNu;
      MTTO_uu(i,jk) := df(MTTO_u(i,j), uk);
      END;
    END;
  END;

%Create O_ux - derivative of O_u wrt x 
%This is a multi-dimensional matrix kth elements stacked sideways.
Matrix MTTO_ux(MTTNyy,MTTNu*MTTNx);
FOR k := 1:MTTNx DO
  BEGIN
  xk := MTTx(k,1);	
  FOR j := 1:MTTNu DO
    BEGIN
    FOR i := 1:MTTNyy DO
      BEGIN
      jk := j+(k-1)*MTTNu;
      MTTO_ux(i,jk) := df(MTTO_u(i,j), xk);
      END;
    END;
  END;



%%Create the _obs.r file
OUT "$1_obs.r";

MTTYY := MTTYY;
MTTO_x := MTTO_x;
MTTO_u := MTTO_u;
MTTO_uu := MTTO_uu;
MTTO_ux := MTTO_ux;

write ";END;";

SHUT "$1_obs.r";
quit;

EOF

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

#! /bin/sh

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

# Bourne shell script: ode2odes_m

# Transforms descriptor matrix rep to step response

# Copyright (c) P.J.Gawthrop, 1996.

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.5  1996/08/18 12:01:26  peter
## Unified format of time responses.
##
## Revision 1.4  1996/08/16 13:04:46  peter
## Fixed problem with more than one output (y vector).
##
## Revision 1.3  1996/08/16 06:36:03  peter
## Removed u from default arg list.
##
## Revision 1.2  1996/08/15 16:24:43  peter
## Uses T in place of t to avoid name clash within function.
##
## Revision 1.1  1996/08/15 11:56:38  peter
## Initial revision
##
###############################################################

echo Creating $1_odes.m
rm -f ode2odes_m.log

if [ "$2" = "" ]; 
then
  PARAMS='T=[0:0.1:10]; x0=zeros(nx,1);'
  echo Using default parameter $PARAMS
else
  PARAMS=$2;
fi

PARAMS="$PARAMS ;"


$MATRIX << EOF > ode2odes_m.log 2>mtt_error

  %Read in parameters
  $1_numpar;

  [nx,ny,nu,nz,nyz] = $1_def;
  t=0;	%Just in case it appears in the parameter list.

  $PARAMS

  %Defaults
  if exist('T')==0
    T=[0:0.1:10]
  end;

  if exist('x0')==0
    x0 = zeros(nx,1);
  end;

  [n,m]=size(T);
  if m>n
    T=T';
  end;

if nx>0
  x = lsode('$1_ode', x0, T);
  write_matrix([T,x], '$1_odes');
else
  x = zeros(size(T));
end;

if ny>0
  i=0;
  for tt=T'
    i=i+1;
    y(i,:) = $1_odeo(x(i,:),tt)';
  end;
  write_matrix([T,y], '$1_odeso');
end;

EOF

err_length=$(wc -c <mtt_error)

# Test for errors and print if any
if [ $err_length != "0" ]
then
  echo MTT has failed with the following errors '...'
  cat mtt_error
  exit 1
else
  exit 0
fi

#! /bin/sh

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

# Bourne shell script: ode_r2m
# Reduce ODE to simulab  ODE  
# P.J.Gawthrop 14 June 1991, 12 Jan 1994, April 1994, Jan 95.
# Copyright (c) P.J.Gawthrop 1991, 1994, 1995, 1996

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.1  1996/08/18 12:03:49  peter
## Initial revision
##
###############################################################

#Inform user
echo Creating $1_ode.m
echo Creating $1_odeo.m


# Remove the old log file
rm -f ode_r2m.log

# Use reduce to accomplish the transformation
reduce >ode_r2m.log << EOF

%Read the reduce definitions file
in "$1_def.r";

%Read the reduce  ODE   file
in "$1_ode.r";

%Set up the number of argument variables to zero in case the user has forgotten
MTTNVar := 0;

%Read the parameter file
in "$1_sympar.r";


ON NERO;        % Suppress zero elements

%Define the common part of the functions.

PROCEDURE common;
BEGIN
  IF MTTNvar>0 THEN
  BEGIN
    write "% Read in the parameters";
    write "global ...;;";
    FOR i := 1:MTTNvar DO
    BEGIN
       IF i<MTTNvar THEN write MTTVar(i,1), " ..."
                    ELSE write MTTVar(i,1), ";"
    END;
  END;
  
  write "% Read in the input";
  write "u = $1_input(t)";

  write "% Read in the arguments";
  write "$1_args";
  
  write "% Set up the State variables";
  FOR i := 1:MTTNx DO
  BEGIN
    write "MTTx", i, " = x(", i, ");";
  END;
  
  write "% Set up the Input variables";
  IF MTTNu>0 THEN
  FOR i := 1:MTTNu DO
  BEGIN
    write "MTTu", i, " = u(", i, ");";
  END;
END;
  
% Firstly do the dx = f(x,t) function.
OUT "$1_ode.m";
  
write "function MTTdX = $1_ode(x,t);";
write "% dX = $1_ode(x,t);";
write "%ODE in Simulab form for system $1;;";
write "%File $1_ode.m;;";
write "%Generated by MTT;;";

common();

%Fortran switches - one line expressions
OFF echo;
ON fort$
cardno!* := 1$
fortwidth!* := 100$
OFF period$

MTTdx := MTTdx;

SHUT "$1_ode.m";
OFF fort;
% Now do the y = g(x,t) function.
OUT "$1_odeo.m";
  
write "function MTTy = $1_odeo(x,t);";
write "% dX = $1_odeo(x,t);";
write "%ODE in Simulab form for system $1;;";
write "%File $1_odeo.m;;";
write "%Generated by MTT;;";

common();

%Fortran switches - one line expressions
OFF echo;
ON fort$
cardno!* := 1$
fortwidth!* := 100$
OFF period$

MTTy := MTTy;

SHUT "$1_odeo.m";

#! /bin/sh

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

# Bourne shell script: rbg2abg_m
#
###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.4  1996/08/24 17:57:41  peter
## Removed `touch mtt_info.txt'
##
## Revision 1.3  1996/08/24 14:58:56  peter
## Fixed missprint - mtt_err --> mtt_error
##
## Revision 1.2  1996/08/24 14:33:07  peter
## Error handling included.
##
## Revision 1.1  1996/08/05 12:22:16  peter
## Initial revision
##
###############################################################

# Raw bond graph to structured acausal bond graph: mfile format
# The structured BG is described by four matrices:
#   junctions: describes the junctions -
#     Column 1   Index of component (see components)
#     Column 2.. Indices of bonds on junction (see jbonds)
#                  rows may be padded with zeros
#
#   jbonds:    describes the true bonds -
#     Column 1     Index of component (see components)
#     Column 2     Direction
#                    1 away from junction
#                   -1 towards junction
#     Column 3     Causality 
#                     1 stroke at arrow end
#                    -1 stroke at other end
#                     0 no stroke (ie acausal)
#     Columns 4:5  x,y coordinates of non-arrow end (in Fig)
#     Columns 6:7  x,y coordinates of arrow end (in Fig)
#
#   mbonds:    describes the modulation bonds -
#     Column 1     Index of component at non-arrow end
#     Column 2     Index of component at arrow end
#     Columns 4:5  x,y coordinates of non-arrow end (in Fig)
#     Columns 6:7  x,y coordinates of arrow end (in Fig)
#
#   components: describes the comonents -
#     Column 1     Type of component
#                    0  0-junction
#                    1  1-junction
#                    2  R component
#                    3  C component
#                    4  I component
#                    11 TF component
#                    12 GY component
#                    13 A  component
#     Column 2     Flag for explicit junction
#                    4  Explicit component (4 means text type)
#                    0  Implicit component
#     Column 3-14    fig file information (see fig documentation)
#
# P.J.Gawthrop May 1996
# Copyright (c) P.J.Gawthrop, 1996.


# Remove the old log file
rm -f rbg2abg_m.log
rm -f $1_abg.m

#Inform user
echo Creating $1_abg.m

# Use matrix manipulation to accomplish the transformation
$MATRIX > rbg2abg_m.log  2>mtt_error.txt << EOF

  infofile = 'mtt_info.txt';
  %Convert from the fig version of the bonds to a structured version
  [rbonds,rstrokes,rcomponents,rports,n_ports] = $1_rbg;
  [bonds,components] = rbg2abg(rbonds,rstrokes,rcomponents,rports,infofile);

  %Write the function m-file for the causal bond graph
  filename = '$1_abg.m';
  c = '%';
  fprintf(filename, 'function [bonds,components,n_ports] = $1_abg\n');
  fprintf(filename, '%s [bonds,components] = $1_abg\n', c);
  fprintf(filename, '%s Acausal bond graph created by MTT on %s\n\n', ...
                     c, date);
  mat2mfile(bonds, 'bonds', filename);
  mat2mfile(components, 'components', filename);
  fprintf(filename, 'n_ports = %1.0f;\n', n_ports);

EOF

cat mtt_info.txt

mtt_error mtt_error.txt

#! /bin/sh

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

# Bourne shell script: rbg_fig2m
# Raw  bond-graph conversion from fig to matlab
# P.J.Gawthrop May 1996
# Modified June 1996 to handle xfig version 3
# Copyright (c) P.J.Gawthrop, 1996.

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


# Inform user
echo "Creating $1_rbg.m"
echo "Creating $1_cmp.m"
echo "Creating $1_fig.fig"

rm -f $1_fig.fig

# The following horrible sed thing is to replace tab by ----
# this is because I can't get sed to recognise tabs even if FS=" ".
sed 's/	/---- /'<$1_abg.fig > $1_abg.mtt

# This is the main transformation using gawk
cat $1_lbl.txt $1_abg.mtt |  gawk -f $MTTPATH/rbg_fig2m.awk $1

rm -f $1_abg.mtt

#! /bin/sh

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

# Bourne shell script: sm2sr_m

# Transforms descriptor matrix rep to step response


# Copyright (c) P.J.Gawthrop, 1996.

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.10  1996/08/19 15:34:14  peter
## Now based on state - not descriptor matrices.
##
## Revision 1.9  1996/08/18 12:01:45  peter
## Unified format of time responses.
##
## Revision 1.8  1996/08/15 16:23:39  peter
## Uses T in place of t to avoid name clash within function.
##
## Revision 1.7  1996/08/15 11:54:32  peter
## Now has optional initial condition.
##
## Revision 1.6  1996/08/14 09:15:02  peter
## Now encapsulates the data into an m file
## Name changed
##
## Revision 1.5  1996/08/12 20:19:06  peter
## Fiddled about with parameter passing - still not right
##
## Revision 1.4  1996/08/11 19:48:04  peter
## Parameter passing added.
##
## Revision 1.3  1996/08/11 10:38:30  peter
## Minor debugging changes.
##
## Revision 1.2  1996/08/11 10:06:50  peter
## Now allows e format and NaN numbers.
##
## Revision 1.1  1996/08/11 09:45:46  peter
## Initial revision
##
###############################################################

echo Creating $1_sr.m
echo Creating $1_sro.m
rm -f sm2sr_m.log
rm -f $1_sr.m
rm -f $1_sro.m
rm -f mtt_error

if [ "$2" = "" ]; 
then
  PARAMS='T=[0:0.1:10];u0=ones(nu,1);x0=zeros(nx,1);'
  echo Using default parameter $PARAMS
else
  PARAMS=$2;
fi

PARAMS="$PARAMS ;"

$MATRIX << EOF > sm2sr_m.log 2>mtt_error

  % Read in the numeric system parameters
  $1_numpar;

  [nx,ny,nu,nz,nyz] = $1_def;
  t=0;	%Just in case it appears in the parameter list.

  $PARAMS

  %Defaults
  if exist('T')==0
    T=[0:0.1:10]
  end;

  if exist('u0')==0
    u0 = ones(nu,1);
  end;

  if exist('x0')==0
    x0 = zeros(nx,1);
  end;


  [n,m]=size(T);
  if m>n
    T=T';
  end;

  [A,B,C,D] = $1_sm;
  [y,x] = sm2sr(A,B,C,D,T,u0,x0);

  if nx>0
    write_matrix([T x], '$1_sr');
  end;
  if ny>0
    write_matrix([T y], '$1_sro');
  end;

EOF

err_length=$(wc -c <mtt_error)

# Test for errors and print if any
if [ $err_length != "0" ]
then
  echo MTT has failed with the following errors '...'
  cat mtt_error
  exit 1
else
  exit 0
fi

#! /bin/sh

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

# Bourne shell script: sm_r2m
# Reduce constained state-space A,B,C,D and E matrices to Matlab
# P.J.Gawthrop July 1990, Feb 1991, May 1991, Sep 91, Dec 93, Jan 94, May 94.
# Copyright (c) P.J.Gawthrop, 1989, 1990, 1991, 1993, 1994, 1996.

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.3  1996/08/19 18:38:43  peter
## Read in the numerical params etc ...
##
## Revision 1.2  1996/08/19 16:25:18  peter
## Fixed bug when Nx=0.
##
## Revision 1.1  1996/08/19 15:14:51  peter
## Initial revision
##
###############################################################


# Inform user
echo Creating $1_sm.m

# Remove the old log file
rm -f sm_r2m.log

# Use reduce to accomplish the transformation
reduce >sm_r2m.log << EOF

%Read in the definitions file
IN "$1_def.r";

%Read the reduce state-space A,B,C and D matrices file
IN "$1_sm.r";

%Set up the number of argument variables to zero in case the user has forgotten
MTTNVar := 0;

%Read in the parameter file
IN "$1_sympar.r";


OUT "$1_sm.m";
%Headings - Matlab style
%(Note. The ;; are deleted by for2mat)

write "function [MTTA,MTTB,MTTC,MTTD] = $1_sm;;";
write "%function [MTTA,MTTB,MTTC,MTTD] = $1_sm;;";

write "%Linearised state matrices for system $1";
write "%File $1_sm.m";
write "%Generated by MTT";
IF MTTNvar>0 THEN
BEGIN
  write "% Set the parameters";
  write "global ...;;";
  FOR i := 1:MTTNvar DO
  BEGIN
     IF i<MTTNvar THEN write MTTVar(i,1), "..."
                  ELSE write MTTVar(i,1), ";"
  END;
END;


%Fortran switches - one line expressions
OFF echo;
ON fort$
cardno!* := 1$
fortwidth!* := 100$
OFF period$

write "MTTA = zeros(", MTTNx, ",", MTTNx, ");";
write "MTTB = zeros(", MTTNx, ",", MTTNu, ");";
write "MTTC = zeros(", MTTNy, ",", MTTNx, ");";
write "MTTD = zeros(", MTTNy, ",", MTTNu, ");";

ON NERO;	% Suppress zero elements.
MTTA := MTTA;
MTTB := MTTB;
MTTC := MTTC;
MTTD := MTTD;

OFF FORT;
write "";

SHUT "$1_sm.m";

#! /bin/sh

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

# Bourne shell script: sspar2ss_r
# Reduce steady steady-state parameter file to reduce steady-state:
#	X0, and Y0 and U0;
# P.J.Gawthrop 10 July 1990, 22 April 1991, Dec 1993, April 1994
# Copyright (c) P.J.Gawthrop, 1990, 1994.

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

#Inform user
echo Creating $1_ss.r

# Remove the old log file
rm -f sspar2ss_r.log

# Use reduce to accomplish the transformation
reduce >sspar2ss_r.log << EOF

%Read the definition file
in "$1_def.r";

%Read the constitutive relationships (May be used in rss file). 
in "$1_cr.r";

%Read the steady-state values 
in "$1_sspar.r";

%Read in the parameter file
in "$1_sympar.r";


OFF Echo;
OFF Nat;
OUT "$1_ss.r";


%Write out the steady-state values (if any)
write "%Steady-state values";
IF MTTNx>0 THEN
BEGIN
  write "matrix MTTX0(", MTTNx, ",1)";
END;
MTTX0 := MTTX;

IF MTTNu>0 THEN
BEGIN
  write "matrix MTTU0(", MTTNu, ",1)";
END;
MTTU0 := MTTU;

IF MTTNy>0 THEN
BEGIN
  write "matrix MTTY0(", MTTNy, ",1)";
END;
MTTY0 := MTTY;

%Write them out by component - but only if values have been set

%(This assumes that if element 1 not set then the rest aren't either)
clear MTTx1;
IF MTTNx>0 THEN
IF MTTX0(1,1) NEQ MTTX1 THEN
FOR i := 1:MTTNx DO
BEGIN
  write "MTTX", i, " := ", MTTX0(i,1);
END;

clear MTTu1;
IF MTTNu>0 THEN
IF MTTu0(1,1) NEQ MTTu1 THEN
FOR i := 1:MTTNu DO
BEGIN
  write "MTTu", i, " := ", MTTu0(i,1);
END;

%Derivatives are zero in the steady-state
IF MTTNz>0 THEN
FOR i := 1:MTTNz DO
BEGIN
  write "MTTdz", i, " := 0";
END;

write "END;";

SHUT "$1_ss.r";
quit;
EOF

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

#! /bin/sh

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

# Bourne shell script: mknumpar
# Creates the deafault numerical parameters file for the system

# Copyright (c) P.J.Gawthrop 1996

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.6  1996/08/20 08:21:24  peter
## *** empty log message ***
##
## Revision 1.5  1996/08/16 14:27:25  peter
## Make a dummy function when no parameters (MTTNVar=0).
##
## Revision 1.4  1996/08/15 12:46:59  peter
## Removed a spuriuos echo.
##
## Revision 1.3  1996/08/12 20:20:59  peter
## Explicit parameters via out put list.
## Parameter conversion stuff deleted.
##
## Revision 1.2  1996/08/11 10:41:52  peter
## Fixed various bugs.
##
## Revision 1.1  1996/08/11 08:48:07  peter
## Initial revision
##
###############################################################

# Inform user
echo Creating $1_numpar.m

# Remove the old log file
rm -f mknumpar.log


# Use reduce to accomplish the transformation
reduce >mknumpar.log << EOF

%Set up the number of argument variables to zero in case the user has forgotten
MTTNVar := 0;

%Read the symbolic parameters file
IN "$1_sympar.r";

OUT "$1_numpar.m";

write "% Script $1_numpar;";
write "% Created by MTT";

IF MTTNvar>0 THEN
BEGIN
  write "global ...;;";
  FOR i := 1:MTTNvar DO
  BEGIN
     IF i<MTTNvar THEN write MTTVar(i,1), " ..."
                  ELSE write MTTVar(i,1), ";";
  END
END; 


IF MTTNvar>0 THEN
BEGIN
  FOR i := 1:MTTNvar DO
  BEGIN
    write MTTVar(i,1), "  = 1;	% Default value";
  END;
END;

if MTTNVar=0 THEN
  write "foo=0";

SHUT "$1_numpar.m";

#! /bin/sh

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

# Bourne shell script: tex2doc
# Encapsulates a tex file in a document.

# P.J.Gawthrop August 1996
# Copyright (c) P.J.Gawthrop, 1996.

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

# Inform user
echo "Creating $1_$2.doc"

# This is the main transformation using gawk
gawk -f $MTTPATH/tex2doc.awk $1_$2.txt  > $1_$2.doc

function  structure = AE_eqn(bond_number,bonds,direction,cr,args, ...
                            structure,eqnfile);
% AE_eqn. Equation generation for effort amplifier.

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<6
  eqnfile = 'stdout';
end;

% There are 2 ports; extract the information
e_1 = bonds(1,1);
f_1 = bonds(1,2);
e_2 = bonds(2,1);
f_2 = bonds(2,2);


% Flow on port 1 is always zero
   fprintf(eqnfile, '%s := 0;\n', ...
      varname(bond_number(1), -1));
  
  
LHS_cause = 1;
RHS_cause = 1;

if e_2 == -1 				% Unicausal: e_2 := e_1
  LHS_number = bond_number(2);
  RHS_number = bond_number(1);
else 	                                % Bicausal: e_1 := e_2      
  LHS_number = bond_number(1);
  RHS_number = bond_number(2);
end

oneeqn(LHS_number,LHS_cause,RHS_number,RHS_cause,cr,args,eqnfile);

 

function structure =  C_eqn(bond_number,bonds,direction,cr,args, ...
                            structure,eqnfile);
% C_eqn: equation generation for C component


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


			
if nargin<7
  eqnfile = 'stdout';
end;

% Unicausal version
CorI = 1;
structure = cieqn(bond_number,bonds,direction,cr,args, structure, ...
                 CorI, eqnfile);

function [bonds,status] = FMR_cause(bonds);

% Causality for Flow-modulated R component

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


status = -1;

% There are 2 ports; extract the information
e_1 = bonds(1,1);
f_1 = bonds(1,2);
e_2 = bonds(2,1);
f_2 = bonds(2,2);

% e_2 must be effort causality (imposes zero effort)
 if e_2 == 1 				% Conflict
   status = 1;
 else 					% Do the rest of the causality
   if e_2==0 				% Set to the fixed causality
     e_2 = -1;
   end;
   
   number_set = sum(sum([e_1 f_1 f_2]~=zeros(1,3) ));
   if number_set<2 			% Under causal
     status = -1;
   elseif number_set==2 		% Set the causality
     if f_2 == -1 			% Unicausal
       if e_1 == 0
 	e_1 = f_1;
       else
	 f_1 = e_1;
       end;
     elseif f_2 == 1 			% Bicausal
       e_1 = 1;
       f_1 = -1;
     elseif f_2 == 0
       if e_1==f_1
	 f_2 = -1;                      % Unicausal
       else
	 f_2 = 1;                      % Bicausal
       end;
     end;
     status = 0;
   elseif number_set==3 		% Check the causality
     if ( (f_2==-1)&(e_1~=f_1) )|( (f_2==1)&(e_1==f_1) )
       status = 1;
     else
       status = 0;
     end;
   end;
 end;
 
 bonds = [e_1 f_1
           e_2 f_2];

function structure =  SS_eqn(bond_number,bonds,direction,cr,args, ...
                            structure,filenum);

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %% Revision 1.3  1996/08/09 14:08:04  peter
% %% Empty effort and flow attributes replaced by 'external'.
% %%
% %% Revision 1.2  1996/08/08 16:38:19  peter
% %% Changed to filenumer type of fprintf
% %%
% %% Revision 1.1  1996/08/08 16:19:08  peter
% %% Initial revision
% %%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% The string 'cr' contains the effort attributes;  
% 'arg' contains the flow attributes. 
% external indicates an external variable
% internal indicates an internal variable
% anything else is the name of a constant.
% Blank is taken to be external

effort_attribute = cr;
flow_attribute = args;

if length(effort_attribute)==0
  effort_attribute = 'external';
end;

if length(flow_attribute)==0
  flow_attribute = 'external';
end;

inputs = structure(3);
outputs = structure(4);
zero_outputs = structure(5);

% Effort
if strcmp(effort_attribute, 'external')
  if bonds(1,1)==-1 % Source
    inputs = inputs+1;
    fprintf(filenum, '%s := MTTu(%1.0f,1);\n', ...
        varname(bond_number,1),inputs);
  else % Sensor
    outputs = outputs+1;
    fprintf(filenum, 'MTTy(%1.0f,1) := %s;\n', ...
        outputs, varname(bond_number,1));
  end;
elseif strcmp(effort_attribute, 'internal')
  % Do nothing
else % named constant
  if bonds(1,1)==-1 % Source
    fprintf(filenum, '%s := %s;\n', ...
	varname(bond_number,1), effort_attribute);
  else % Sensor
    if strcmp(effort_attribute, 'zero') %Zero output
      zero_outputs = zero_outputs + 1;
      fprintf(filenum, 'MTTyz%1.0f := %s;\n', ...
	  zero_outputs, varname(bond_number,1));
      fprintf(filenum, '%s := MTTUi%1.0f;\n', ...
	  varname(bond_number,-1), zero_outputs);
    else
      mtt_info([effort_attribute, ' not appropriate for an output ']);
    end;
  end;
end;

% Flow
if strcmp(flow_attribute, 'external')
  if bonds(1,2)==1 % Source
    inputs = inputs+1;
    fprintf(filenum, '%s := MTTu(%1.0f,1);\n', varname(bond_number,-1),inputs);
  else % Sensor
    outputs = outputs+1;
    fprintf(filenum, 'MTTy(%1.0f,1) := %s;\n', outputs, ...
        varname(bond_number,-1));
  end;
elseif strcmp(flow_attribute, 'internal')
  % Do nothing
else % Named constant
  if bonds(1,2)==1 % Source
    fprintf(filenum, '%s := %s;\n', ...
	varname(bond_number,-1), flow_attribute);
  else % Sensor
    if strcmp(flow_attribute, 'zero') %Zero output
      zero_outputs = zero_outputs + 1;
      fprintf(filenum, 'MTTyz%1.0f := %s;\n', ...
	  zero_outputs, varname(bond_number,-1));
      fprintf(filenum, '%s := MTTUi%1.0f;\n', ...
	  varname(bond_number,1), zero_outputs);
    else
      mtt_info([effort_attribute, ' not appropriate for an output ']);
    end;
  end;
end;

  
structure(3) = inputs;
structure(4) = outputs;
structure(5) = zero_outputs;  

  
  

function  structure = TF_eqn(bond_number,bonds,direction,cr,args, ...
                            structure,eqnfile);

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

			
if nargin<6
  eqnfile = 'stdout';
end;

% Needs port labels etc...

% Effort
if bonds(1,1)==1
  oneeqn(bond_number(2),1,bond_number(1),1,cr,args,eqnfile);
else
  oneeqn(bond_number(1),1,bond_number(2),1,cr,args,eqnfile);
end;

% Flow
if bonds(1,2)==-1
  oneeqn(bond_number(2),-1,bond_number(1),-1,cr,args,eqnfile);
else
  oneeqn(bond_number(1),-1,bond_number(2),-1,cr,args,eqnfile);
end;

function [bonds,status] = juncause(bonds,jun,cause)
% [bonds,status] = juncause(bonds,jun,cause)

% Causality for  either effort or flow on  either zero or one junctions

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


j =  (3-cause)/2; % j is 1 for effort, 2 for flow

[n_bonds,junk] = size(bonds);
[causing_bond, n,other_bonds,m] = getindex(bonds(:,j),jun);
if n>1 % over causal
  status = 1;
elseif n==1 %causal
  status = 0;
  bonds(other_bonds(:,1),j) = -jun*ones(m,1);
else  % undercausal - try other way
  [causing_bond, n,other_bonds,m] = getindex(bonds(:,j),-jun);
  if n==n_bonds % over causal
    status = 1;
  elseif n==n_bonds-1 %causal
    status = 0;
    bonds(other_bonds(:,1),j) = jun*ones(m,1);
  else  % undercausal
    status = -1;
  end;
end;

function juneqn(bond_number,bonds,direction,cr,args,jun,filenumber)

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %% Revision 1.1  1996/08/08 16:38:50  peter
% %% Initial revision
% %%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Equations for  both  effort and flow on  either zero or one junctions

% Same causality as junction
column =  (3-jun)/2; 
[causing_bond,n,other_bonds,m] = getindex(bonds(:,column),jun);
for i=other_bonds'
  oneeqn(bond_number(i), jun, bond_number(causing_bond), jun, '', '',filenumber);
end;


% Opposite causality to junction
column =  (3+jun)/2;
[caused_bond,n,other_bonds,m] = getindex(bonds(:,column),jun);
fprintf(filenumber, '%s\t:= \n',  varname(bond_number(caused_bond),-jun));
for i=other_bonds'
  term_sign = -direction(caused_bond,column)*direction(i,column);
  fprintf(filenumber, '\t\t%s %s\n', sign2name(term_sign), ...
      varname(bond_number(i),-jun));
end;
fprintf(filenumber, ';\n');