#!/bin/sh

root="scr"		# first field of name of output files
E="1.0e-4"		# starting enzyme concentration
S="1.0e-3"		# starting substrate concentration
nojump=`expr 0`		# if 1, no simulation of jump expt after
			# simulate mixing expt
jumponly=`expr 0`	# if 1, simulate mixing and then jump expt, 
			# but no output for mixing simulation
reactionplot=`expr 1`	# plot concentration of all reactan and product
			# species, over full (log) time scale, for both
			# mixing and jump expts, as appropriate
productplot=`expr 1`	# 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

# rate constants for mixing simulation
mixkeu2_eu1=1.0e3	#EU2 -> EU1
mixkeu1_eu2=1.5e4	#EU1 -> EU2
mixkeu1_e_s=0.5e5	#EU1 -> E + S
mixke_s_eu1=4.0e7	#E + S -> EU1

# rate constants for jump simulation
jumpkeu2_eu1=2.0e3	#EU2 -> EU1 -- 2 x faster
jumpkeu1_eu2=1.5e4	#EU1 -> EU2
jumpkeu1_e_s=1.0e5	#EU1 -> E + S -- 2 x faster
jumpke_s_eu1=4.0e7	#E + S -> EU1

usage="\
$0  [-E enzyme_conc] [-S substrate_conc] [-c enzyme_conc substrate_conc] \
	[-prodonly]  [-noprod] [-jumponly] [-nojump]  [-nooutput] \
	[output_file_basename]"

while [ $# -gt 0 ]
do
	case $1 in
	-[Ee])	shift;E="$1";shift;;
	-[Ss])	shift;S="$1";shift;;
	-c*)	shift;S="$1";E="$1";shift;;
	-noreact)	reactionplot=`expr 0`;shift;;
	-noprod)	productplot=`expr 0`;shift;;
	-jumponly)	jumponly=`expr 1`;shift;;
	-nojump)	nojump=`expr 1`;shift;;
	-nooutput)	output=`expr 0`;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
# values for lysozyme, from imoto et al, the enzymes, 1972, and fersht text
# pH = 6.3, 25 deg
$mixkeu2_eu1	EU2 -> EU1
$mixkeu1_eu2	EU1 -> EU2
$mixkeu1_e_s	EU1 -> E + S
$mixke_s_eu1	E + S -> EU1
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]
#0.0	0.0001	.000003	1.0e-10
#0	1.0e6	0	1.0e-10	1.e-9	1.2589253
0	1.0e6	0	1.0e-10	1.e-9	1.2589254117
#0.0	1.0	1.0	1.0e-10
#species	initial concn		absol. tolerance	relat. tolerance
EU2	0	1.e-15	1.e-10
EU1	0	1.e-15	1.e-10
E	$E	1.e-15	1.e-10
S	$S	1.e-15	1.e-10
EOF

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

# and set some variable values
varnames1="`grep time ${rootmix}.r`"
lastline1="`tail -1 ${rootmix}.r`"
firstline1="`cat ${rootmix}.r|awk '/time/ {getline;print;exit}'`"
set $firstline1
first_conc1="$2"
set $lastline1
last_conc1="$2"


# JUMP SIMULATION

# attach JUMP identifier to $root
rootjump=${root}.jump

# create new rootjump.m file,
# with rate constants changed according to jump, e.g., 10 degc delT
cat <<EOF | grep -v "^#" >${rootjump}.m
#rate constant		reactants -> products
#after change from T1 -> T2 (delT = 10 degC, delH = 10 kcal, Q=2)
$jumpkeu2_eu1	EU2 -> EU1
$jumpkeu1_eu2	EU1 -> EU2
$jumpkeu1_e_s	EU1 -> E + S
$jumpke_s_eu1	E + S -> EU1
EOF

# create new rootjump.p file
# from rootmix.p file
# with all concentrations set to those at end of mixing simulation
cat ${rootmix}.p |
awk -v awkvar1="$varnames1" -v awkvar2="$lastline1" '
BEGIN		{
	nspecies=split(awkvar1, species)
	nconcn=split(awkvar2, concn)
}
#10 data points per decade time
NR == 5		{
	print "0	1.0e6	0	1.0e-10	1.e-9	1.2589254117"
	next
}
NR >= 6		{
	for (i=2; i<=nspecies; i++) {
		if ($1 ~ species[i]) {
			$2 = concn[i]
			print $0
			break
		}
	}
	if (i > nspecies) print ERROR
	next
}
		{
	print
}' - >${rootjump}.p


# run it 
#../../reaction/react ${rootjump} >${rootjump}.r
react ${rootjump} >${rootjump}.r

# and set some variable values
varnames2="`grep time ${rootjump}.r`"
lastline2="`tail -1 ${rootjump}.r`"
firstline2="`cat ${rootjump}.r|awk '/time/ {getline;print;exit}'`"
set $firstline2
first_conc2="$2"
set $lastline2
last_conc2="$2"



# 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 product monitored for halftime calculation
# $4 = final concentration of above

label() {
cat ${1}|awk -v awkvar0=$2 -v awkvar1=$3 -v awkvar2=$4 '
/time/ {
	getline
	t = $1
	c = $2
	while(getline > 0) {
	if (flag == 0 && ($2 - awkvar1)/(awkvar2 - awkvar1) > 0.5) {
		flag = 1
		halftime = t + ($1 - t)*(0.5*(awkvar2 - awkvar1) + awkvar1 - c) / ($2 - c)
	}
	line = $0
	t = $1
	c = $2
	}
}
END	{
	nf = split(line,fields)
	printf "%s:  E=%10.3e, S=%10.3e, (E+S)=%10.3e\n", 
		awkvar0,
		fields[nf-1], fields[nf], fields[nf]+fields[nf-1]
	printf "product t(1/2)=%10.3e, tau=%10.3e, 1/tau=%10.3e\n",
		halftime, halftime/0.691, 0.691/halftime
}' -
}


# PLOT DATA

if [ $output -eq 0 ]
then
	exit
fi

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

if [ $jumponly -eq 0 ]
then

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

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


if [ $nojump -eq 1 ]
then
	exit
fi

if [ $productplot -eq 1 ]
then
	echo "\f"
elif [ $reactionplot -eq 1 ]
then
	echo "\f"
else
	echo "\n"
fi


label ${rootjump}.r "JUMP" $first_conc2 $last_conc2
echo
cat ${rootjump}.r|awk '/time/ {print;getline;print}; {line=$0}; END {print line}'
echo

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

if [ $productplot -eq 1 ]
then
	cat ${rootjump}.r|prep_nologx -c 2 $first_conc2 $last_conc2 $factor |ttygraph
	echo "\n"
fi