#!/bin/sh

#MAKE CHANGES - lines beginning so indicate where parameter values 
#	can be changed

root="scr"		# first field of name of output files
E="1.0e-8"		# starting enzyme concentration
S="1.0e-3"		# starting substrate concentration
reactionplot=`expr 0`	# plot concentration of all reactant and product
			# species, over full (log) time scale, for both
			# mixing and jump expts, as appropriate
productplot=`expr 0`	# plot concentration of final product, from 
			# initial to 95%=factor of final concentration
factor="0.95"		# see above
output=`expr 1`		# no print output of any sort, only create the files
			# appropriate for mixing and jump simulations
trackident="P1"		# identifier of species to track, default = product
tracknum=0		# ordinal in output line of species to track 
Sident="S"		# ditto for substrate
Snum=0
Eident="E"		# ditto for enzyme
Enum=0

usage="\
$0  [-E enzyme_conc] [-S substrate_conc] [-c enzyme_conc==substrate_conc] \
	[-[no]react]  [-[no]prod] [-[no]output] \
	[output_file_basename]"

while [ $# -gt 0 ]
do
	case $1 in
	-[Ee])	shift;E="$1";shift;;
	-[Ss])	shift;S="$1";shift;;
	-c*)	shift;E="$1";S="$1";shift;;
	-nor*)	reactionplot=`expr 0`;shift;;
	-nop*)	productplot=`expr 0`;shift;;
	-noo*)	output=`expr 0`;shift;;
	-r*)	reactionplot=`expr 1`;output=`expr 1`;shift;;
	-p*)	productplot=`expr 1`;output=`expr 1`;shift;;
	-o*)	output=`expr 1`;shift;;
	-*)	echo "ERROR in command line argument: $1"
		echo $usage
		exit;;
	*)	root="$1";shift;;
	esac
done



# MIXING SIMULATION

# attach MIX identifier to $root
rootmix=${root}.mix

# create rootmix.m file
cat <<EOF | grep -v "^#" >${rootmix}.m
# rate constant		reactants -> products
#MAKE CHANGES - in the rate constants and mechanism
#values for lysozyme, from banerjee et al, jbc 250, 4355 (1975)
#pH = 6.3, 25 deg
20	EU2 -> EU1
200	EU1 -> EU2
6.25e3	EU1 -> E + S
1.e8	E + S -> EU1
1.e8	E + S -> ESB1
5e3	ESB1 -> E + S
14	ESB1 -> ESB2
5	ESB2 -> ESB1
.46	ESB2 -> ESG
.13	ESG -> ESB2
.038	ESG -> E + P1 + P2
EOF

# create rootmix.p file
cat <<EOF | grep -v "^#" >${rootmix}.p
#iopt
1
#itol	istate	jt	itask
4	1	2	1
#ml	mu	ixpr	mxstep	mxhnil	mxordn	mxords
#iwork1	iwork2	iwork5	iwork6	iwork7	iwork8	iwork9
0	0	0	500	0	0	0
#tcrit	h0	hmax	hmin
#rwork1	rwork5	rwork6	rwork7
0.	0.	0.	0.
#start_time  end_time  {time_step|0}  deflt_tolerance  [first_time  time_factor]
#MAKE CHANGES - in the time range and time step or time factor
#	be careful; note that the format is different for time step vs factor
#model for time step incrementation
#0.0	0.0001	.000003	1.0e-10
#model for time factor incrementation
#0	1.0e6	0	1.0e-10	1.e-9	1.2589254117
0	1.0e8	0	1.0e-10	1.e-9	1.2589254117
#species	initial concn		absol. tolerance	relat. tolerance
#MAKE CHANGES - in the initial concentrations
E	$E	1.e-15	1.e-10
S	$S	1.e-15	1.e-10
EU2	0	1.e-15	1.e-10
EU1	0	1.e-15	1.e-10
ESB1	0	1.e-15	1.e-10
ESB2	0	1.e-15	1.e-10
ESG	0	1.e-15	1.e-10
P1	0	1.e-15	1.e-10
P2	0	1.e-15	1.e-10
EOF

# run it 
react ${rootmix} >${rootmix}.r

# and set some variable values
# set tracknum = ordinal in output line of species tracked (e.g., product)
# set Snum and Enum similarly
varnames1="`grep time ${rootmix}.r`"
((count = 0))
for i in $varnames1
do
	((count = count + 1))
	echo $i | fgrep -xs $trackident && tracknum=$count
	echo $i | fgrep -xs $Sident && Snum=$count
	echo $i | fgrep -xs $Eident && Enum=$count
done
lastline1="`tail -1 ${rootmix}.r`"
firstline1="`cat ${rootmix}.r|awk '/time/ {getline;print;exit}'`"
set $firstline1
eval first_conc1="\$$tracknum"
set $lastline1
eval last_conc1="\$$tracknum"

# label() - function to output a label, derived from output data file
# and the following arguments:
# $1 = file_name, for results of reaction simulation
# $2 = type of simulation (MIX or JUMP)
# $3 = starting concentration of species monitored for halftime calculation
# $4 = final concentration of above
# $5 = tracknum = ordinal of species monitored
# $6 = Snum = ditto for substrate
# $7 = Enum = ditto for substrate

label() {
egrep -v 'e\-3[8-9]|e\-[4-9][0-9]' ${1} |
awk -v awkvar0=$2 -v awkvar1=$3 -v awkvar2=$4 -v awkvar3=$5 -v awkvar4=$6 -v awkvar5=$7 '
/time/ {
	getline
	t = $1
	c = $awkvar3
	while(getline > 0) {
	if (flag == 0 && ($awkvar3 - awkvar1)/(awkvar2 - awkvar1) > 0.1) {
		flag = 1
		rate = $awkvar3 / $1
	}
	if (flag == 1 && ($awkvar3 - awkvar1)/(awkvar2 - awkvar1) > 0.5) {
		flag = 2
		halftime = t + ($1 - t)*(0.5*(awkvar2 - awkvar1) + awkvar1 - c) / ($awkvar3 - c)
	}
	line = $0
	t = $1
	c = $awkvar3
	}
}
END	{
	nf = split(line,fields)
	printf "%s, final concentrations:\n\tE=\t%20.3e\n\tS=\t%20.3e\n\tP=\t%20.3e\n\t(P+S)=\t%20.3e\n", 
		awkvar0,
		fields[awkvar5], fields[awkvar4], fields[awkvar3], 
		fields[awkvar4]+fields[awkvar3]
	printf "product halftimes:\n\tt(1/2)=\t%20.3e\n\ttau=\t%20.3e\n\t1/tau=\t%20.3e\n",
		halftime, halftime/0.691, 0.691/halftime
	printf "product reaction rates:\n\tv_o(0-10%%)=\t%12.3e\n\tv_o(0-10%%)/E_o=\t%12.3e\n",
		rate, rate/fields[awkvar5]
}' -
}


# PLOT DATA

if [ $output -eq 0 ]
then
	exit
fi

label ${rootmix}.r "MIX" $first_conc1 $last_conc1 $tracknum $Snum $Enum
echo
cat ${rootmix}.r|awk '/time/ {print;getline;print}; {line=$0}; END {print line}'
echo

if [ $reactionplot -eq 1 ]
then
	cat ${rootmix}.r | prep | ttygraph
	echo "\n"
fi

if [ $productplot -eq 1 ]
then
	cat ${rootmix}.r |
	prep -c $tracknum $first_conc1 $last_conc1 $factor |
	ttygraph
echo "\n"
fi