/*
The permission is granted to use the software as is and redistribute it in
its original and complete form. If you find this program useful and publish
results obtained by it, please cite the program as:
Michael Whitbeck, Desert Research Institute, 1991, "REACT - Program to solve
kinetic equations for chemical systems". Version 1.00, this code is
available from the authors or from Dr. Whitbeck.
*/
 


This is the REACT package  release 2.


A while back someone asked in sci.chem about chemical kinetics
simulation software.

Well I hacked out a  program that just does the simulation (and parameter
optimization).

It is written in C for the Sun 3|4 computers (should port to other
platforms- I ported an intermediate version to GEMDOS (1040ST) with
no great difficulty).

The code uses  lsoda as the core integrater = livermore solver
for ordinary differential equations, with automatic method switching 
for stiff and nonstiff problems by  Linda R. Petzold  and  Alan C.
Hindmarsh, computing and mathematics research division, l-316
lawrence livermore national laboratory (lsoda.f is available at
netlib@research.att.com, the C port was provided by Hon Wah Tam
Wolfram Research, Inc.  tam@wri.com).

INPUT:
	two files
		root.m where root can be your choice of names
		This file contains a reaction mechanism (free
		format) like:

2.0 	a -> b
0.15 	b -> c 
0.80	b + c -> d + d

the first field is a rate constant in 'F' or 'E' format (e.g. 1.00e-09)
then a list of reactants and products separated by '->' (required).
The specie names can be anything, eg CH3O2, NO2, PABA etc. There
are some restrictions on name lengths (10 chars) and some tokens
are 'special'.

		root.p is a file containing the integration
		parameters again in a free format. It looks like:

0.0  6 .1 1.0E-9
a	1.0	1.0e-10		1.0e-6
b	0.0	1.0e-10		1.0e-8

You simply provide (first line) starting time, stopping time,
output interval, and a solution tolerance value (a default value
for solution accuracy- to be explained in the package).
This line is followed by one or more lines to initialize specie
values (a,b,c,d in the example). No entry is required if the
initial value is zero and the default tolerance is acceptable.
Otherwise give the species name, starting value, absolute
tolerance, relative tolerance (explained in the package).


The program REACT sucks in these two files and spits out a file
named root.r (again root is from your input file name; i.e. 
REACT A produces A.r from A.[m,p]). This file contains the answers!
To facilitate printing and plotting the package contains 3 post
processors:

SELECT	reads root.e, a file containing the species to be printed.
	if root.e contains the lines
	a b d
	c b
	then a table is printed with the headings 
	TIME 	a	b	d
	then a table is printed with the headings
	TIME	c	b
 a sample from the above example might be:
3 reacions involving 4 species
a  b  c  d 
       time              a              b              c
     0.0000e+00      1.0000e+00      0.0000e+00      0.0000e+00 
     1.0000e-01      8.1873e-01      1.7987e-01      1.3925e-03 
     2.0000e-01      6.7032e-01      3.2439e-01      5.1475e-03 
     3.0000e-01      5.4881e-01      4.3989e-01      1.0669e-02 
     4.0000e-01      4.4933e-01      5.3151e-01      1.7433e-02 
     ....
     6.0000e+00      6.1442e-06      2.4471e-01      1.7167e-01 


EXTRACT	same as select but no table headings, useful with unix
	tools like dm, awk, etc.
EXTRACTER	dumps all the data (prints to stdout) in ASCII
	starting with time.


A fifth file root.s contains the list of species that goes with
root.r so that root.m does not have to be reprocessed to make these
tables.

===============================

to make it first compile the lsoda integrater then edit this
makefile to reflect where you are keeping the integrater code then
make react. If that succeeds make select extract extracter.

NOTE: as time goes buy ports may be made to other platforms; check
the makefile for defines for your system!

Send correspondance, bug reports to 

___________________________________________________________
|Mike Whitbeck             | whitbeck@unssun.unr.edu      |
|Desert Research Inst.     | whitbeck@wheeler.wrc.unr.edu |
|POB 60220                 | whitbeck@sanjuan.UUCP        |
|RENO, NV 89506            | 702-673-7348                 |
|__________________________|______________________________|

========================================================================
============================UPDATE======================================
========================================================================
=========================================================================
============additional info for least-squares fitting====================
=========================================================================

REACT can optimize the rate parameters used in the model if it is given a file
root.f containing the measured data for fitting. Add the keyword 'fit' to invoke
this option.


a sample root.f file for react version 2 (pre-release)

(do not enter the <- marker or anything after it)

100            <- max number  of iterations allowed, CAREFUL!
3 3		<- no. parameters to fit    no. observed species to use
1 2 3		<- list of subscripts for k_i's
 a b c		<- list of specie names that are observed
0.01 0.01 0.001	<- step sizes to be used in hunting for the answer
1e-2 1e-2 1e-2 1e-1	<- max. error (fractional) allowed/sought
     0.0000e+00      1.0490e+01      4.1000e+00      0.0000e+00 < d  
     2.0000e+01      7.5018e+00      6.7479e+00      3.4029e-01 < a
     4.0000e+01      6.3005e+00      7.5160e+00      7.7354e-01 < t
     6.0000e+01      5.7445e+00      7.6161e+00      1.2294e+00 < a
     8.0000e+01      5.4238e+00      7.4833e+00      1.6829e+00 < 

the input data for the observables are in the order time and then
concentrations for the measured species in the order given in the list
above.
==========================sample output=======================
equinox:141 src> time react fit tas2 

1) 9.00e-02	a -> b
2) 9.00e-02	b -> a
3) 9.00e-03	b -> c
SIMPLEX Optimization


Program exited after 113 iterations.

The mean is:  3.000023e-02  2.000013e-02  2.999969e-03  5.258171e-09

The estimated fractional error is:  1.747761e-05  2.325398e-05  1.065852e-05  7.440557e-02


The standard deviation is 6.985651e-05

The estimated error of the function is 2.016584e-05

4.030u 0.220s 0:05.61 75.7% 0+235k 3+7io 0pf+0w
==========================sample output=======================
----explanation-----

IN THIS EXAMPLE I specified that k's k_1, k_2, and k_3 were to be optimized
using measured data for [a], [b], and [c]. The initial k's used to
start the model were:  .09, .09, and .009
The model found the best values to be :
  3.000023e-02  2.000013e-02  2.999969e-03
(from the line 'the mean is:')

With a new model first get it running with some rough guesses for the
parameter(s) to be fitted. Next try a loose fit (few iterations low
maxerr criteria) just a few data points closely spaced (e.g. at t= 0,
10, 20, 30, 40). This will give you a better starting point. Improve
the starting guesses for the fitted k's using this output then fit
with as much data as you have.

If convergence is unsatisfactory re-think the model, re-check the data
entry, etc.

Ther may/propbably will be changes and possibly bugs will be
uncovered-  report bugs to 

whitbeck@unssun.unr.edu