diff --git a/.gitignore b/.gitignore
--- a/.gitignore
+++ b/.gitignore
@@ -1,11 +1,10 @@
 syntax: glob
-src/*.o
-src/*.hi
 src/simgi
-src/*.swp
-test/*.o
-test/*.hi
-test/*.swp
+*.o
+*.hi
+*.swp
 average
 RpnStackTest
+EventParserTest
+ReactionParserTest
 check_deviation
diff --git a/ChangeLog b/ChangeLog
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,14 @@
+2009-12-01 Markus Dittrich <haskelladdict@users.sourceforge.net>
+
+      * 0.2 released.
+      * significant overhaul of SGL language spec. This will break
+        backward compatibility with simgi 0.1.
+      * replaced simple RNG with a 64 bit Mersenne Twister implementation.
+      * added expression statements and the ability to describe
+        simulation events via the "events" block.
+      * added ability to specify output format.
+
+
 2009-06-02 Markus Dittrich <haskelladdict@users.sourceforge.net>
 
       * 0.1.1 release (identical to 0.1 with only a few minor
diff --git a/Makefile b/Makefile
--- a/Makefile
+++ b/Makefile
@@ -1,7 +1,7 @@
 # Copyright 2009 Markus Dittrich <haskelladdict@users.sourceforge.net>
 # Distributed under the terms of the GNU General Public License v3
 
-VERSION=0.1.1
+VERSION=0.2
 DESTDIR=
 prefix=/usr
 mandir=$(DESTDIR)$(prefix)/share/man/man1
@@ -17,7 +17,7 @@
 	  src/PrettyPrint.hs src/RpnParser.hs src/RpnCalc.hs \
 	  src/RpnData.hs src/TokenParser.hs 
 
-all: simgi
+all: simgi 
 
 simgi: $(OBJECTS) 
 	ghc -i./Models -i./src $(GHC_FLAGS_RELEASE) --make src/simgi.hs
@@ -33,14 +33,13 @@
 	make -C test
 
 
-
 install: simgi
 	install -d $(docdir)
 	install -d $(bindir)
 	install -d $(htmldir)
 	install -m 0755 src/simgi $(bindir)/
 	install -m 0644 ChangeLog COPYING AUTHORS $(docdir)/
-	install -m 0644 doc/*.pdf doc/*.rst $(docdir)/
+	install -m 0644 doc/*.pdf $(docdir)/
 	install -m 0644 doc/*.html $(htmldir)/
 
 
diff --git a/Models/brusselator.sgl b/Models/brusselator.sgl
--- a/Models/brusselator.sgl
+++ b/Models/brusselator.sgl
@@ -6,21 +6,25 @@
 ------------------------------------------------------}
 
 def parameters
-  time       = 50.0
-  outputIter = 50000
-  outputFreq = 200
-  systemVol  = nil  -- interpret rates as propensities
-  outputFile = "brusselator_output.dat"
+  time         = 50.0
+  outputBuffer = 50000
+  outputFreq   = 200
+  systemVol    = nil  -- interpret rates as propensities
+  outputFile   = "brusselator_output.dat"
 end
 
 def molecules
-  x 1000
-  y 2000
+  x = 1000
+  y = 2000
 end
 
 def reactions
-  nil    -> x         { 5000    }
-  x      -> y         { 50.0    }
-  2x + y -> 3x        { 0.00005 }
-  x      -> nil       { 5.0     }
-end 
+  nil    -> x         | 5000    |
+  x      -> y         | 50.0    |
+  2x + y -> 3x        | 0.00005 |
+  x      -> nil       | 5.0     |
+end
+
+def output
+  [x,y]
+end
diff --git a/Models/oregonator.sgl b/Models/oregonator.sgl
--- a/Models/oregonator.sgl
+++ b/Models/oregonator.sgl
@@ -6,23 +6,27 @@
 ------------------------------------------------------}
 
 def parameters
-  time       = 50.0
-  outputIter = 50000
-  outputFreq = 500
-  systemVol  = nil  -- interpret rates as propensities
-  outputFile = "oregonator_output.dat"
+  time         = 50.0
+  outputBuffer = 50000
+  outputFreq   = 500
+  systemVol    = nil  -- interpret rates as propensities
+  outputFile   = "oregonator_output.dat"
 end
 
 def molecules
-  x 500
-  y 1000
-  z 2000  
+  x = 500
+  y = 1000
+  z = 2000  
 end
 
 def reactions
-  y      -> x         { 2.0   }
-  x + y  -> nil       { 0.1   }
-  x      -> 2x + z    { 104.0 }
-  2x     -> nil       { 0.016 }
-  z      -> y         { 26.0  }
-end 
+  y      -> x         | 2.0   |
+  x + y  -> nil       | 0.1   |
+  x      -> 2x + z    | 104.0 |
+  2x     -> nil       | 0.016 |
+  z      -> y         | 26.0  |
+end
+
+def output
+  [x,y,z]
+end
diff --git a/Models/volterra.sgl b/Models/volterra.sgl
--- a/Models/volterra.sgl
+++ b/Models/volterra.sgl
@@ -6,21 +6,24 @@
 ------------------------------------------------------}
 
 def parameters
-  time       = 50.0
-  outputIter = 50000
-  outputFreq = 200
-  systemVol  = nil  -- interpret rates as propensities
-  outputFile = "volterra_output.dat"
+  time         = 50.0
+  outputBuffer = 50000
+  outputFreq   = 200
+  systemVol    = nil  -- interpret rates as propensities
+  outputFile   = "volterra_output.dat"
 end
 
 def molecules
-  x 1000
-  y 2000
+  x = 1000
+  y = 2000
 end
 
 def reactions
-  x      -> 2x       { 10.0 }
-  x + y  -> 2y       { 0.01 }
-  y      -> nil      { 10.0 }
+  x      -> 2x       | 10.0 |
+  x + y  -> 2y       | 0.01 |
+  y      -> nil      | 10.0 |
 end 
 
+def output
+  [x,y]
+end
diff --git a/README b/README
--- a/README
+++ b/README
@@ -1,6 +1,7 @@
 DESCRIPTION:
 ------------
 
+
 simgi is a stochastic simulator using the Gillespie algorithm [1,2].
 simgi is released under the GPL v3.
 
diff --git a/doc/Makefile b/doc/Makefile
--- a/doc/Makefile
+++ b/doc/Makefile
@@ -21,11 +21,12 @@
 
 
 clean_old:
-	rm -f *.html *.pdf
+	rm -f *.html *.pdf 
 
 
 .PHONY: clean
 
 clean:
+	rm -f *.aux *.log *.html *.pdf *.tex *.out
 	
 	
diff --git a/doc/simgi.html b/doc/simgi.html
--- a/doc/simgi.html
+++ b/doc/simgi.html
@@ -3,14 +3,14 @@
 <html xmlns="http://www.w3.org/1999/xhtml" xml:lang="en" lang="en">
 <head>
 <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
-<meta name="generator" content="Docutils 0.5: http://docutils.sourceforge.net/" />
+<meta name="generator" content="Docutils 0.6: http://docutils.sourceforge.net/" />
 <title>simgi - A Stochastic Gillespie Simulator for Molecular Systems</title>
 <meta name="author" content="Markus Dittrich" />
 <style type="text/css">
 
 /*
 :Author: David Goodger (goodger@python.org)
-:Id: $Id: html4css1.css 5196 2007-06-03 20:25:28Z wiemann $
+:Id: $Id: html4css1.css 5951 2009-05-18 18:03:10Z milde $
 :Copyright: This stylesheet has been placed in the public domain.
 
 Default cascading style sheet for the HTML output of Docutils.
@@ -158,12 +158,33 @@
 hr.docutils {
   width: 75% }
 
-img.align-left {
-  clear: left }
+img.align-left, .figure.align-left{
+  clear: left ;
+  float: left ;
+  margin-right: 1em }
 
-img.align-right {
-  clear: right }
+img.align-right, .figure.align-right {
+  clear: right ;
+  float: right ;
+  margin-left: 1em }
 
+.align-left {
+  text-align: left }
+
+.align-center {
+  clear: both ;
+  text-align: center }
+
+.align-right {
+  text-align: right }
+
+/* reset inner alignment in figures */
+div.align-right {
+  text-align: left }
+
+/* div.align-center * { */
+/*   text-align: left } */
+
 ol.simple, ul.simple {
   margin-bottom: 1em }
 
@@ -217,8 +238,7 @@
 pre.address {
   margin-bottom: 0 ;
   margin-top: 0 ;
-  font-family: serif ;
-  font-size: 100% }
+  font: inherit }
 
 pre.literal-block, pre.doctest-block {
   margin-left: 2em ;
@@ -296,7 +316,7 @@
 <tr class="field"><th class="docinfo-name">email:</th><td class="field-body">haskelladdict at users dot sourceforge dot net</td>
 </tr>
 <tr><th class="docinfo-name">Version:</th>
-<td>0.1.1 (06/02/2009)</td></tr>
+<td>0.2 (12/01/2009)</td></tr>
 </tbody>
 </table>
 <div class="section" id="contents">
@@ -309,28 +329,22 @@
 <li><a class="reference internal" href="#simgi-model-generation-language-sgl">Simgi Model Generation Language (SGL)</a></li>
 <li><a class="reference internal" href="#example-input-files">Example Input Files</a></li>
 <li><a class="reference internal" href="#bugs">Bugs</a></li>
-<li><a class="reference internal" href="#references">References</a></li>
 </ol>
 </div>
 <div class="section" id="introduction">
 <h1>Introduction</h1>
 <p><strong>simgi</strong> is a fairly simple and straightforward stochastic simulator
-based on Gillspie's <a class="footnote-reference" href="#id7" id="id1">[1]</a> direct method. <strong>simgi</strong> is implemented in
+based on Gillespie's <a class="footnote-reference" href="#id4" id="id1">[1]</a> direct method. <strong>simgi</strong> is implemented in
 pure Haskell, command line driven and comes with a flexible simulation
 description language called <a class="reference internal" href="#simgi-model-generation-language-sgl">Simgi Model Generation Language (SGL)</a>.
-More information is available from the <a class="reference external" href="http://sourceforge.net/projects/simgi">project summary page</a>.</p>
+<strong>simgi</strong> uses a fast 64 bit implementation of the Mersenne Twister
+algorithm as random number source.</p>
 </div>
 <div class="section" id="status">
 <h1>Status</h1>
-<p>The 0.1 release of <strong>simgi</strong> provides a fully functional simulator
-but should still be treated as an alpha version since several parts
-of the code are currently not fully optimal. This is particularly
-true for the random number generator which presently leverages the
-StdGen instance of RandomGen <a class="footnote-reference" href="#id8" id="id2">[2]</a> and is probably not sufficient for
-large simulations in terms of random number quality. Later revisions
-of simgi will have a more sophisticated random number generator.
-Nevertheless, for small systems (such as the examples in the
-<em>Models/</em> directory) the current implementation should be sufficient.</p>
+<p>The 0.2 release of <strong>simgi</strong> provides a fully functional simulator
+which has been tested on several model systems some of which were
+fairly large.</p>
 </div>
 <div class="section" id="download">
 <h1>Download</h1>
@@ -342,6 +356,7 @@
 <ul class="simple">
 <li><a class="reference external" href="http://haskell.org/ghc/">&gt;=ghc-6.10</a></li>
 <li><a class="reference external" href="http://gmplib.org/">&gt;=gmp-4.3</a></li>
+<li><a class="reference external" href="http://hackage.haskell.org/package/mersenne-random-pure64">&gt;=mersenne-random-pure64</a></li>
 </ul>
 <p>To compile the documentation (not required), you will also need</p>
 <ul class="simple">
@@ -350,15 +365,16 @@
 </ul>
 <p>Building of <strong>simgi</strong> can be done either via</p>
 <ul class="simple">
-<li>the standard <tt class="docutils literal"><span class="pre">make,</span> <span class="pre">make</span> <span class="pre">check,</span> <span class="pre">make</span> <span class="pre">install</span></tt></li>
+<li>the standard <tt class="docutils literal">make, make check, make install</tt></li>
 <li>or via cabal</li>
 </ul>
 </div>
 <div class="section" id="simgi-model-generation-language-sgl">
 <h1>Simgi Model Generation Language (SGL)</h1>
-<p>simgi simulations are described via <a class="reference internal" href="#simgi-model-generation-language-sgl">Simgi Model Generation Language
-(SGL)</a>. The corresponding simulation files typically have an <em>.sgl</em>
-extension, but this is not enforced by the <strong>simgi</strong> simulation engine.</p>
+<p><strong>simgi</strong> simulations are described via <a class="reference internal" href="#simgi-model-generation-language-sgl">Simgi Model Generation Language
+(SGL)</a>. The corresponding simulation input files typically have an <em>.sgl</em>
+extension, but this is not enforced by the <strong>simgi</strong> simulation
+engine.</p>
 <p>A SGL file consists of zero or more descriptor blocks of the form</p>
 <pre class="literal-block">
 def &lt;block name&gt;
@@ -367,120 +383,190 @@
 
 end
 </pre>
-<p>The (but see <a class="footnote-reference" href="#id9" id="id3">[3]</a>) formatting of the input files is very flexible. In
-particular, neither newlines <a class="footnote-reference" href="#id10" id="id4">[4]</a> nor extraneous whitespace matter.
-Hence, the above SDL block could also be written on a single line.
+<p>The formatting of the input files is very flexible. In
+particular, neither newlines <a class="footnote-reference" href="#id5" id="id2">[2]</a> nor extraneous whitespace matter.
+Hence, the above SGL block could have also been written on a single line.
 However, it is strongly recommended to stick to a consistent and
 &quot;visually simple&quot; layout to aid in &quot;comprehending&quot; the underlying
-model.</p>
+model. Also, it is important to point out that <strong>simgi</strong>'s parser is
+case sensitive.</p>
 <p><strong>Comments</strong> can be added to the SGL file and are parsed according to
 the Haskell language specs</p>
 <ul class="simple">
 <li>simple line comments begin with a <tt class="docutils literal"><span class="pre">--</span></tt> token and treat everything
-until the next newline as a comment, including valid SDL commands.
-Hence, SDL blocks containing line comments need to be separated by
+until the next newline as a comment, including valid SGL commands.
+Hence, SGL blocks containing line comments need to be separated by
 newlines in order to be parsed correctly.</li>
-<li>block comments begin with a <tt class="docutils literal"><span class="pre">{-</span></tt> token and end with a <tt class="docutils literal"><span class="pre">-}</span></tt>
+<li>block comments begin with a <tt class="docutils literal">{-</tt> token and end with a <tt class="docutils literal"><span class="pre">-}</span></tt>
 token. Everything within a comment block is ignored by the parser
 and block comments can be nested.</li>
 </ul>
-<p>Currently, the SDL specs define the following block types with their
+<p><strong>Expression Statements</strong> are an important and useful part of SGL.
+<tt class="docutils literal">Expression statements</tt> are enclosed in curly braces and can contain
+any mathematical expression involving doubles, the simulation time
+(via the keyword <tt class="docutils literal">TIME</tt>), as well as the values of any variable or
+molecule count. The values of time, molecule counts and variables
+are evaluated at run time and represent the instantaneous values at the
+time at which the expression is evaluated.
+<tt class="docutils literal">Expressions statements</tt> can contain any
+arithmetic expression involving the standard operators &quot;+&quot;, &quot;-&quot;, &quot;*&quot;, &quot;/&quot;, &quot;^&quot;
+(exponentiation), and the mathematical functions <tt class="docutils literal">sqrt, exp, log, log2, log10, sin,
+cos, tan, asin, acos, atan, sinh, cosh, tanh, asinh, acosh, atanh, acosh, atanh,
+erf, erfc, abs</tt>.</p>
+<p>Internally, <tt class="docutils literal">expression statements</tt> are converted into a compute stack
+in RPN format which is evaluated at run-time. Even though this
+procedure is fairly efficient, there is some numerical overhead
+incurred at each iteration and the use of complicated rate
+expressions should therefore be avoided if possible.</p>
+<p>Below is a list of all SGL blocks available for describing simulations.
+Presently, the order of blocks matters and should be exactly the same
+in which they are described below. Several SGL blocks are
+optional and are marked as such below.</p>
+<p>Currently, the SGL specs define the following block types with their
 respective block commands and block content:</p>
-<p><strong>parameter block:</strong> <em>&lt;block name&gt; = parameters</em></p>
+<p><strong>parameter block:</strong> <tt class="docutils literal">&lt;block name&gt;</tt> = <em>parameters</em></p>
 <blockquote>
 <p>The purpose of the parameter block is to describe the global
 simulation parameters. The following parameters are currently
 supported:</p>
 <dl class="docutils">
-<dt><em>time = &lt;double&gt;</em></dt>
+<dt><em>time</em> = <tt class="docutils literal">Double</tt></dt>
 <dd>Maximum simulation time in seconds. Default is 0.0 s.</dd>
-<dt><em>outputIter = &lt;Integer&gt;</em></dt>
+<dt><em>outputBuffer</em> = <tt class="docutils literal">Integer</tt></dt>
 <dd><p class="first">Output will be kept in memory and written to the output file and
-stdout every <em>outputIter</em> iterations. Larger values should result
-in faster simulations but require more system memory.
+stdout every <em>outputBuffer</em> iterations. Larger values should
+result in faster simulations but require more system memory.
 Default is to write output every 10000 iterations.</p>
-<p class="last">Note: <em>outputIter</em> only affects how often output is written to
-the output file, not how much is being accumulated during a
+<p class="last">Note: <em>outputBuffer</em> only affects how often output is written to
+the output file, not how much output is actually generated during a
 simulation (see outputFreq parameter).</p>
 </dd>
-<dt><em>outputFreq = &lt;Integer&gt;</em></dt>
-<dd>Frequency with which output is generated. Default is 1000.</dd>
-<dt><em>systemVol = &lt;double&gt;</em></dt>
+<dt><em>outputFreq</em> = <tt class="docutils literal">Integer</tt></dt>
+<dd>Iteration frequency with which output is generated. Default is every 1000
+iterations. Please note that output is written to the output file in batches of
+<em>outputBuffer</em>.</dd>
+<dt><em>systemVol</em> = <tt class="docutils literal">Double</tt></dt>
 <dd>Volume of the simulation system in liters. This is needed to
 properly compute the reaction rates in molar units. If rates
 should rather be interpreted as reaction propensities (like
-in <a class="footnote-reference" href="#id7" id="id5">[1]</a>) please set <em>systemVol = nil</em>. Default is a system
+in <a class="footnote-reference" href="#id4" id="id3">[1]</a>) please set <em>systemVol = nil</em>. Default is a system
 volume of 1.0 liter.</dd>
-<dt><em>outputFile = &lt;quoted string&gt;</em>:</dt>
+<dt><em>outputFile</em> = <tt class="docutils literal">Quoted String</tt></dt>
 <dd>Name of the output file. This is the only required parameter
 in the parameter section. If not given, the simulation will
 terminate.</dd>
 </dl>
 </blockquote>
-<p><strong>molecule block:</strong> <em>&lt;block name&gt; = molecules</em></p>
+<p><strong>variable block:</strong> <tt class="docutils literal">&lt;block name&gt;</tt> = <em>variables</em></p>
 <blockquote>
 <p>This block consist of a list of pairs of the form</p>
 <pre class="literal-block">
-&lt;String&gt; &lt;Integer&gt;
+String = &lt;variable expression&gt;
 </pre>
+<p>where <tt class="docutils literal">String</tt> is the variable name, and <tt class="docutils literal">&lt;variable expression&gt;</tt>
+is either a <tt class="docutils literal">Double</tt> or an <tt class="docutils literal">expression statement</tt> as defined above.
+Variables can be used in any other <tt class="docutils literal">expression statement</tt> in the
+SGL file including reaction rate definitions. Please make sure to
+not define a variable in terms of itself to avoid infinite recursion.</p>
+</blockquote>
+<p><strong>molecule block:</strong> <tt class="docutils literal">&lt;block name&gt;</tt> = <em>molecules</em></p>
+<blockquote>
+<p>This block consist of a list of pairs of the form</p>
+<pre class="literal-block">
+String = Integer
+</pre>
 <p>giving the name of each molecule and the number of molecules
 present initially. For example, the following molecule definition
-block defines molecules <tt class="docutils literal"><span class="pre">A</span></tt> and <tt class="docutils literal"><span class="pre">B</span></tt> with initial numbers of
+block defines molecules <tt class="docutils literal">A</tt> and <tt class="docutils literal">B</tt> with initial numbers of
 100 and 200, respectively</p>
 <pre class="literal-block">
 def molecules
-  A 100
-  B 200
+  A = 100
+  B = 200
 end
 </pre>
+<p><strong>NOTE</strong>: Please do not use any of the predefined mathematical
+functions or defined variables (including <tt class="docutils literal">TIME</tt>) as
+molecule names since this will lead to undefined behavior.</p>
 </blockquote>
-<p><strong>reaction block</strong>: <em>&lt;block name&gt; = reactions</em></p>
+<p><strong>reaction block</strong>: <tt class="docutils literal">&lt;block name&gt;</tt> = <em>reactions</em></p>
 <blockquote>
 <p>This block describes the reactions between molecules defined in
 the molecule block. Reactions are specified via</p>
 <pre class="literal-block">
-reactants -&gt; product { rate expression }
+&lt;reactants&gt; -&gt; &lt;products&gt;  | &lt;rate expression&gt; |
 </pre>
-<p>Here, <tt class="docutils literal"><span class="pre">reactants</span></tt> and <tt class="docutils literal"><span class="pre">products</span></tt> are of the form</p>
+<p>Here, <tt class="docutils literal">&lt;reactants&gt;</tt> and <tt class="docutils literal">&lt;products&gt;</tt> are of the form</p>
 <pre class="literal-block">
-&lt;Integer&gt; &lt;String&gt; + &lt;Integer&gt; &lt;String&gt; + .....
+Integer String + Integer String + .....
 </pre>
-<p>In this expression, <tt class="docutils literal"><span class="pre">&lt;String&gt;</span></tt> is the reactant or product name
-as defined in the molecule block and <tt class="docutils literal"><span class="pre">&lt;Integer&gt;</span></tt> an optional
-integer specifying the stoichiometry. If <tt class="docutils literal"><span class="pre">&lt;Integer&gt;</span></tt> is not
+<p>In this expression, <tt class="docutils literal">String</tt> is a molecule name
+as defined in the molecule block and <tt class="docutils literal">Integer</tt> an optional
+integer specifying the stoichiometry. If <tt class="docutils literal">Integer</tt> is not
 explicitly given, it is assumed to be 1.</p>
-<p>The reaction rate can either be a fixed value of type <tt class="docutils literal"><span class="pre">&lt;Double&gt;</span></tt>
-or else an mathematical expression involving <tt class="docutils literal"><span class="pre">&lt;Double&gt;</span></tt>,
-molecule names, and the current simulation time. Hence, <strong>simgi</strong>
-rate expressions can be arbitrary complex functions of the
-instantaneous simulation time and the instantaneous numbers of any
-molecule in the model. The parser will interpret any string in the
-rate expression as a molecule name in a case sensitive fashion,
-a mathematical operator or function (see <a class="footnote-reference" href="#id11" id="id6">[5]</a> for supported
-functions), or the special variable TIME which refers to the
-current simulation time. Hence, do <strong>not</strong> use any of the
-mathematical keywords as a molecule name; this leads to undefined
-behavior.</p>
-<p>Here is an example reaction block for the two molecules <tt class="docutils literal"><span class="pre">A</span></tt> and
-<tt class="docutils literal"><span class="pre">B</span></tt> defined above:</p>
+<p>The <tt class="docutils literal">&lt;rate expression&gt;</tt> can either be a fixed value of type
+<tt class="docutils literal">Double</tt> or an <tt class="docutils literal">expression statement</tt> as defined above.</p>
+<p>Below is an example reaction block for the two molecules <tt class="docutils literal">A</tt> and
+<tt class="docutils literal">B</tt> defined above:</p>
 <pre class="literal-block">
 define reactions
-  2A + B -&gt; A  { 10.0e-5 }
-  B      -&gt; A  { 2.0e-5 * A * exp(-0.5*TIME) }
+  2A + B -&gt; A  | 10.0e-5 |
+  B      -&gt; A  | {2.0e-5 * A * exp(-0.5*TIME)} |
 end
 </pre>
-<p>In the first reaction, 2 <tt class="docutils literal"><span class="pre">A</span></tt> molecules react with one <tt class="docutils literal"><span class="pre">B</span></tt> to
-yield another <tt class="docutils literal"><span class="pre">A</span></tt> at a rate of 10.0e-5 1/(Mol s). The second
-reaction describes a decay of <tt class="docutils literal"><span class="pre">B</span></tt> back to <tt class="docutils literal"><span class="pre">A</span></tt> at a rate
-that is computed based on the instantaneous number of <tt class="docutils literal"><span class="pre">A</span></tt>
+<p>In the first reaction, 2 <tt class="docutils literal">A</tt> molecules react with one <tt class="docutils literal">B</tt> to
+yield another <tt class="docutils literal">A</tt> at a rate of 10.0e-5. The second
+reaction describes a decay of <tt class="docutils literal">B</tt> back to <tt class="docutils literal">A</tt> at a rate
+that is computed based on the instantaneous number of <tt class="docutils literal">A</tt>
 molecules present and which decays exponentially with simulation
 time.</p>
-<p>Internally, rate expressions are converted into a compute stack
-in RPN format which is evaluated at run-time. Even though this
-procedure is fairly efficient, there is some numerical overhead
-incurred at each iteration and the use of complicated rate
-expressions should therefore be avoided if possible.</p>
 </blockquote>
+<p><strong>event block</strong>: <tt class="docutils literal">&lt;block name&gt;</tt> = <em>events</em></p>
+<blockquote>
+<p>An event block allows one to specify events which will occur during
+the simulation. Each event consists of a <tt class="docutils literal">&lt;trigger expression&gt;</tt> and
+an associated set of <tt class="docutils literal">&lt;action expressions&gt;</tt>.
+Events are specified via</p>
+<pre class="literal-block">
+{ &lt;trigger expression&gt; } =&gt; { &lt;action expression&gt; }
+</pre>
+<p>Here, <tt class="docutils literal">trigger expression</tt> is of the form</p>
+<pre class="literal-block">
+&lt;trigger primitive&gt; [ &lt;boolean operator&gt; &lt;trigger primitive&gt;]
+</pre>
+<p>with <tt class="docutils literal">&lt;trigger primitive&gt;</tt> defined by</p>
+<pre class="literal-block">
+&lt;expression statement&gt; relational operator &lt;expression statement&gt;
+</pre>
+<p>Each <tt class="docutils literal">&lt;trigger primitive&gt;</tt> contains two <tt class="docutils literal">expression statements</tt>
+as defined above and a <tt class="docutils literal">relational operator</tt> which can be
+any of <tt class="docutils literal">&gt;=</tt>, <tt class="docutils literal">&lt;=</tt>, <tt class="docutils literal">==</tt>, <tt class="docutils literal">&gt;</tt>, and <tt class="docutils literal">&lt;</tt>. Hence, each
+<tt class="docutils literal">&lt;trigger primitive&gt;</tt> evaluates to either <tt class="docutils literal">true</tt> or <tt class="docutils literal">false</tt>.</p>
+<p>Several <tt class="docutils literal">&lt;trigger primitives&gt;</tt> can be chained together via the
+<tt class="docutils literal">&lt;boolean operators&gt;</tt> <tt class="docutils literal">&amp;&amp;</tt> and <tt class="docutils literal">||</tt> to yield a final boolean
+value of <tt class="docutils literal">true</tt> or <tt class="docutils literal">false</tt>.</p>
+<p>If the <tt class="docutils literal">&lt;trigger expression&gt;</tt> evaluates to true during an
+iteration, the associated <tt class="docutils literal">&lt;action expressions&gt;</tt> is executed
+during the same timestep.</p>
+<p><tt class="docutils literal">&lt;action expression&gt;</tt> consists of a semi-colon separated list of
+assignments</p>
+<pre class="literal-block">
+String = &lt;assignment expression&gt; [; String = &lt;assignment expression&gt;]
+</pre>
+<p>where <tt class="docutils literal">String</tt> is a molecule or variable name and
+<tt class="docutils literal">&lt;expression&gt;</tt> either a <tt class="docutils literal">Double</tt> or an <tt class="docutils literal">expression statement</tt>.</p>
+<p><strong>NOTE</strong>: Since molecule counts are integer values assignments
+to molecule counts in <tt class="docutils literal">&lt;action expression&gt;</tt> will be converted
+to an integer value via <tt class="docutils literal">floor</tt>.</p>
+</blockquote>
+<p><strong>output block</strong>: <tt class="docutils literal">&lt;block name&gt;</tt> = <em>output</em></p>
+<blockquote>
+<p>This block consists of a simple list of variable and molecule
+names that will be streamed to the output file in the same order:</p>
+<pre class="literal-block">
+[ name1, name2, name3, .... ]
+</pre>
+</blockquote>
 </div>
 <div class="section" id="example-input-files">
 <h1>Example Input Files</h1>
@@ -496,37 +582,16 @@
 <h1>Bugs</h1>
 <p>Please report all bugs and feature requests to
 &lt;haskelladdict at users dot sourceforge dot net&gt;.</p>
-</div>
-<div class="section" id="references">
-<h1>References</h1>
-<table class="docutils footnote" frame="void" id="id7" rules="none">
-<colgroup><col class="label" /><col /></colgroup>
-<tbody valign="top">
-<tr><td class="label">[1]</td><td><em>(<a class="fn-backref" href="#id1">1</a>, <a class="fn-backref" href="#id5">2</a>)</em> Daniel T. Gillespie (1977). &quot;Exact Stochastic Simulation of Coupled Chemical Reactions&quot;. The Journal of Physical Chemistry 81 (25): 2340-2361</td></tr>
-</tbody>
-</table>
-<table class="docutils footnote" frame="void" id="id8" rules="none">
-<colgroup><col class="label" /><col /></colgroup>
-<tbody valign="top">
-<tr><td class="label"><a class="fn-backref" href="#id2">[2]</a></td><td><a class="reference external" href="http://hackage.haskell.org/packages/archive/random/1.0.0.1/doc/html/System-Random#globalrng.html">http://hackage.haskell.org/packages/archive/random/1.0.0.1/doc/html/System-Random#globalrng.html</a></td></tr>
-</tbody>
-</table>
-<table class="docutils footnote" frame="void" id="id9" rules="none">
-<colgroup><col class="label" /><col /></colgroup>
-<tbody valign="top">
-<tr><td class="label"><a class="fn-backref" href="#id3">[3]</a></td><td>Since <strong>simgi</strong> currently is an alpha version there may be fairly drastic changes to the SDL specs in future releases until the first beta release.</td></tr>
-</tbody>
-</table>
-<table class="docutils footnote" frame="void" id="id10" rules="none">
+<table class="docutils footnote" frame="void" id="id4" rules="none">
 <colgroup><col class="label" /><col /></colgroup>
 <tbody valign="top">
-<tr><td class="label"><a class="fn-backref" href="#id4">[4]</a></td><td>An exception to this rule are line comments starting with <tt class="docutils literal"><span class="pre">--</span></tt> which ingnore everything until the next newline.</td></tr>
+<tr><td class="label">[1]</td><td><em>(<a class="fn-backref" href="#id1">1</a>, <a class="fn-backref" href="#id3">2</a>)</em> Daniel T. Gillespie (1977). &quot;Exact Stochastic Simulation of Coupled Chemical Reactions&quot;. The Journal of Physical Chemistry 81 (25): 2340-2361</td></tr>
 </tbody>
 </table>
-<table class="docutils footnote" frame="void" id="id11" rules="none">
+<table class="docutils footnote" frame="void" id="id5" rules="none">
 <colgroup><col class="label" /><col /></colgroup>
 <tbody valign="top">
-<tr><td class="label"><a class="fn-backref" href="#id6">[5]</a></td><td>Rate expressions can contain any arithmetic expression involving the standard operators &quot;+&quot;, &quot;-&quot;, &quot;*&quot;, &quot;/&quot;, &quot;^&quot; (exponentiation), and the mathematical functions <tt class="docutils literal"><span class="pre">sqrt,</span> <span class="pre">exp,</span> <span class="pre">log,</span> <span class="pre">log2,</span> <span class="pre">log10,</span> <span class="pre">sin,</span> <span class="pre">cos,</span> <span class="pre">tan,</span> <span class="pre">asin,</span> <span class="pre">acos,</span> <span class="pre">atan,</span> <span class="pre">sinh,</span> <span class="pre">cosh,</span> <span class="pre">tanh,</span> <span class="pre">asinh,</span> <span class="pre">acosh,</span> <span class="pre">atanh,</span> <span class="pre">acosh,</span> <span class="pre">atanh,</span> <span class="pre">erf,</span> <span class="pre">erfc,</span> <span class="pre">abs</span></tt>.</td></tr>
+<tr><td class="label"><a class="fn-backref" href="#id2">[2]</a></td><td>An exception to this rule are line comments starting with <tt class="docutils literal"><span class="pre">--</span></tt> which ingnore everything until the next newline.</td></tr>
 </tbody>
 </table>
 </div>
diff --git a/doc/simgi.pdf b/doc/simgi.pdf
Binary files a/doc/simgi.pdf and b/doc/simgi.pdf differ
diff --git a/doc/simgi.rst b/doc/simgi.rst
--- a/doc/simgi.rst
+++ b/doc/simgi.rst
@@ -6,7 +6,7 @@
 
 :email: haskelladdict at users dot sourceforge dot net
 
-:Version: 0.1.1 (06/02/2009)
+:Version: 0.2 (12/01/2009)
 
 
 Contents
@@ -19,30 +19,25 @@
 5) `Simgi Model Generation Language (SGL)`_
 6) `Example Input Files`_
 7) Bugs_
-8) References_
 
+
 Introduction
 ------------
 
 **simgi** is a fairly simple and straightforward stochastic simulator 
-based on Gillspie's [1]_ direct method. **simgi** is implemented in 
+based on Gillespie's [1]_ direct method. **simgi** is implemented in 
 pure Haskell, command line driven and comes with a flexible simulation
 description language called `Simgi Model Generation Language (SGL)`_.
-More information is available from the `project summary page <http://sourceforge.net/projects/simgi>`_.
+**simgi** uses a fast 64 bit implementation of the Mersenne Twister
+algorithm as random number source.
 
 
 Status 
 ------
 
-The 0.1 release of **simgi** provides a fully functional simulator 
-but should still be treated as an alpha version since several parts 
-of the code are currently not fully optimal. This is particularly 
-true for the random number generator which presently leverages the 
-StdGen instance of RandomGen [2]_ and is probably not sufficient for 
-large simulations in terms of random number quality. Later revisions 
-of simgi will have a more sophisticated random number generator. 
-Nevertheless, for small systems (such as the examples in the 
-*Models/* directory) the current implementation should be sufficient. 
+The 0.2 release of **simgi** provides a fully functional simulator 
+which has been tested on several model systems some of which were
+fairly large. 
 
 
 Download
@@ -58,6 +53,7 @@
 
 - `>=ghc-6.10 <http://haskell.org/ghc/>`_
 - `>=gmp-4.3  <http://gmplib.org/>`_  
+- `>=mersenne-random-pure64 <http://hackage.haskell.org/package/mersenne-random-pure64>`_
 
 To compile the documentation (not required), you will also need
 
@@ -75,9 +71,10 @@
 Simgi Model Generation Language (SGL)
 -------------------------------------
 
-simgi simulations are described via `Simgi Model Generation Language 
-(SGL)`_. The corresponding simulation files typically have an *.sgl* 
-extension, but this is not enforced by the **simgi** simulation engine. 
+**simgi** simulations are described via `Simgi Model Generation Language 
+(SGL)`_. The corresponding simulation input files typically have an *.sgl* 
+extension, but this is not enforced by the **simgi** simulation 
+engine. 
 
 A SGL file consists of zero or more descriptor blocks of the form
 
@@ -89,132 +86,218 @@
 
   end
 
-The (but see [3]_) formatting of the input files is very flexible. In
-particular, neither newlines [4]_ nor extraneous whitespace matter. 
-Hence, the above SDL block could also be written on a single line. 
+The formatting of the input files is very flexible. In
+particular, neither newlines [2]_ nor extraneous whitespace matter. 
+Hence, the above SGL block could have also been written on a single line. 
 However, it is strongly recommended to stick to a consistent and 
 "visually simple" layout to aid in "comprehending" the underlying
-model.
+model. Also, it is important to point out that **simgi**'s parser is 
+case sensitive.
 
 **Comments** can be added to the SGL file and are parsed according to 
 the Haskell language specs
 
 - simple line comments begin with a ``--`` token and treat everything 
-  until the next newline as a comment, including valid SDL commands. 
-  Hence, SDL blocks containing line comments need to be separated by 
+  until the next newline as a comment, including valid SGL commands. 
+  Hence, SGL blocks containing line comments need to be separated by 
   newlines in order to be parsed correctly.
 - block comments begin with a ``{-`` token and end with a ``-}`` 
   token. Everything within a comment block is ignored by the parser 
   and block comments can be nested.
 
-Currently, the SDL specs define the following block types with their 
+**Expression Statements** are an important and useful part of SGL.
+``Expression statements`` are enclosed in curly braces and can contain
+any mathematical expression involving doubles, the simulation time 
+(via the keyword ``TIME``), as well as the values of any variable or 
+molecule count. The values of time, molecule counts and variables
+are evaluated at run time and represent the instantaneous values at the
+time at which the expression is evaluated.
+``Expressions statements`` can contain any 
+arithmetic expression involving the standard operators "+", "-", "*", "/", "^" 
+(exponentiation), and the mathematical functions ``sqrt, exp, log, log2, log10, sin, 
+cos, tan, asin, acos, atan, sinh, cosh, tanh, asinh, acosh, atanh, acosh, atanh, 
+erf, erfc, abs``.
+
+Internally, ``expression statements`` are converted into a compute stack
+in RPN format which is evaluated at run-time. Even though this
+procedure is fairly efficient, there is some numerical overhead
+incurred at each iteration and the use of complicated rate 
+expressions should therefore be avoided if possible.
+
+
+Below is a list of all SGL blocks available for describing simulations.
+Presently, the order of blocks matters and should be exactly the same
+in which they are described below. Several SGL blocks are 
+optional and are marked as such below. 
+
+Currently, the SGL specs define the following block types with their 
 respective block commands and block content:
 
 
-**parameter block:** *<block name> = parameters*
+**parameter block:** ``<block name>`` = *parameters* 
 
   The purpose of the parameter block is to describe the global 
   simulation parameters. The following parameters are currently
   supported:
 
-  *time = <double>*
+  *time* = ``Double``
     Maximum simulation time in seconds. Default is 0.0 s.
 
-  *outputIter = <Integer>*
+  *outputBuffer* = ``Integer``
     Output will be kept in memory and written to the output file and 
-    stdout every *outputIter* iterations. Larger values should result 
-    in faster simulations but require more system memory. 
+    stdout every *outputBuffer* iterations. Larger values should 
+    result in faster simulations but require more system memory. 
     Default is to write output every 10000 iterations.
 
-    Note: *outputIter* only affects how often output is written to 
-    the output file, not how much is being accumulated during a 
+    Note: *outputBuffer* only affects how often output is written to 
+    the output file, not how much output is actually generated during a 
     simulation (see outputFreq parameter).
 
-  *outputFreq = <Integer>*
-    Frequency with which output is generated. Default is 1000.
+  *outputFreq* = ``Integer``
+    Iteration frequency with which output is generated. Default is every 1000
+    iterations. Please note that output is written to the output file in batches of
+    *outputBuffer*.
 
-  *systemVol = <double>*
+  *systemVol* = ``Double``
     Volume of the simulation system in liters. This is needed to 
     properly compute the reaction rates in molar units. If rates 
     should rather be interpreted as reaction propensities (like 
     in [1]_) please set *systemVol = nil*. Default is a system
     volume of 1.0 liter.
 
-  *outputFile = <quoted string>*:
+  *outputFile* = ``Quoted String``
     Name of the output file. This is the only required parameter 
     in the parameter section. If not given, the simulation will 
     terminate.
 
 
 
-**molecule block:** *<block name> = molecules*
 
+**variable block:** ``<block name>`` = *variables*
+
   This block consist of a list of pairs of the form ::
 
-     <String> <Integer>
+     String = <variable expression>
 
+  where ``String`` is the variable name, and ``<variable expression>``
+  is either a ``Double`` or an ``expression statement`` as defined above.
+  Variables can be used in any other ``expression statement`` in the
+  SGL file including reaction rate definitions. Please make sure to
+  not define a variable in terms of itself to avoid infinite recursion.
+
+
+**molecule block:** ``<block name>`` = *molecules*
+
+  This block consist of a list of pairs of the form ::
+
+     String = Integer
+
   giving the name of each molecule and the number of molecules
   present initially. For example, the following molecule definition 
   block defines molecules ``A`` and ``B`` with initial numbers of 
   100 and 200, respectively ::
 
     def molecules
-      A 100
-      B 200
+      A = 100
+      B = 200
     end
 
+  **NOTE**: Please do not use any of the predefined mathematical
+  functions or defined variables (including ``TIME``) as 
+  molecule names since this will lead to undefined behavior.
 
-**reaction block**: *<block name> = reactions*
 
+
+
+**reaction block**: ``<block name>`` = *reactions*
+
   This block describes the reactions between molecules defined in 
   the molecule block. Reactions are specified via ::
 
-     reactants -> product { rate expression }
+     <reactants> -> <products>  | <rate expression> |
 
-  Here, ``reactants`` and ``products`` are of the form ::
+  Here, ``<reactants>`` and ``<products>`` are of the form ::
 
-     <Integer> <String> + <Integer> <String> + .....
+     Integer String + Integer String + .....
 
-  In this expression, ``<String>`` is the reactant or product name 
-  as defined in the molecule block and ``<Integer>`` an optional 
-  integer specifying the stoichiometry. If ``<Integer>`` is not 
+  In this expression, ``String`` is a molecule name 
+  as defined in the molecule block and ``Integer`` an optional 
+  integer specifying the stoichiometry. If ``Integer`` is not 
   explicitly given, it is assumed to be 1.
 
-  The reaction rate can either be a fixed value of type ``<Double>`` 
-  or else an mathematical expression involving ``<Double>``, 
-  molecule names, and the current simulation time. Hence, **simgi** 
-  rate expressions can be arbitrary complex functions of the 
-  instantaneous simulation time and the instantaneous numbers of any
-  molecule in the model. The parser will interpret any string in the 
-  rate expression as a molecule name in a case sensitive fashion, 
-  a mathematical operator or function (see [5]_ for supported 
-  functions), or the special variable TIME which refers to the 
-  current simulation time. Hence, do **not** use any of the 
-  mathematical keywords as a molecule name; this leads to undefined
-  behavior.
+  The ``<rate expression>`` can either be a fixed value of type 
+  ``Double`` or an ``expression statement`` as defined above.
   
-  Here is an example reaction block for the two molecules ``A`` and 
+  Below is an example reaction block for the two molecules ``A`` and 
   ``B`` defined above::
 
     define reactions
-      2A + B -> A  { 10.0e-5 }
-      B      -> A  { 2.0e-5 * A * exp(-0.5*TIME) }
+      2A + B -> A  | 10.0e-5 |
+      B      -> A  | {2.0e-5 * A * exp(-0.5*TIME)} |
     end
    
   In the first reaction, 2 ``A`` molecules react with one ``B`` to 
-  yield another ``A`` at a rate of 10.0e-5 1/(Mol s). The second 
+  yield another ``A`` at a rate of 10.0e-5. The second 
   reaction describes a decay of ``B`` back to ``A`` at a rate 
   that is computed based on the instantaneous number of ``A`` 
   molecules present and which decays exponentially with simulation
   time.
 
-  Internally, rate expressions are converted into a compute stack
-  in RPN format which is evaluated at run-time. Even though this
-  procedure is fairly efficient, there is some numerical overhead
-  incurred at each iteration and the use of complicated rate 
-  expressions should therefore be avoided if possible.
 
+  
+**event block**: ``<block name>`` = *events*
 
+  An event block allows one to specify events which will occur during 
+  the simulation. Each event consists of a ``<trigger expression>`` and 
+  an associated set of ``<action expressions>``. 
+  Events are specified via ::
+
+     { <trigger expression> } => { <action expression> }
+
+  Here, ``trigger expression`` is of the form ::
+
+     <trigger primitive> [ <boolean operator> <trigger primitive>]
+
+  with ``<trigger primitive>`` defined by ::
+
+     <expression statement> relational operator <expression statement>
+
+  Each ``<trigger primitive>`` contains two ``expression statements``
+  as defined above and a ``relational operator`` which can be
+  any of ``>=``, ``<=``, ``==``, ``>``, and ``<``. Hence, each
+  ``<trigger primitive>`` evaluates to either ``true`` or ``false``.
+
+  Several ``<trigger primitives>`` can be chained together via the 
+  ``<boolean operators>`` ``&&`` and ``||`` to yield a final boolean
+  value of ``true`` or ``false``.
+
+  If the ``<trigger expression>`` evaluates to true during an
+  iteration, the associated ``<action expressions>`` is executed 
+  during the same timestep.
+
+  ``<action expression>`` consists of a semi-colon separated list of  
+  assignments ::
+
+    String = <assignment expression> [; String = <assignment expression>]
+
+ 
+  where ``String`` is a molecule or variable name and 
+  ``<expression>`` either a ``Double`` or an ``expression statement``.
+
+  **NOTE**: Since molecule counts are integer values assignments
+  to molecule counts in ``<action expression>`` will be converted
+  to an integer value via ``floor``.
+
+
+**output block**: ``<block name>`` = *output*
+
+  This block consists of a simple list of variable and molecule
+  names that will be streamed to the output file in the same order::
+
+    [ name1, name2, name3, .... ]
+
+
+
 Example Input Files
 -------------------
 
@@ -234,15 +317,8 @@
 <haskelladdict at users dot sourceforge dot net>. 
 
 
-References
-----------
-
 .. [1] Daniel T. Gillespie (1977). "Exact Stochastic Simulation of Coupled Chemical Reactions". The Journal of Physical Chemistry 81 (25): 2340-2361
 
-.. [2] http://hackage.haskell.org/packages/archive/random/1.0.0.1/doc/html/System-Random#globalrng.html
+.. [2] An exception to this rule are line comments starting with ``--`` which ingnore everything until the next newline.
 
-.. [3] Since **simgi** currently is an alpha version there may be fairly drastic changes to the SDL specs in future releases until the first beta release.
 
-.. [4] An exception to this rule are line comments starting with ``--`` which ingnore everything until the next newline.
-
-.. [5] Rate expressions can contain any arithmetic expression involving the standard operators "+", "-", "*", "/", "^" (exponentiation), and the mathematical functions ``sqrt, exp, log, log2, log10, sin, cos, tan, asin, acos, atan, sinh, cosh, tanh, asinh, acosh, atanh, acosh, atanh, erf, erfc, abs``.
diff --git a/simgi.cabal b/simgi.cabal
--- a/simgi.cabal
+++ b/simgi.cabal
@@ -1,5 +1,5 @@
 Name:          simgi
-Version:       0.1.1
+Version:       0.2
 License:       GPL
 license-file:  COPYING
 copyright:     (C) 2009 Markus Dittrich
@@ -9,16 +9,16 @@
                molecular systems using Gillespie's method.
 Author:        <haskelladdict@users.sourceforge.net>
 Maintainer:    <haskelladdict@users.sourceforge.net>
-stability:     alpha
+stability:     beta 
 build-type:    Simple
 Homepage:      http://simgi.sourceforge.net/
 cabal-version: >= 1.6
 extra-source-files: README
 
 Executable simgi
-  Build-Depends:  base, containers >= 0.1.0.0,
+  Build-Depends:  base >= 2 && <= 4, containers >= 0.1.0.0,
                   parsec == 2.1.*, mtl >= 1.1.0.2, haskell98,
-                  random >= 1.0.0.1
+                  random >= 1.0.0.1, mersenne-random-pure64 >= 0.2
   ghc-options:    -O2
   Main-Is:        simgi.hs
   hs-source-dirs: src
diff --git a/src/CommandLine.hs b/src/CommandLine.hs
--- a/src/CommandLine.hs
+++ b/src/CommandLine.hs
@@ -20,66 +20,101 @@
 
 -- | main gsim driver
 module CommandLine ( process_commandline 
-                   , Options(..)
+                   , SimgiOpt(..)
                    ) where
 
+
 -- imports
+import Data.Word
 import Prelude
 import System
 import System.Console.GetOpt
 
 -- local imports
+import GenericModel
 import Messages
 
 
+
 -- | main driver for command line processing
-process_commandline :: [String] -> IO ((CmdlRequest,String), [String])
-process_commandline args = 
+process_commandline :: ModelState -> [String] 
+                    -> IO (ModelState, [String])
+process_commandline state args = 
 
   let 
     (actions, nonOpts, _) = getOpt RequireOrder options args
   in
-    foldl (>>=) ( return defaultOptions ) actions >>= \opts ->
+    foldl (>>=) ( return [] ) actions >>= \opts ->
 
     let 
-      Options { cmdlRequest = request 
-              , cmdlString   = pattern
-              } = opts
+      newState = process_options opts state 
     in
-      return ((request,pattern),nonOpts) 
+      return (newState,nonOpts) 
 
 
 
--- | possible options for commandline
-data CmdlRequest = None | Help
+-- | process all user provided commandline options and adjust
+-- the model state accordingly
+process_options :: [SimgiOpt] -> ModelState -> ModelState
+process_options [] state     = state
+process_options ( (SimgiOpt { cmdlRequest = req
+                            , cmdlString  = val }):xs) state =  
 
+  case req of 
+    Seed -> process_options xs ( state { seed = parse_seed val } )
+    _    -> process_options xs state -- ignore unknown requests
 
--- | data structure for keeping track of 
--- selected command line options
-data Options = Options {
+  
+  where
+    parse_seed aSeed = floor (read aSeed :: Double) :: Word64 
+
+
+
+-- | data structure describing all commandline options we know of
+data CmdlRequest = Seed 
+
+
+
+-- | data structure for keeping track of a specific user provided
+-- command line switch 
+data SimgiOpt = SimgiOpt {
   cmdlRequest :: CmdlRequest,
   cmdlString  :: String
 }
 
 
--- | default selections
-defaultOptions :: Options
-defaultOptions = Options {
-  cmdlRequest = None,
-  cmdlString  = ""
-}
 
-
 -- | available command line flags
-options :: [OptDescr (Options -> IO Options)]
+options :: [OptDescr ([SimgiOpt] -> IO [SimgiOpt])]
 options = [
   Option ['v'] ["version-info"] (NoArg version_info) 
-         "show version information"]
+         "show version information",
+  Option ['h'] ["help"] (NoArg help_msg) "show help message",
+  Option ['s'] ["seed"] (ReqArg seed_value "SEED") "seed value"
+ ]
 
 
+
 -- | extractor function for version info
-version_info :: Options -> IO Options
+version_info :: [SimgiOpt] -> IO [SimgiOpt]
 version_info _ =
   do
     show_version
     exitWith ExitSuccess
+
+
+
+-- | extractor function for help message
+help_msg :: [SimgiOpt] -> IO [SimgiOpt]
+help_msg _ =
+  do
+    usage
+    exitWith ExitSuccess
+
+
+
+-- | extract the seed value 
+seed_value :: String -> [SimgiOpt] -> IO [SimgiOpt]
+seed_value arg opt = 
+  return ( (SimgiOpt { cmdlString = arg, cmdlRequest = Seed }):opt)
+
diff --git a/src/Engine.hs b/src/Engine.hs
--- a/src/Engine.hs
+++ b/src/Engine.hs
@@ -19,25 +19,30 @@
 --------------------------------------------------------------------}
 
 -- | the main compute Engine
-module Engine ( create_initial_output
+module Engine ( compute_trigger
+              , create_initial_output
               , create_initial_state
+              , execute_actions
               , gillespie_driver
               , module GenericModel
               ) where
 
+
 -- imports
 import Control.Monad.State
 import qualified Data.Map as M
 import Prelude
 import Text.Printf
+import System.Random(randomR)
+import qualified System.Random.Mersenne.Pure64 as MT
 import System.IO
 
+
 -- local imports
 import ExtraFunctions
 import GenericModel
 import RpnCalc
 
--- import Debug.Trace
 
 -- | main simulation driver
 -- the simulator either stops when
@@ -46,59 +51,82 @@
 --    zero t_max is treated as being infinity 
 gillespie_driver :: Handle -> Double -> Integer -> ModelState -> IO ()
 gillespie_driver handle simTime dmpIter state =  
-  let (output, outState)  = runState run_gillespie $ state 
-      (curTime, newState) = update_state dmpIter outState
+  let 
+    (output, outState)  = runState run_gillespie $ state 
+    (curTime, newState) = update_state dmpIter outState
+    reversedOutput      = reverse output
   in
-    (write_output handle . reverse $ output)
-    >> if curTime >= simTime
-         then return ()
-         else gillespie_driver handle simTime dmpIter newState
+    -- write output to console and the output file
+    (write_info $ head reversedOutput)
+    >> (write_data handle reversedOutput)
 
+    -- next iteration if we're not at the end
+    >> if curTime < simTime
+         then gillespie_driver handle simTime dmpIter newState
+         else return ()
 
 
+
 -- | updates the state for the next iteration
 update_state :: Integer -> ModelState -> (Double,ModelState)    
-update_state dataDumpIter state@(ModelState { currentTime = t 
-                                            , maxIter     = it
-                                            }) =
-  (t, state { maxIter = it + dataDumpIter, outputList = [] })
+update_state dataDumpIter 
+             state@(ModelState { currentTime      = t 
+                               , outputBufferSize = it
+                               }) 
+  = (t, state { outputBufferSize = it + dataDumpIter, outputCache = [] })
 
 
+
 -- | actual compute loop
 run_gillespie :: GillespieState [Output]
 run_gillespie = get
 
-  >>= \inState@(ModelState { molCount    = in_mols
-                           , reactions   = in_reacts
-                           , randNums    = (r1:r2:randRest)
-                           , currentTime = t
-                           , currentIter = it
-                           , maxTime     = t_max
-                           , maxIter     = it_max
-                           , outputFreq  = freq
-                           , outputList  = output
+  >>= \inState@(ModelState { molCount         = in_mols
+                           , reactions        = in_reacts
+                           , randGen          = rGen
+                           , events           = molEvents
+                           , currentTime      = t
+                           , currentIter      = it
+                           , maxTime          = t_max
+                           , outputBufferSize = it_max
+                           , outputFreq       = freq
+                           , outputRequest    = outputVars
+                           , outputCache      = output
+                           , variables        = theVars
                            }) ->
 
+
     -- compute and update the next state
-    let out_rates  = compute_rates in_reacts in_mols t []
-        a_0        = sum out_rates
-        tau        = (-1.0/a_0) * log(r1)
-        t_new      = t+tau
-        mu         = get_mu (a_0*r2) out_rates
-        out_mols   = adjust_molcount in_mols in_reacts mu
-        new_output = generate_output freq it t_new out_mols output
-        newState   = inState { molCount    = out_mols
-                             , rates       = out_rates
-                             , randNums    = randRest
-                             , currentTime = t_new
-                             , currentIter = it+1
-                             , outputList  = new_output
-                             }
+    let 
+      -- generate two random numbers
+      (r1,rGen1) = randomR (0.0 :: Double, 1.0) rGen
+      (r2,rGen2) = randomR (0.0 :: Double, 1.0) rGen1
+
+      -- update state
+      symbols    = SymbolTable in_mols theVars
+      out_rates  = compute_rates symbols in_reacts t []
+      a_0        = sum out_rates
+      tau        = (-1.0/a_0) * log(r1)
+      t_new      = t+tau
+      mu         = get_mu (a_0*r2) out_rates
+      out_mols   = update_molcount in_mols in_reacts mu
+      newSymbols = (symbols { molSymbols = out_mols })
+      evt_syms   = handle_events molEvents newSymbols t_new
+      new_output = generate_output freq it t_new newSymbols outputVars output
+      newState   = inState { molCount    = (molSymbols evt_syms)
+                           , rates       = out_rates
+                           , randGen     = rGen2
+                           , currentTime = t_new
+                           , currentIter = it+1
+                           , outputCache  = new_output
+                           , variables   = (varSymbols evt_syms)
+                           }
     in
 
     -- this prevents simulation from getting stuck
     -- FIXME: We need to come up with mechanism to propagate
-    -- error message corresponding to cases such as this one to the user!
+    -- error message corresponding to cases such as this one 
+    -- to the user!
     if (is_equal tau 0.0)
       then let finalState = newState { currentTime = t_max } in
            put finalState >> return output
@@ -108,103 +136,232 @@
           else put newState >> run_gillespie
 
 
+
+-- | handle all user defined events and return the adjusted
+-- number of molecules
+-- WARNING: We should probably check the Event Stack before we use
+-- it to compute stuff; at least make sure molecule exist
+handle_events :: [Event] -> SymbolTable -> Double -> SymbolTable
+handle_events [] symbols     _ = symbols
+handle_events (x:xs) symbols t = 
+  let
+    newSymbols = handle_single_event x symbols t
+  in
+    handle_events xs newSymbols t
+
+
+
+-- | handle a single user event
+handle_single_event :: Event -> SymbolTable -> Double -> SymbolTable
+handle_single_event evt symbols t =
+  let
+    triggers     = evtTrigger evt
+    actions      = evtActions evt
+    triggerVal   = compute_trigger symbols t triggers
+  in 
+    if triggerVal
+      then execute_actions actions symbols t 
+      else symbols
+
+
+
+-- | compute the value of a trigger
+compute_trigger :: SymbolTable -> Double 
+                -> ([EventTriggerPrimitive],[EventTriggerCombinator]) -> Bool
+compute_trigger _       _ ([],_) = False -- this is should never happen
+compute_trigger symbols t ((x:xs),combs) = compute_trigger_h (eval_trigger x) xs combs
+
+  where
+    compute_trigger_h acc []     _      = acc
+    compute_trigger_h acc _      []     = acc
+    compute_trigger_h acc (y:ys) (c:cs) = 
+
+      case c of
+        AndCombinator -> compute_trigger_h (acc && (eval_trigger y)) ys cs
+        OrCombinator  -> compute_trigger_h (acc || (eval_trigger y)) ys cs
+
+                                          
+    eval_trigger e = (trigRelation e) (leftTrigger e) (rightTrigger e)
+    leftTrigger    = rpn_compute symbols t . trigLeftExpr 
+    rightTrigger   = rpn_compute symbols t . trigRightExpr
+
+
+
+-- | handle all actions associated with a user event
+execute_actions :: [EventAction] -> SymbolTable -> Double 
+                -> SymbolTable
+execute_actions [] symbols _     = symbols
+execute_actions (x:xs) symbols t =
+  let
+    newSymbol = execute_single_action x symbols t
+  in
+    execute_actions xs newSymbol t
+
+
+
+-- | handle a single event triggered action
+execute_single_action :: EventAction -> SymbolTable -> Double
+                      -> SymbolTable
+execute_single_action eventAction symbols t =
+  let
+    aName  = evtName eventAction
+    action = evtAct eventAction
+  in
+    case action of
+      Constant c   -> adjust_count aName c 
+
+      Function rpn -> let
+                        newCount = rpn_compute symbols t rpn
+                      in
+                        adjust_count aName newCount 
+  
+      where
+        -- NOTE: presently, converting double -> int is done
+        -- via floor. Is this a good policy (once documented
+        -- properly)?
+        to_int :: Double -> Int
+        to_int = floor 
+
+        -- adjust either a molecule count or the value of a
+        -- variable
+        adjust_count key val = case M.member key (molSymbols symbols) of
+
+          True  -> symbols { molSymbols = 
+                             M.insert key (to_int val) (molSymbols symbols) }
+          False -> symbols { varSymbols = 
+                             M.insert key (Constant val) (varSymbols symbols) }
+
+
+
 -- | generate a new Output data structure based on the current
 -- molecule counts
-generate_output :: Integer -> Integer -> Double -> MoleculeMap 
-                -> [Output] -> [Output]
-generate_output afreq it t amol outlist  
+generate_output :: Integer -> Integer -> Double -> SymbolTable
+                -> [String] -> [Output] -> [Output]
+generate_output afreq it t symTable outVars outlist  
   | mod it afreq /= 0  = outlist
-  | otherwise        = new_out:outlist
+  | otherwise          = new_out:outlist
 
     where
-      new_out = Output { iteration = it
-                       , time      = t
-                       , mols      = amol
+      currentOutputList = grab_output_data outVars t symTable
+
+      new_out = Output { iteration  = it
+                       , time       = t
+                       , outputData = currentOutputList
                        }
 
 
 
+-- | given a list of variable or molecule names, goes through the
+-- symbol table, grabs the current values associated with the variables,
+-- and returns them as a list
+grab_output_data :: [String] -> Double -> SymbolTable -> [Double]
+grab_output_data vars aTime symbols = 
+  foldr (\x acc -> (get_val_from_symbolTable x aTime symbols):acc) [] vars
+                                
+
+
 -- | depending on which reaction happened adjust the number of 
 -- molecules in the system
-adjust_molcount :: MoleculeMap -> [Reaction] -> Int -> MoleculeMap
-adjust_molcount theMap rs mID =
+update_molcount :: MoleculeMap -> [Reaction] -> Int -> MoleculeMap
+update_molcount theMap rs mID =
 
-  let (Reaction { react = react_in }) = rs !! mID
+  let 
+    (Reaction { reaction = react_in }) = rs !! mID
   in
     adjustMap react_in theMap 
 
   where
     adjustMap :: [(String,Int)] -> MoleculeMap -> MoleculeMap
     adjustMap [] m = m
-    adjustMap ((k,a):rands) m = let val   = (M.!) m k
+    adjustMap ((k,a):changes) m = let 
+                                    val   = (M.!) m k
                                     m_new = M.insert k (a+val) m
-                                in
-                                  adjustMap rands m_new
+                                  in
+                                    adjustMap changes m_new
 
 
+
 -- | pick the \mu value for the randomly selected next reaction 
 -- reaction to happen
 get_mu :: Double -> [Double] -> Int
 get_mu val = length . takeWhile ( <val ) . scanl1 (+) 
 
-                 
--- | compute the current value for the reaction probabilities based on 
--- the number of molecules and reaction rates
-compute_rates :: [Reaction] -> MoleculeMap -> Double 
+
+
+-- | compute the current value for the reaction probabilities based 
+-- on the number of molecules and reaction rates
+compute_rates :: SymbolTable -> [Reaction] -> Double 
               -> RateList -> RateList
-compute_rates [] _ _ rts = reverse rts
-compute_rates ((Reaction {rate = c_in, aList = a_in }):rs) 
-  theMap theTime rts = 
+compute_rates _ [] _ rts = reverse rts
+compute_rates symbols ((Reaction {rate = c_in, actors = a_in }):rs) 
+  theTime rts = 
   
   case c_in of
-    (Constant aRate)    -> compute_rates rs theMap theTime
+    (Constant aRate)    -> compute_rates symbols rs theTime
        ((a_new aRate): rts)
-    (Function rateFunc) -> compute_rates rs theMap theTime
-       ((a_new . (rpn_compute theMap theTime) $ rateFunc):rts)
+    (Function rateFunc) -> compute_rates symbols rs theTime
+       ((a_new . (rpn_compute symbols theTime) $ rateFunc):rts)
  
   where
     mult  = product $ map (\(a,f) -> f . fromIntegral $ 
-            (M.!) theMap a) a_in 
+            (M.!) (molSymbols symbols) a) a_in 
     a_new = (*) mult 
+    
 
 
+
 -- | initialize the output data structure
 create_initial_output :: ModelState -> Output
-create_initial_output (ModelState { molCount = initialMols }) = 
+create_initial_output (ModelState { molCount      = initialMols 
+                                  , variables     = initialVars
+                                  , outputRequest = outVars
+                                  }) =
   
   Output { iteration = 1
          , time      = 0.0
-         , mols      = initialMols
+         , outputData = initialOutput
          }
 
 
+  where
+    symbols = SymbolTable initialMols initialVars
+    initialOutput = grab_output_data outVars 0.0 symbols
+
+
+
 -- | set up the initial state
-create_initial_state:: ModelState -> [Double] -> Output -> ModelState
-create_initial_state state rand output = 
+create_initial_state:: ModelState -> Output -> ModelState
+create_initial_state state@(ModelState { seed = theSeed}) out = 
 
   state { rates       = defaultRateList
-        , randNums    = rand
+        , randGen     = MT.pureMT theSeed
         , currentTime = 0.0
         , currentIter = 1
-        , outputList  = [output]
+        , outputCache = [out]
         }
 
 
--- | basic routine writing the simulation output to stdout
-write_output :: Handle -> [Output] -> IO ()
-write_output _ [] = return ()
-write_output handle ((Output {iteration = it, time = t, mols = m}):xs) = 
-  let header = (printf "%-10d %18.15g" it t) :: String
-      counts = create_count_string m
+
+-- | routine for writing basic accounting info to stdout
+write_info :: Output -> IO ()
+write_info (Output {iteration = it, time = t}) = 
+    putStrLn $ printf "iteration: %-10d  --> time: %6.5g s" it t 
+
+
+
+-- | basic routine writing the simulation output to the 
+-- file handle corresponding to the output file
+write_data :: Handle -> [Output] -> IO ()
+write_data _ [] = return ()
+write_data handle ((Output {iteration = it, time = t, outputData = out}):xs) = 
+
+  let 
+    header = (printf "%-10d %18.15g  " it t) :: String
+    counts = create_count_string out
   in
-    -- write molecule data to output file
     hPutStrLn handle (header ++ counts)
-
-    -- write current iteration to stdout
-    >> (putStrLn $ "iteration | time   --->   " ++ header)
-    >> write_output handle xs
+    >> write_data handle xs
 
   where
-    create_count_string :: MoleculeMap -> String
-    create_count_string = foldr (\x a -> (printf "%10d " x) ++ a) 
-                          "" . M.elems 
+    create_count_string :: [Double] -> String
+    create_count_string = foldr (\x a -> (printf "%18.15f  " x) ++ a) ""
diff --git a/src/ExtraFunctions.hs b/src/ExtraFunctions.hs
--- a/src/ExtraFunctions.hs
+++ b/src/ExtraFunctions.hs
@@ -26,9 +26,9 @@
                       , fact 
                       , is_equal
                       , is_equal_with
+                      , maybe_to_int
+                      , maybe_to_positive_int
                       , real_exp 
-                      , to_int
-                      , to_positive_int
                       ) where
 
 
@@ -43,20 +43,25 @@
 avogadroNum = 6.0221415e23
 
 
+
 -- | use glibc DBL_EPSILON
 dbl_epsilon :: Double
 dbl_epsilon = 2.2204460492503131e-16
 
+
+
 -- | comparison function for doubles via dbl_epsion
 is_equal :: Double -> Double -> Bool
 is_equal x y = abs(x-y) <= abs(x) * dbl_epsilon
 
 
+
 -- | comparison function for doubles via threshold
 is_equal_with :: Double -> Double -> Double -> Bool
 is_equal_with x y th = abs(x-y) <= abs(x) * th
 
 
+
 -- | function checking if a Double can be interpreted as a non
 -- negative Integer. We need this since all parsing of numbers 
 -- is done with Doubles but some functions only work for 
@@ -65,22 +70,24 @@
 -- Integer via floor and the compare if the numbers are identical.
 -- If yes, the number seems to be an Integer and we return it,
 -- otherwise Nothing
-to_positive_int :: Double -> Maybe Integer
-to_positive_int x = 
+maybe_to_positive_int :: Double -> Maybe Integer
+maybe_to_positive_int x = 
   case (is_equal (fromInteger . floor $ x) x) && (x > 0.0) of
     True  -> Just $ floor x
     False -> Nothing
 
 
+
 -- | function checking if a Double can be interpreted as an
 -- Integer. See is_positive_int for more detail
-to_int :: Double -> Maybe Integer
-to_int x = 
+maybe_to_int :: Double -> Maybe Integer
+maybe_to_int x = 
   case is_equal (fromInteger . floor $ x) x of
     True  -> Just $ floor x
     False -> Nothing
 
 
+
 -- | helper function for defining real powers
 -- NOTE: We use glibc's pow function since it is more
 -- precise than implementing it ourselves via, e.g.,
@@ -92,10 +99,12 @@
 real_exp a x = realToFrac $ c_pow (realToFrac a) (realToFrac x)
 
 
+
 -- | factorial function
 fact :: Integer -> Integer
 fact 0 = 1
 fact n = n * fact (n-1)
+
 
 
 -- | error function erf(x)
diff --git a/src/GenericModel.hs b/src/GenericModel.hs
--- a/src/GenericModel.hs
+++ b/src/GenericModel.hs
@@ -20,36 +20,85 @@
 
 -- | data structures needed for defining a stochastic model
 module GenericModel ( defaultRateList
+                    , Event(..)
+                    , EventAction(..)
+                    , EventTriggerCombinator(..)
+                    , EventTriggerPrimitive(..)
                     , GillespieState
                     , initialModelState
+                    , MathExpr(..)
                     , ModelState(..)
                     , MoleculeMap
                     , Output(..)
-                    , Rate(..)
+                    , Rate
                     , RateList
                     , Reaction(..)
+                    , SymbolTable(..)
+                    , VariableMap
+                    , VariableValue
                     ) where
 
 
 -- imports
 import Control.Monad.State
 import qualified Data.Map as M
+import Data.List((\\))
+import Data.Word
 import Prelude
+import qualified System.Random.Mersenne.Pure64 as MT
 
 
+
 -- local imports
 import RpnData
 
 
+--import Debug.Trace 
+
+
 -- | A MoleculeMap keeps track of the current number of molecules
 type MoleculeMap = M.Map String Int
 
 
--- | data type for reaction rates which can either be a Double
--- or an RpnStack describing a function to compute the rate
--- at run time
-data Rate = Constant Double | Function RpnStack
+-- | A VariableMap holds all definied variables and their 
+-- current value.
+-- NOTE: variables may change their each iteration since
+-- they may be time dependent.
+type VariableMap = M.Map String MathExpr
 
+
+-- | SymbolTable holds all names we know about such as molecule
+-- names, variable names ...
+data SymbolTable = SymbolTable { molSymbols :: MoleculeMap
+                               , varSymbols :: VariableMap
+                               }
+
+
+-- | generic data type for a mathematical expression. This could
+-- either be a constant or an expression inside an RpnStack
+data MathExpr = Constant Double | Function RpnStack
+
+
+-- | make MathExpr an instance of Eq
+-- We allow only comparison of constants with each other
+-- and RPN stacks with each other
+instance Eq MathExpr where
+ 
+  (Constant x)  == (Constant y)  =   x  == y
+  (Function f1) == (Function f2) =   f1 == f2
+  _             == _             =   False
+
+
+-- | data type for variable values which are of type MathExpr
+-- i.e. they can either be a number or a function involving, e.g.,
+-- TIME or molecule counts
+type VariableValue = Double
+
+
+-- | data type for reaction rates which are of type MathExpr
+type Rate = MathExpr
+
+
 -- | List of reactions and corresponding rates
 type RateList    = [Double]
 
@@ -57,63 +106,173 @@
 defaultRateList = [] 
 
 
--- | for each elementary reaction i we need to provide 
---   1) the reaction rate c_i or rate function 
---   2) the reaction order (first, second, ...)
---   2) aList describing which molecular species are participating
---      in a reaction (needed for computing h_mu in Gillespie's 
---      notation) and a function mapping a molecule count to the
---      proper h_mu value (needed e.g. for cases where we have
---      2X terms where h_my would be 0.5*X*(X-1).
---   3) a list of tuple (i,j) describing that the count of molecule
---      i changes by j should this reaction take place
+-- | an actor is a description of a molecular species participating
+-- in a reaction (needed for computing h_mu in Gillespie's 
+-- notation) and a function mapping a molecule count to the
+-- proper h_mu value (needed e.g. for cases where we have
+-- 2X terms where h_mu would be 0.5*X*(X-1).
+type Actor = (String, Double -> Double)
+
+
+
+-- | for each elementary reaction i we need to keep track of
+--   
+--   rate  : the reaction rate c_i or rate function 
+--   actors: a list of Actors
+--   react : a list of tuples (i,j) describing that the reaction
+--           changes the count of molecule i by j 
+--
 data Reaction = Reaction { rate       :: Rate
-                         , aList      :: [(String,Double -> Double)]
-                         , react      :: [(String,Int)]
+                         , actors     :: [Actor]
+                         , reaction   :: [(String,Int)]
                          }
 
 
+-- | make Reaction an instance of Eq so we can compare them
+-- (used in our unit tests)
+instance Eq Reaction where
+
+  x == y  =  compare_reactions x y
+
+    where
+    -- | compare two reactions and return True if they
+    -- are equal and false otherwise
+    compare_reactions :: Reaction -> Reaction -> Bool
+    compare_reactions 
+      (Reaction { rate     = rate1
+                , actors   = actors1
+                , reaction = reaction1 })
+      (Reaction { rate     = rate2
+                , actors   = actors2 
+                , reaction = reaction2 }) =
+
+      let
+        rateComp  = rate1 == rate2
+        actorComp = compare_actors actors1 actors2
+
+        -- since reaction is assembled from Data.Map in the
+        -- input parser we can't use == here
+        reactComp = ( reaction1 \\ reaction2 ) == []
+      in
+        rateComp && actorComp && reactComp
+
+
+    -- | compare two actors
+    compare_actors :: [Actor] -> [Actor] -> Bool
+    compare_actors [] []  = True
+    compare_actors xs ys  = and $ zipWith compare_actor_elem xs ys
+
+
+    -- | compare an two actor elements
+    compare_actor_elem :: Actor -> Actor -> Bool
+    compare_actor_elem act1 act2 =
+      
+      let
+        testNum  = 133.0 :: Double
+        nameComp = fst act1 == fst act2
+        funcComp = (snd act1) testNum == (snd act2) testNum
+      in
+        nameComp && funcComp
+
+
+
+
+-- | data type describing an action triggered during an event
+-- It consists of a String tracking the molecule affected
+-- as well as a mathematic expression describing the new molecule
+-- count for this molecule
+data EventAction = EventAction { evtName   :: String
+                               , evtAct    :: MathExpr
+                               }
+
+
+
+-- | data type describing an expression that triggers a 
+-- user event
+data EventTriggerPrimitive = EventTriggerPrimitive 
+  { trigLeftExpr  :: RpnStack
+  , trigRelation  :: Double -> Double -> Bool
+  , trigRightExpr :: RpnStack
+  }
+
+
+-- | combinators that can be used to combine EventTriggerPrimitives
+data EventTriggerCombinator = AndCombinator | OrCombinator
+
+
+-- | data type keeping track of possible events occuring during
+-- the simulation. Each event consist of a
+--
+--   <trigger>: list of expressions each evaluating to a bool.
+--              event is triggered if all expression evaluate to true
+--              FIXME: In the future we should support more complex 
+--                      boolen operations involving &&, ||, etc.
+-- 
+--   <action>: a list of semicolon separated expressions of the form 
+--
+--              mol/var = <numerical expression>
+--
+--             changing the value of mol/var by <numerical expression>
+--
+data Event = Event { evtTrigger :: ([EventTriggerPrimitive], [EventTriggerCombinator])
+                   , evtActions :: [EventAction]
+                   }
+
+
+
 -- | Our model state
-data ModelState = ModelState { molCount    :: MoleculeMap
-                             , rates       :: RateList
-                             , reactions   :: [Reaction]
-                             , randNums    :: [Double]
-                             , systemVol   :: Double
-                             , currentTime :: Double
-                             , currentIter :: Integer
-                             , maxTime     :: Double
-                             , maxIter     :: Integer
-                             , outputFreq  :: Integer
-                             , outputList  :: [Output]
-                             , outfileName :: String
+data ModelState = ModelState { molCount      :: MoleculeMap
+                             , rates         :: RateList
+                             , reactions     :: [Reaction]
+                             , randNums      :: [Double]
+                             , seed          :: Word64
+                             , randGen       :: MT.PureMT
+                             , events        :: [Event]
+                             , systemVol     :: Double
+                             , currentTime   :: Double
+                             , currentIter   :: Integer
+                             , maxTime       :: Double
+                             , outputBufferSize :: Integer
+                             , outputFreq    :: Integer
+                             , outputRequest :: [String]
+                             , outputCache   :: [Output]
+                             , outfileName   :: String
+                             , variables     :: VariableMap
                              }
 
 type GillespieState a = State ModelState a
 
 
+
 -- | data structure for keeping track of our output
-data Output = Output { iteration :: Integer
-                     , time      :: Double
-                     , mols      :: MoleculeMap 
+data Output = Output { iteration  :: Integer
+                     , time       :: Double
+                     , outputData :: [Double]
                      }
   deriving(Show)
 
 
+
 -- | initial model state to be partially filled by the 
 -- parser from the input deck
 initialModelState :: ModelState
-initialModelState = ModelState { molCount    = M.empty
-                               , rates       = []
-                               , reactions   = []
-                               , randNums    = []
-                               , systemVol   = 1.0
-                               , currentTime = 0.0
-                               , currentIter = 0
-                               , maxTime     = 0.0
-                               , maxIter     = 10000
-                               , outputFreq  = 1000
-                               , outputList  = []
-                               , outfileName = ""
+initialModelState = ModelState { molCount         = M.empty
+                               , rates            = []
+                               , reactions        = []
+                               , randNums         = []
+                               , events           = []
+                               , seed             = 1
+                               , randGen          = MT.pureMT 1
+                               , systemVol        = 1.0
+                               , currentTime      = 0.0
+                               , currentIter      = 0
+                               , maxTime          = 0.0
+                               , outputBufferSize = 10000
+                               , outputFreq       = 1000
+                               , outputRequest    = []
+                               , outputCache      = []
+                               , outfileName      = ""
+                               , variables        = M.empty
                                }
 
 
diff --git a/src/InputCheck.hs b/src/InputCheck.hs
--- a/src/InputCheck.hs
+++ b/src/InputCheck.hs
@@ -37,45 +37,82 @@
 -- should this ever become more extensive we should probably consider
 -- using Control.Monad.Error
 check_input :: ModelState -> Either String Bool 
-check_input (ModelState { molCount    = theMols
-                        , reactions   = theReactions
-                        , maxIter     = iterCount
-                        , outputFreq  = outFreq
-                        , outfileName = fileName
+check_input (ModelState { molCount         = theMols
+                        , reactions        = theReactions
+                        , outputBufferSize = iterCount
+                        , outputFreq       = outFreq
+                        , outfileName      = fileName
+                        , variables        = theVars
+                        , events           = theEvents
                         }) 
-  = check_molecules (defined_mols theMols) 
-      (react_mols theReactions)
+  = check_molecules (M.keys theMols) (react_mols theReactions)
     >> check_positive_outfreq outFreq
     >> check_positive_itercount iterCount
     >> check_filename fileName
-    >> check_reaction_rate_functions (defined_mols theMols) 
-         (rate_mols theReactions)
+    >> check_variable_names defined_names
+         (extract_variable_names (M.elems theVars)) "variables"
+    >> check_variable_names defined_names 
+         (extract_variable_names_from_rates theReactions) "reactions"
+    >> check_variable_names defined_names
+         (extract_variable_names_from_events theEvents) "events"
 
  where
   -- | extract all reaction participants
-  react_mols = L.nub . L.concat . map (map (fst) . react) 
+  react_mols = L.nub . L.concat . map (map (fst) . reaction) 
 
 
-  -- | extract all definied molecules
-  defined_mols = M.keys 
+  -- | extract all definied names (molecules, variables, ...)
+  defined_names = (M.keys theMols) ++ (M.keys theVars)
 
+  
 
-  -- | extract all molecules appearing in reaction rate
-  -- functions
-  rate_mols theRates = 
-    let
-      stacks    = foldr extract_rate_func [] . map rate $ theRates
-    in
-      L.nub . concat . map (foldr extract_rate_vars []) $ stacks
+-- | extract all molecules/variables appearing in reaction rate
+-- functions
+extract_variable_names_from_rates :: [Reaction] -> [String]
+extract_variable_names_from_rates = extract_variable_names . map rate
 
-      where
-        extract_rate_func (Function a) acc = a:acc
-        extract_rate_func _            acc = acc
 
-        extract_rate_vars (Variable a) acc = a:acc
-        extract_rate_vars _            acc = acc
 
+-- | extract all variable/molecule names appearing in a list of
+-- rpn stacks corresponding to some rate of variable definition
+-- expression    
+extract_variable_names :: [MathExpr] -> [String]
+extract_variable_names inputList =
+  let
+    stacks    = foldr extract_rate_func [] inputList
+  in
+    L.nub . concat . map (foldr extract_rate_vars []) $ stacks
 
+    where
+      extract_rate_func (Function a) acc = (toList a):acc
+      extract_rate_func _            acc = acc
+
+      extract_rate_vars (Variable a) acc = a:acc
+      extract_rate_vars _            acc = acc
+
+
+
+-- | extract all variable/molecule names from expressions inside
+-- events, i.e. insider triggers and actions
+extract_variable_names_from_events :: [Event] -> [String]
+extract_variable_names_from_events theEvents = 
+
+  L.nub $ allActionNames theEvents ++
+    (extract_variable_names $ triggerExps theEvents ++ actionExps theEvents)
+
+  where
+    allTriggers = concat . foldr ((:) . fst . evtTrigger) []
+    triggerExps = foldr (\x acc -> 
+                    ((Function $ trigLeftExpr x):(Function $ trigRightExpr x):acc)) 
+                    [] . allTriggers
+                    
+    allActions  = concat . foldr ((:) . evtActions) []
+    actionExps  = foldr ((:) . evtAct) [] . allActions
+
+    allActionNames = foldr ((:) . evtName) [] . allActions
+
+
+
 -- | make sure the user specified an output file name
 check_filename :: String -> Either String Bool
 check_filename name 
@@ -83,6 +120,7 @@
   | otherwise   = Right True
 
 
+
 -- | make sure all molecules in reactions are defined
 check_molecules :: [String] -> [String] -> Either String Bool
 check_molecules defMols reactMols = 
@@ -114,15 +152,15 @@
 
 -- | make sure the user defined reaction rate function reference
 -- only existing molecule names
-check_reaction_rate_functions :: [String] -> [String] 
-                              -> Either String Bool
-check_reaction_rate_functions defMols rateMols =
+check_variable_names :: [String] -> [String] -> String -> Either String Bool
+check_variable_names defMols rateMols checkType =
   let 
-    no_mol = rateMols L.\\ defMols
+    noMol = rateMols L.\\ defMols
   in
-    case null no_mol of
+    case null noMol of
       True  -> Right True
       False -> Left $
-        "Error: The following molecules defined in reaction "
-        ++ "rates do not exist: " 
-        ++ (L.concat $ L.intersperse "," no_mol)
+        "Error: The following molecules or variables defined in the "
+        ++ checkType
+        ++ " block do not exist:\n -->  " 
+        ++ (L.concat $ L.intersperse ", " noMol)
diff --git a/src/InputParser.hs b/src/InputParser.hs
--- a/src/InputParser.hs
+++ b/src/InputParser.hs
@@ -19,7 +19,10 @@
 --------------------------------------------------------------------}
 
 -- | input file parser 
-module InputParser ( input_parser ) where
+module InputParser ( input_parser
+                   , parse_events 
+                   , parse_reaction
+                   ) where
 
 -- imports
 import Control.Monad
@@ -40,34 +43,209 @@
 -- | main parser entry point
 input_parser :: CharParser ModelState ModelState
 input_parser = whiteSpace 
-               *> many ( choice [ try parse_parameter_def
+{-               *> many ( choice [ try parse_parameter_def
                                 , try parse_molecule_def
                                 , try parse_reaction_def
-                                ])
+                                , try parse_event_def
+                                ])-}
+               *> parse_parameter_def
+               *> optional (try parse_variable_def)
+               *> parse_molecule_def
+               *> parse_reaction_def
+               *> optional (try parse_event_def)
+               *> optional parse_output_def 
                *> eof
                >> getState
             <?> "main parser"
 
 
+
+-- | parser for variable definitions
+parse_variable_def :: CharParser ModelState ()
+parse_variable_def = join ( updateState <$> insert_variables <$>
+                            parse_def_block "variables" (many parse_variable) ) 
+                  <?> "variable definition block" 
+
+  where
+    insert_variables :: [(String, MathExpr)] -> ModelState -> ModelState
+    insert_variables theVars state = state { variables = M.fromList theVars }
+
+
+-- | parser for a single variable definition
+parse_variable :: CharParser ModelState (String, MathExpr)
+parse_variable = tuple_it <$> ((try parse_variable_name) )
+                 <*> (symbol "=" *> parse_variable_definition)
+              <?> "variable definition"
+  
+  where
+    tuple_it one two = (one, two)
+
+
+-- | parser for variable name
+parse_variable_name :: CharParser ModelState String
+parse_variable_name = identifier 
+                   <?> "variable name"
+
+
+-- | parse the definition for a variable
+parse_variable_definition :: CharParser ModelState MathExpr
+parse_variable_definition =  (try parse_constant_expression)
+                             <|> (braces parse_function_expression)
+                         <?> "variable value"
+
+
+-- | parser for output definitions
+parse_output_def :: CharParser ModelState ()
+parse_output_def = join ( updateState <$> insert_output_request <$>
+                         parse_def_block "output" (parse_output_list) ) 
+               <?> "event definitions" 
+
+  where
+    insert_output_request :: [String] -> ModelState -> ModelState
+    insert_output_request outDataList state = state { outputRequest = outDataList }
+
+
+-- | parse the list with variables or molecules to be punched to the 
+-- output file
+parse_output_list :: CharParser ModelState [String]
+parse_output_list = brackets (commaSep parse_variable_name)
+
+
+-- | parser for event definitions
+parse_event_def :: CharParser ModelState ()
+parse_event_def = join ( updateState <$> insert_events <$>
+                         parse_def_block "events" (many parse_events) ) 
+               <?> "event definitions" 
+
+  where
+    insert_events :: [Event] -> ModelState -> ModelState
+    insert_events newEvents state = state { events = newEvents }
+
+
+
+-- | parser for individual events
+parse_events :: CharParser ModelState Event
+parse_events = Event <$> (parse_trigger) <*> (reservedOp "=>" *> parse_actions)
+            <?> "reaction event"
+
+
+
+-- | parser for an event trigger
+parse_trigger :: CharParser ModelState 
+                 ([EventTriggerPrimitive], [EventTriggerCombinator])
+parse_trigger = braces parse_trigger_expressions
+             <?> "event trigger block"
+
+
+
+-- | parse a list of trigger expressions
+parse_trigger_expressions :: CharParser ModelState 
+                             ([EventTriggerPrimitive], [EventTriggerCombinator])
+parse_trigger_expressions = combine_it <$> parse_single_trigger_expression 
+                            <*> (many parse_boolean_trigger_expression)
+                         <?> "event trigger"
+  
+  where
+    combine_it e = foldr (\(x,y) (u,v) -> (x:u,y:v) ) ([e],[])
+
+
+
+-- | parse a single trigger expression
+parse_single_trigger_expression :: CharParser ModelState EventTriggerPrimitive
+parse_single_trigger_expression = 
+  EventTriggerPrimitive <$> parse_infix_to_rpn <*> parse_relational
+                        <*> parse_infix_to_rpn
+                               <?> "single event trigger expression"
+
+
+
+-- | parse a single trigger expression prefixed with a && or ||
+parse_boolean_trigger_expression :: CharParser ModelState 
+                                    (EventTriggerPrimitive, EventTriggerCombinator)
+parse_boolean_trigger_expression = 
+  tuple_it <$> parse_boolean_combinator <*> parse_single_trigger_expression
+                                <?> "boolean trigger expression"
+
+  where
+    tuple_it a b = (b,a)
+
+
+
+-- | parse a boolean combinator (&& or ||)
+parse_boolean_combinator :: CharParser ModelState EventTriggerCombinator
+parse_boolean_combinator = try parse_AND <|> parse_OR
+                        <?> "boolean combinator"
+ 
+
+
+-- | parse an && combinator
+parse_AND :: CharParser ModelState EventTriggerCombinator
+parse_AND = symbol "&&" *> (pure AndCombinator)
+          <?> "&&"
+  
+
+
+-- | parse an || combinator
+parse_OR :: CharParser ModelState EventTriggerCombinator
+parse_OR = symbol "||" *> (pure OrCombinator)
+        <?> "||"
+
+
+
+-- | parse a relational expression and return its associated
+-- binary function
+parse_relational :: CharParser ModelState (Double -> Double -> Bool)
+parse_relational =  try ( reservedOp ">=" >> pure (>=) )
+                <|> try ( reservedOp "<=" >> pure (<=) )
+                <|> try ( reservedOp "==" >> pure (==) )
+                <|> ( reservedOp ">" >> pure (>) )
+                <|> ( reservedOp "<" >> pure (<) )
+                <?> "relational expression"
+
+
+
+-- | parser for an event action
+parse_actions :: CharParser ModelState [EventAction]
+parse_actions = braces parse_action_expressions
+            <?> "event action block"
+
+
+
+-- | parser for a list of action expressions
+parse_action_expressions :: CharParser ModelState [EventAction]
+parse_action_expressions = 
+  parse_single_action_expression `sepEndBy` semi 
+                       <?> "event action expression"
+
+
+-- | parser for a single event action expression
+parse_single_action_expression :: CharParser ModelState EventAction
+parse_single_action_expression = EventAction <$> 
+  (molname) <*> (reservedOp "=" *> parse_function_expression)
+                              <?> "single event action expression"
+
+
+
 -- | parser for simulation parameters
 parse_parameter_def :: CharParser ModelState ()
-parse_parameter_def = between (reserved "def" *> reserved "parameters")
-                              (reserved "end")
-                              (parse_parameters `sepBy` whiteSpace) 
+parse_parameter_def = parse_def_block "parameters" (many parse_parameters) 
                       *> pure ()
                    <?> "parameter definitions"
 
 
+
 -- | parse the individual parameters
 parse_parameters :: CharParser ModelState ()
 parse_parameters = parse_time
                 <|> parse_outputFile
-                <|> parse_outputIter
+                <|> parse_outputBuffer
                 <|> parse_outputFreq
                 <|> parse_systemVol
-                <?> "time, outputIter, systemVol, outputFreq, outputFile"
+                <?> "time, outputBuffer, systemVol, outputFreq,\
+                    \outputFile"
 
 
+
 -- | parse the simulation time specs
 parse_time :: CharParser ModelState ()
 parse_time = join (updateState <$> insert_time 
@@ -78,6 +256,7 @@
     insert_time t state = state { maxTime = t }
 
 
+
 -- | parse the value of the simulated system volume
 parse_systemVol :: CharParser ModelState ()
 parse_systemVol = join (updateState <$> insert_volume
@@ -94,6 +273,7 @@
     insert_volume vol state = state { systemVol = vol }
 
 
+
 -- | parse the name of the output file 
 -- accepts paths but will NOT create any of the parents
 parse_outputFile :: CharParser ModelState ()
@@ -102,25 +282,28 @@
                           *> parse_filename ))
 
   where
-    insert_filename name state = state { outfileName = name }
+    insert_filename aName state = state { outfileName = aName }
 
 
+
 -- | parse a filename
 parse_filename :: CharParser ModelState String
 parse_filename = stringLiteral
 
 
 
+
 -- | parse the output iteration specification if present
-parse_outputIter :: CharParser ModelState ()
-parse_outputIter = join (updateState <$> insert_outputIter
-                        <$> (reserved "outputIter" *> reservedOp "="
+parse_outputBuffer :: CharParser ModelState ()
+parse_outputBuffer = join (updateState <$> insert_outputBuffer
+                     <$> (reserved "outputBuffer" *> reservedOp "="
                              *> integer ))
 
   where
-    insert_outputIter i state = state { maxIter = i }
+    insert_outputBuffer i state = state { outputBufferSize = i }
 
 
+
 -- | parse the output iteration specification if present
 parse_outputFreq :: CharParser ModelState ()
 parse_outputFreq = join (updateState <$> insert_outputFreq
@@ -131,12 +314,11 @@
     insert_outputFreq i state = state { outputFreq = i }
 
 
+
 -- | parser for molecule definitions
 parse_molecule_def :: CharParser ModelState ()
 parse_molecule_def = join ( updateState <$> insert_molecules <$> 
-  between (reserved "def" *> reserved "molecules")
-          (reserved "end")
-          (parse_molecules `sepBy` whiteSpace) )
+                            parse_def_block "molecules" (many parse_molecules)) 
                   <?> "molecule definitions"
 
   where
@@ -145,9 +327,11 @@
       state { molCount = M.fromList theMols }
 
 
+
 -- | parse a molecule name and the number of molecules of this type
 parse_molecules :: CharParser ModelState (String,Int)
-parse_molecules = make_molecule <$> (try molname) <*> integer
+parse_molecules = make_molecule <$> (try molname) <*> (symbol "=" *> integer)
+                <?> "molecule expression"
   where
     make_molecule mol aCount = (mol,fromInteger aCount)
 
@@ -156,24 +340,14 @@
 -- A molecule name can consist of letters and numbers but has to 
 -- start with a letter. The following keywords are reserved
 molname :: CharParser ModelState String
-molname = not_end ((:) <$> letter <*> many (alphaNum <?> ""))
-        <?> "molecule name" 
-
-
--- | short checker making sure we don't scan beyond the "end" statement
--- of a block 
-not_end :: CharParser ModelState String -> CharParser ModelState String
-not_end p = p >>= \name -> case name /= "end" of
-                             True  -> pure name
-                             False -> pzero
+molname = identifier
+       <?> "molecule name" 
 
 
 -- | parser for reaction definitions
 parse_reaction_def :: CharParser ModelState ()
 parse_reaction_def = join ( updateState <$> insert_reactions <$>
-                     between (reserved "def" *> reserved "reactions")
-                             (reserved "end")
-                             (parse_reaction `sepBy` whiteSpace) )
+                            parse_def_block "reactions" (many parse_reaction) ) 
                   <?> "reaction definitions"
   
   where
@@ -181,6 +355,7 @@
     insert_reactions reacts state = state { reactions = reacts }
 
 
+
 -- | parser for a single reaction specification of the type
 -- aA + bB + cC + .... -> n1P1 + n2P2 + ......   : rate :
 -- NOTE: In order to convert the reaction rates (if requested
@@ -190,7 +365,7 @@
 parse_reaction = setup_reaction 
                  <$> (parse_react_prod <* reservedOp "->") 
                  <*> parse_react_prod 
-                 <*> parse_rate
+                 <*> parse_rate_expression
                  <*> (getState 
                       >>= \(ModelState {systemVol = vol}) -> pure vol)  
   where
@@ -199,13 +374,13 @@
       let 
         action  = create_react r p
         hFactor = create_hFact r 
-        theRate = if (vol < 0.0)    -- no rate conversion for systemVol = nil
-                    then cin
+        theRate = if (vol < 0.0) -- no rate conversion for 
+                    then cin     -- systemVol = nil
                     else convert_rate cin (M.size r) vol
       in 
-        Reaction { rate       = theRate
-                 , aList      = hFactor
-                 , react      = action
+        Reaction { rate     = theRate
+                 , actors   = hFactor
+                 , reaction = action
                  }
 
 
@@ -222,21 +397,25 @@
       case order of
         1 -> theFunc
         _ -> let mult = 1.0/(avogadroNum * volume^(order-1)) in
-               Function $ stack ++ [Number mult,BinFunc (*)]
+               Function . RpnStack $ (toList stack) 
+                                     ++ [Number mult,BinFunc (*)]
 
 
 
     -- | create the list holding the molecule number changes for 
-    -- this reaction
+    -- this reactioni. If the net change in molecule number is
+    -- zero ( a nop) we remove the action.
     create_react r p = let 
                          reacts = M.map (*(-1)) r 
                        in
-                         M.assocs $ M.unionWith (+) reacts p
+                         M.assocs . snd . M.partition (==0)
+                          $ M.unionWith (+) reacts p
+      
 
 
     -- | create the list containing the h factors
-    -- WARNING/FIXME: Currently, things are ill defined if the number of
-    -- molecules for species A is below the stoichiometric reactant
+    -- WARNING/FIXME: Currently, things are ill defined if the number 
+    -- of molecules for species A is below the stoichiometric reactant
     -- coefficient; i.e. if #A = 2 then 3A -> ... does not make sense
     create_hFact :: (M.Map String Int) -> [(String, Double -> Double)]
     create_hFact     = create_hFact_h [] . M.assocs  
@@ -256,34 +435,14 @@
               generate_lambda n x   = (x-n+1) * generate_lambda (n-1) x   
 
 
--- | parse rate parses a reaction rate. This can either be a simple
--- constant of a full blown infix math expression.
--- Reaction rates must be enclosed by colons ":"
-parse_rate :: CharParser ModelState Rate
-parse_rate = (try parse_constant_rate) <|> parse_rate_function
-          <?> "rate constant or rate function"
-
-
--- | parser for a simple rate constant expression
-parse_constant_rate :: CharParser ModelState Rate
-parse_constant_rate = Constant <$> braces parse_number
-                   <?> "rate constant" 
-
-
--- | parser for a rate function
-parse_rate_function :: CharParser ModelState Rate
-parse_rate_function = Function <$> braces parse_infix_to_rpn
-                   <?> "rate function" 
-
-
 -- | parse list of reactants/products of reaction
 -- we expect to parse a stream that looks like
 -- n_1 R1 + n_2 R2 + n_3 R3 + .....
 -- If n_i is missing we assume it is 1.0
 parse_react_prod :: CharParser ModelState (M.Map String Int)
 parse_react_prod = (reserved "nil" *> pure (M.empty))
-                <|> (M.fromList <$> ((make_tuple <$> option 1 integer 
-                    <*> (try molname <* whiteSpace)) 
+                <|> (M.fromList <$> 
+                        ((make_tuple <$> option 1 integer <*> try molname) 
                          `sepBy` reservedOp "+") )
                 <?> "reactant or product list"
   
@@ -291,6 +450,7 @@
     make_tuple x y = (y, fromInteger x)
 
 
+
 -- | parse a number, can be used with 'many' and other parser
 -- combinators; integers are automatically promoted to double
 parse_number :: CharParser ModelState Double
@@ -302,6 +462,7 @@
                       Right x -> x
 
 
+
 -- | parse a positive number, can be used with 'many' and other 
 -- parser combinators; integers are automatically promoted to double
 parse_positive_number :: CharParser ModelState Double
@@ -315,3 +476,45 @@
                                 else pzero
 
   <?> "unsigned integer or double"         
+
+
+
+-- | parser for a def block structure
+parse_def_block :: String -> CharParser ModelState a 
+                -> CharParser ModelState a
+parse_def_block blockName parser = 
+
+  between (reserved "def" *> reserved blockName )
+          (reserved "end")
+          (parser)
+  <?> "parameter definitions"
+
+
+
+-- | parse a simple rate expression
+-- FIXME: We can not re-use parse_expression below 
+-- since currently the order of function/constant parsing 
+-- has to be reversed otherwise rates are always parsed
+-- as trivial Functions. 
+parse_rate_expression :: CharParser ModelState MathExpr
+parse_rate_expression = (try (lineToken parse_constant_expression))
+                        <|> (lineToken (braces parse_function_expression)) 
+                     <?> "constant or function expression"
+
+  where
+    lineToken = between (symbol "|") (symbol "|") 
+
+
+-- | parser for a simple rate constant expression
+parse_constant_expression :: CharParser ModelState MathExpr
+parse_constant_expression = Constant <$> parse_number
+                         <?> "rate constant" 
+
+
+
+-- | parser for a rate function
+parse_function_expression :: CharParser ModelState MathExpr
+parse_function_expression = Function <$> parse_infix_to_rpn
+                         <?> "rate function" 
+
+
diff --git a/src/Messages.hs b/src/Messages.hs
--- a/src/Messages.hs
+++ b/src/Messages.hs
@@ -24,22 +24,57 @@
                 , usage
                 ) where
 
+
 -- imports
 import Prelude
 
 
+-- local imports
+import GenericModel
+
+
 -- | show version info
 show_version :: IO ()
-show_version = putStrLn "This is simgi v0.1.1 (C) 2009 Markus Dittrich"
+show_version = putStrLn "This is simgi v0.2 (C) 2009 Markus Dittrich"
 
 
+
 -- | show a brief startup message
-startup_message :: IO ()
-startup_message = show_version 
-  >> putStrLn "\nstarting simulation ..... here we go\n"
+startup_message :: ModelState -> IO ()
+startup_message state = show_version
+  >> (putStrLn $ "\n-------- Simulation parameters ----------")
+  >> (putStrLn $ "rng seed              : " ++ (show simSeed))
+  >> (putStrLn $ "max time              : " ++ (show simTime) ++ " s")
+  >> (putStrLn $ "system volume         : " ++ finalSimVol)
+  >> (putStrLn $ "data output frequency : " ++ (show simFreq))
+  >> (putStrLn $ "log output frequency  : " ++ (show simLogFreq))
+  >> (putStrLn $ "output filename       : " ++ simOutFile)
+  >> (putStrLn $ "-----------------------------------------")
+  >> putStrLn "\nstarting simulation ...\n"
 
 
+  where
+    (ModelState { seed             = simSeed 
+                , maxTime          = simTime
+                , systemVol        = simVol
+                , outputBufferSize = simLogFreq
+                , outputFreq       = simFreq
+                , outfileName      = simOutFile 
+                }) = state
 
+    -- for nil volumes we want to display nil not -1.0
+    finalSimVol = if (simVol < 0) 
+                    then "nil"
+                    else (show simVol) ++ " m^3"
+
+
 -- | provide brief usage info
 usage :: IO ()
-usage = putStrLn "Usage: simgi <input file>"
+usage = putStrLn "Usage: simgi [options] <input file>\n\n\
+        \Currently supported options are:\n\n\
+        \\t -s --seed <seed value> \n\
+        \\t       specify value for the rng starting seed.\n\n\
+        \\t -v --version-info \n\
+        \\t       print version info.\n\n\
+        \\t -h --help \n\
+        \\t       print this help message."
diff --git a/src/RpnCalc.hs b/src/RpnCalc.hs
--- a/src/RpnCalc.hs
+++ b/src/RpnCalc.hs
@@ -21,7 +21,9 @@
 -- | RpnCalc defines the data structures and a calculator engine
 -- for computing mathematical expressions that have been parsed
 -- into reverse polish notations
-module RpnCalc ( rpn_compute ) where
+module RpnCalc ( rpn_compute
+               , get_val_from_symbolTable 
+               ) where
 
 
 -- imports 
@@ -32,14 +34,16 @@
 import GenericModel
 import RpnData
 
+-- import Debug.Trace
 
 -- | computes an expressions based on an rpn stack
+-- molecule names are looked up in a MoleculeMap
 -- NOTE: This function expects the RPNstack to be sanitized
 -- with respect to the variables, i.e., all variables in
--- the stack are assumed to exist in the MoleculeMap
-rpn_compute :: MoleculeMap -> Double -> RpnStack -> Double
-rpn_compute _      _    [(Number x)] = x
-rpn_compute molMap theTime xs           = num 
+-- the stack are assumed to exist in the VariableMap
+rpn_compute :: SymbolTable -> Double -> RpnStack -> Double
+rpn_compute _      _   (RpnStack [(Number x)]) = x
+rpn_compute symbols theTime (RpnStack xs)       = num 
 
   where
     (Number num) = head . foldl evaluate [] $ xs
@@ -56,9 +60,18 @@
     evaluate ys (Time) = (Number theTime):ys
 
     -- extract molecule variable
-    evaluate ys (Variable x) = (Number $ replace_var x):ys
-      where
-        replace_var :: String -> Double
-        replace_var = fromIntegral . (M.!) molMap
+    evaluate ys (Variable x) = (Number $ get_val_from_symbolTable x theTime symbols):ys 
 
     evaluate ys item = item:ys
+
+
+-- | retrieve the value of a given symbol (either variable or molecule count) from
+-- the symbol table 
+get_val_from_symbolTable :: String -> Double -> SymbolTable -> Double
+get_val_from_symbolTable var aTime symbols =
+  
+  case M.lookup var (molSymbols symbols) of
+    Just value -> fromIntegral value
+    Nothing    -> case (M.!) (varSymbols symbols) var of
+                    Constant c -> c
+                    Function s -> rpn_compute symbols aTime s
diff --git a/src/RpnData.hs b/src/RpnData.hs
--- a/src/RpnData.hs
+++ b/src/RpnData.hs
@@ -22,7 +22,7 @@
 -- for computing mathematical expressions that have been parsed
 -- into reverse polish notations
 module RpnData ( RpnItem(..)
-               , RpnStack
+               , RpnStack(..)
                ) where
 
 
@@ -38,9 +38,58 @@
                | UnaFunc (Double -> Double) 
                | BinFunc (Double -> Double -> Double)
 
-type RpnStack = [RpnItem]
 
+-- | RpnStack describes a computation stored in a stack of
+-- RpnItems
+newtype RpnStack = RpnStack { toList :: [RpnItem] }
 
+
+-- | make RpnStack an instance of Eq. 
+-- We compare two RpnStacks by computing them and replacing
+-- each occuring Variable by a default value. This should
+-- be ok in all but some pathological cases.
+instance Eq RpnStack where
+  
+  x == y  = (internal_rpn_compute x) == (internal_rpn_compute y)
+
+
+-- | computes an expressions based on an rpn stack
+-- NOTE: This function is intended to be used for
+-- the Eq instance of RpnItem only!
+-- names are replaced by defaultVar. 
+
+-- | default variable used to replace all encountered
+-- variable. Its value is arbitrary
+defaultVar :: Double
+defaultVar = 47.0
+
+internal_rpn_compute :: RpnStack -> Double
+internal_rpn_compute (RpnStack [(Number x)]) = x
+internal_rpn_compute (RpnStack xs)           = num 
+
+  where
+    (Number num) = head . foldl evaluate [] $ xs
+
+    -- evaluate unary function (sin, cos, ..)
+    evaluate ((Number x):ys) (UnaFunc f) = 
+      (Number $ f x):ys
+
+    -- evaluate binary function (*,+,..)
+    evaluate ((Number x):(Number y):ys) (BinFunc f) =
+      (Number $ f y x):ys
+
+    -- extrace current time
+    evaluate ys (Time) = (Number defaultVar):ys
+
+    -- extract molecule variable
+    evaluate ys (Variable x) = (Number defaultVar):ys
+
+    evaluate ys item = item:ys
+
+
+
+-- | a pretty lame show instance
+-- use only for debugging purposes
 instance Show RpnItem where
   show (Time)       = "TIME"
   show (Number x)   = show x
diff --git a/src/RpnParser.hs b/src/RpnParser.hs
--- a/src/RpnParser.hs
+++ b/src/RpnParser.hs
@@ -37,12 +37,13 @@
 -- | parses a mathematical infix expression and converts
 -- into a stack in rpn
 parse_infix_to_rpn :: CharParser ModelState RpnStack 
-parse_infix_to_rpn = add_term -- >>= \inp -> trace (show inp) (return inp)
+parse_infix_to_rpn = RpnStack <$> add_term
                   <?> "infix math expression"
 
 
+
 -- | parser for expressions chained via "+" 
-add_term :: CharParser ModelState RpnStack
+add_term :: CharParser ModelState [RpnItem]
 add_term = concat . insert_adds <$> 
              sub_term `sepBy` (reservedOp "+")
         <?> "addition term"
@@ -52,8 +53,9 @@
     insert_adds (x:xs) = x:foldr (\y a -> y:[BinFunc (+)]:a) [] xs
 
 
+
 -- | parser for expressions chained via "-"
-sub_term :: CharParser ModelState RpnStack
+sub_term :: CharParser ModelState [RpnItem]
 sub_term = concat . insert_subs <$> 
              div_term `sepBy` (reservedOp "-")
         <?> "subtraction term"
@@ -63,8 +65,9 @@
     insert_subs (x:xs) = x:foldr (\y a -> y:[BinFunc (-)]:a) [] xs
 
 
+
 -- | parser for expressions chained via "*" 
-div_term :: CharParser ModelState RpnStack
+div_term :: CharParser ModelState [RpnItem]
 div_term = concat . insert_divs <$> 
              mul_term `sepBy` (reservedOp "/")
         <?> "division term"
@@ -74,8 +77,9 @@
     insert_divs (x:xs) = x:foldr (\y a -> y:[BinFunc (/)]:a) [] xs
 
 
+
 -- | parser for expressions chained via "/"
-mul_term :: CharParser ModelState RpnStack
+mul_term :: CharParser ModelState [RpnItem]
 mul_term = concat . insert_muls <$> 
              exp_term `sepBy` (reservedOp "*")
         <?> "product term"
@@ -85,9 +89,10 @@
     insert_muls (x:xs) = x:foldr (\y a -> y:[BinFunc (*)]:a) [] xs
 
 
+
 -- | parser for potentiation operations "^"
-exp_term :: CharParser ModelState RpnStack
-exp_term = concat . insert_exps <$> (whiteSpace *> factor) `sepBy` (reservedOp "^")
+exp_term :: CharParser ModelState [RpnItem]
+exp_term = concat . insert_exps <$> factor `sepBy` (reservedOp "^")
         <?> "exponent"
 
   where
@@ -95,9 +100,10 @@
     insert_exps (x:xs) = x:foldr (\y a -> y:[BinFunc real_exp]:a) [] xs
 
 
+
 -- | parser for individual factors, i.e, numbers,
 -- variables or operations
-factor :: CharParser ModelState RpnStack
+factor :: CharParser ModelState [RpnItem]
 factor = try parse_single_number  -- need try due to the unary "-"
       <|> try signed_parenthesis  -- (otherwise we get stuck)
       <|> parse_functions
@@ -105,9 +111,10 @@
       <?> "token or variable"         
 
 
+
 -- | parse all operations of type (Double -> Double)
 -- we currently know about
-parse_functions :: CharParser ModelState RpnStack
+parse_functions :: CharParser ModelState [RpnItem]
 parse_functions = (msum $ extract_ops builtinFunctions)
                <?> "builtin unary function"
   where
@@ -124,10 +131,11 @@
     insert_func f xs = xs ++ [UnaFunc f]
 
 
+
 -- | parse a potentially signed expression enclosed in parenthesis.
 -- In the case of parenthesised expressions we 
 -- parse -(...) as (-1.0)*(...)
-signed_parenthesis :: CharParser ModelState RpnStack
+signed_parenthesis :: CharParser ModelState [RpnItem]
 signed_parenthesis = push_parens <$> parse_sign <*> parens add_term
                   <?> "signed parenthesis"
 
@@ -135,17 +143,19 @@
     push_parens sign xs = xs ++ [Number sign,BinFunc (*)]
 
 
+
 -- | parse a single number; integers are automatically promoted 
 -- to double
 -- NOTE: Due to the notFollowedBy this parser can not be used
 -- with 'many' and other parser combinators.
-parse_single_number :: CharParser ModelState RpnStack
-parse_single_number = push <$> (parse_number <* notFollowedBy alphaNum)
+parse_single_number :: CharParser ModelState [RpnItem]
+parse_single_number = push <$> (parse_number) 
                    <?> "signed integer or double"
   where
     push x = [Number x]
 
 
+
 -- | parse a number, can be used with 'many' and other parser
 -- combinators; integers are automatically promoted to double
 parse_number :: CharParser ModelState Double
@@ -157,17 +167,17 @@
                            Right x -> sign * x
 
 
+
 -- | parse the sign of a numerical expression
 parse_sign :: CharParser ModelState Double
-parse_sign = option 1.0 ( whiteSpace *> char '-' *> pure (-1.0) )
+parse_sign = option 1.0 ( char '-' *> pure (-1.0) )
           <?> "sign"
 
 
+
 -- | this is how valid variable names have to look like
-parse_variable :: CharParser ModelState RpnStack
-parse_variable = 
-  push <$> parse_sign 
-  <*> ((:) <$> letter <*> many (alphaNum <?> "") <* whiteSpace)
+parse_variable :: CharParser ModelState [RpnItem]
+parse_variable = push <$> parse_sign <*> identifier 
               <?> "variable"
   where
     -- in case of a unary minus we also push the necessary
diff --git a/src/TokenParser.hs b/src/TokenParser.hs
--- a/src/TokenParser.hs
+++ b/src/TokenParser.hs
@@ -22,15 +22,20 @@
 module TokenParser ( module Control.Applicative
                    , module Text.ParserCombinators.Parsec
                    , braces
+                   , brackets
                    , builtinFunctions
+                   , colon
                    , comma
+                   , commaSep
                    , charLiteral
                    , float
+                   , identifier
                    , integer
                    , parens
                    , keywords
                    , lexer
                    , naturalOrFloat
+                   , operator
                    , operators
                    , reservedOp
                    , reserved
@@ -48,7 +53,6 @@
 import Text.ParserCombinators.Parsec hiding (many,optional, (<|>)) 
 import qualified Text.ParserCombinators.Parsec.Token as PT
 import Text.ParserCombinators.Parsec.Language (haskellDef
-                                              , opLetter
                                               , reservedOpNames
                                               , reservedNames )
 
@@ -64,6 +68,7 @@
   (<*>) = ap
 
 
+
 -- |Alternative instance for MonadPlus
 instance Alternative (GenParser s a) where
   empty = mzero
@@ -98,36 +103,47 @@
                    ]
 
 
+
 -- | all other keywords that are not regular functions
 keywords :: [String]
 keywords = [ "def", "molecules", "reactions", "time", "outputIter"
-           , "nil", "outputFreq", "outputFile", "systemVol"
+           , "nil", "outputFreq", "outputFile", "systemVol", "seed"
+           , "end", "variables", "time"
            ]
 
 operators :: [String]
-operators = ["+","->","::","=","{","}"]
+operators = ["+", "->", "::", "=", "{", "}", ">=", "==", "<="
+            , "<", ">", "*", "/", "-", "&&"]
 
 
+
 -- | function generating a token parser based on a 
 -- lexical parser combined with a language record definition
 lexer :: PT.TokenParser st
 lexer  = PT.makeTokenParser 
          ( haskellDef { reservedOpNames = operators
-                      , opLetter      = oneOf "*+/^"
                       , reservedNames = keywords 
                             ++ map fst builtinFunctions
                       } )
 
 
+
 -- | token parser for parenthesis
 parens :: CharParser st a -> CharParser st a
 parens = PT.parens lexer
 
+
+
 -- | token parser for parenthesis
 braces :: CharParser st a -> CharParser st a
 braces = PT.braces lexer
 
 
+-- | token parser for brackets 
+brackets :: CharParser st a -> CharParser st a
+brackets = PT.brackets lexer
+
+
 -- | token parser for Integer
 integer :: CharParser st Integer
 integer = PT.integer lexer
@@ -158,27 +174,53 @@
 reservedOp = PT.reservedOp lexer
 
 
+
 -- | token parser for keywords
 reserved :: String -> CharParser st ()
 reserved = PT.reserved lexer
 
 
+
 -- | token parser for whitespace
 whiteSpace :: CharParser st ()
 whiteSpace = PT.whiteSpace lexer
 
 
+
+-- | token parser for colon
+colon:: CharParser st String
+colon = PT.colon lexer
+
+
+
 -- | token parser for semicolon
 semi :: CharParser st String
 semi = PT.semi lexer
 
 
+
 -- | token parser for comma
 comma :: CharParser st String
 comma = PT.comma lexer
 
 
+-- | token parser for comma separated list of items
+commaSep :: CharParser st a -> CharParser st [a]
+commaSep = PT.commaSep lexer
+
+
 -- | token parser for symbol
 symbol :: String -> CharParser st String
 symbol = PT.symbol lexer
+
+
+-- | token parser for symbol
+identifier :: CharParser st String
+identifier = PT.identifier lexer
+
+
+-- | token parser for symbol
+operator:: CharParser st String
+operator = PT.operator lexer
+
 
diff --git a/src/simgi.hs b/src/simgi.hs
--- a/src/simgi.hs
+++ b/src/simgi.hs
@@ -24,7 +24,7 @@
 -- imports
 import Prelude
 import System.IO
-import System.Random
+--import System.Random
 import System.Environment
 
 -- local imports
@@ -43,9 +43,10 @@
 main = 
 
   -- process command line arguments
-  getArgs >>= process_commandline 
-  >>= \((_,_),files) -> 
+  getArgs >>= process_commandline initialModelState
+  >>= \(state, files) -> 
 
+
   -- reject anything but a single input file
   if length files /= 1 
     then usage
@@ -56,28 +57,29 @@
       >>= \content -> 
 
         -- parse input file
-        case runParser input_parser initialModelState "" content of
+        case runParser input_parser state "" content of
           Left er           -> putStrLn (show er)
           Right parsedState -> 
             case check_input parsedState of
               Left err -> putStrLn err
-              Right _  -> newStdGen
-                >>= \gen -> 
+              Right _  -> 
 
                 -- set up simuation
                 let 
-                  rands         = randomRs (0,1) gen :: [Double] 
                   initialOutput = create_initial_output parsedState
                   initialState  = create_initial_state parsedState 
-                                  rands initialOutput
+                                   initialOutput
                   totalTime     = maxTime parsedState 
-                  dataDumpIter  = maxIter parsedState
+                  dataDumpIter  = outputBufferSize parsedState
                   outFile       = outfileName parsedState
                 in
 
                   -- open output file
                   openFile outFile WriteMode 
-                  >>= \handle -> startup_message
+                  >>= \handle -> 
+
+                  -- print initial startup info
+                  startup_message initialState
 
                   -- ready to run the simulation
                   >> gillespie_driver handle totalTime dataDumpIter 
diff --git a/test/Makefile b/test/Makefile
--- a/test/Makefile
+++ b/test/Makefile
@@ -3,13 +3,37 @@
 default: check
 
 
-.PHONY: check clean
-check:
-	make -C rpnstack_test
-	make -C reversible_test
+.PHONY: check check_event_parser check_rpnstack check_reversible \
+	check_reaction_parser check_irreversible clean 
 
 
-clean:
-	make -C rpnstack_test clean
-	make -C reversible_test clean
+check_rpnstack:
+	make -C $@
+
+
+check_reversible:
+	make -C $@
+
+check_irreversible:
+	make -C $@	
+
+check_event_parser:
+	make -C $@
+
+
+check_reaction_parser:
+	make -C $@
+
+check: check_rpnstack check_event_parser check_reaction_parser \
+	check_reversible check_irreversible 
+
+
+
+clean: 
+	make -C check_rpnstack clean
+	make -C check_reversible clean
+	make -C check_irreversible clean
+	make -C check_event_parser clean
+	make -C check_reaction_parser clean
+
 	
diff --git a/test/check_event_parser/EventParserTest.hs b/test/check_event_parser/EventParserTest.hs
new file mode 100644
--- /dev/null
+++ b/test/check_event_parser/EventParserTest.hs
@@ -0,0 +1,248 @@
+{-----------------------------------------------------------------
+ 
+  (c) 2009 Markus Dittrich 
+ 
+  This program is free software; you can redistribute it 
+  and/or modify it under the terms of the GNU General Public 
+  License Version 3 as published by the Free Software Foundation. 
+ 
+  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 Version 3 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, Inc., 59 Temple Place - Suite 330, 
+  Boston, MA 02111-1307, USA.
+
+--------------------------------------------------------------------}
+
+-- | this routine tests some aspects of the event block parsing
+-- routines
+module Main where
+
+
+-- imports 
+import Control.Monad.Writer
+import qualified Data.Map as M
+import Prelude
+import System.Exit
+
+-- local imports
+import Engine
+import GenericModel
+import PrettyPrint
+import InputParser
+import TestHelpers
+import TokenParser
+
+
+-----------------------------------------------------------------
+-- 
+-- test data specifications
+--
+-----------------------------------------------------------------
+
+-- | expected values for parser events. The first
+-- entry in the tuple is of type Bool giving the expected
+-- value of the evaluated trigger expression. The second
+-- is a Map of expected molecule counts after the actions
+-- have been evaluated
+type ExpectedOutput = (Bool, M.Map String Int)
+
+
+
+-- | set up the test data
+-- Format : (parse expression, expected result)
+type TestCase = (String, ExpectedOutput)
+
+
+
+-- | simple tests without access to local variables
+--- NOTE: for now we simply test if parsing succeeds
+simpleEventParseTests :: [TestCase]
+simpleEventParseTests = 
+  [ ("{ z == 100 } => { x = x }", 
+      (True, M.fromList [("x",1000),("y",2000),("z",100)]))
+
+  , ("{ z== 100} => { x = x+y; y = y+1 }", 
+      (True, M.fromList [("x",3000),("y",2001),("z",100)]))
+
+  , ("{z==121 } => { y=10 }", 
+      (False, M.fromList [("x",1000),("y",10),("z",100)]))
+
+  , ("{z==101} => { x = x/2.0 }", 
+      (False, M.fromList [("x",500),("y",2000),("z",100)]))
+
+  , ("{ z == 100}=>{ z= 121; x = 500; y=2000;}", 
+      (True, M.fromList [("x",500),("y",2000),("z",121)]))
+
+  , ("{ y >= 121 }=>{z=1050}", 
+      (True, M.fromList [("x",1000),("y",2000),("z",1050)]))
+
+  , ("{ z <= 100} => { x= sqrt(4)*100; y = x^2/10}", 
+      (True, M.fromList [("x",200),("y",4000),("z",100)]))
+
+  , ("{ x >= 100} => { z= z; z = 2*z; z = 3*z;}", 
+      (True, M.fromList [("x",1000),("y",2000),("z",600)]))
+
+  , ("{ z < 100} => { x= 10; y = z}", 
+      (False, M.fromList [("x",10),("y",100),("z",100)]))
+
+  , ("{ x > 100} => { y= x+z}", 
+      (True, M.fromList [("x",1000),("y",1100),("z",100)]))
+
+  , ("{ y == 2000 } => { x=1001 }", 
+      (True, M.fromList [("x",1001),("y",2000),("z",100)]))
+
+  , ("{ z <= 302 } => { x = 10; y = 10; z= 200; z=10 }", 
+      (True, M.fromList [("x",10),("y",10),("z",10)]))
+
+  , ("{ y == 100 } => { x = 10; y = 20; }", 
+      (False, M.fromList [("x",10),("y",20),("z",100)]))
+
+  , ("{ z == 100 } => { x = 10; y = 20 }", 
+      (True, M.fromList [("x",10),("y",20),("z",100)]))
+
+  , ("{ z == 100 } => { x = 10; y = 20 +      z }", 
+      (True, M.fromList [("x",10),("y",120),("z",100)]))
+
+  , ("{ x > 100 } => { x = 10*x+y; y = 20+z^2 }", 
+      (True, M.fromList [("x",12000),("y",10020),("z",100)]))
+
+  , ("{ z == 100 } => { x = 1e6*exp(-TIME); y = x; }",
+      (True, M.fromList [("x",4),("y",4),("z",100)]))
+
+  , ("{ y == 100 } => { x = 2; x = y^x; y = x*exp(-TIME) }", 
+      (False, M.fromList [("x",4000000),("y",17),("z",100)]))
+
+  , ("{ z == 100 } => { x = 10+z; y = sqrt(log((x^2))); }",
+      (True, M.fromList [("x",110),("y",3),("z",100)])) 
+  ]
+
+
+----------------------------------------------------------------
+-- Molecule maps and definitions for specific simulation times
+----------------------------------------------------------------
+
+-- | testmap containing a set of molecules and their 
+-- concentrations
+testMap_1 :: MoleculeMap
+testMap_1 = M.fromList [("x",1000),("y",2000),("z",100)]
+
+-- | a simulation time
+time_1 :: Double
+time_1 = 12.345
+
+
+
+
+----------------------------------------------------------------
+--
+--  main driver routines
+--
+----------------------------------------------------------------
+
+
+-- | main test driver
+main :: IO ()
+main = putStrLn "\n\n\nTesting Input Parser"
+
+  -- check event parser
+  >> (putStr $ color_string Cyan "\nEvent parse tests:\n")
+  >> let eventParseOut = execWriter $ event_parser_test_driver 
+           testMap_1 time_1 simpleEventParseTests
+     in
+  examine_output eventParseOut >>= \eventParseStatus ->
+
+
+  -- check event triggers and actions
+  (putStr $ color_string Cyan "\n\nEvent trigger/action tests:\n")
+  >> let eventActionOut = execWriter $ event_action_test_driver 
+           testMap_1 time_1 simpleEventParseTests
+     in
+  examine_output eventActionOut >>= \eventActionStatus ->
+
+
+  -- evaluate status and return
+  let status = eventParseStatus && eventActionStatus 
+  in
+    if status == True then
+      exitWith ExitSuccess
+    else
+      exitWith $ ExitFailure 1 
+
+
+
+-- | this driver parses the event expression and checks that
+-- there are no errors (there shouldn't be any)
+event_parser_test_driver :: MoleculeMap -> Double -> [TestCase] 
+            -> Writer [TestResult] ()
+event_parser_test_driver _   _    []     = return ()
+event_parser_test_driver mol t (x:xs) =
+
+  let expr     = fst x
+      expected = snd x
+  in
+
+    -- parse expression
+    case runParser parse_events testModelState "" expr of
+      Left er -> tell [TestResult False expr (show expected) (show er)]
+      Right _ -> 
+        tell [TestResult True expr (show expected) ("good parse")]
+        >> event_parser_test_driver mol t xs
+          
+
+
+-- | this driver parses the event expression and checks if
+-- the trigger and its actions action check out, i.e., we apply 
+-- them to our current state and see if we get the proper number 
+-- of molecules.
+-- NOTE: We do not thread the state through all tests but always
+-- start with our default, otherwise it becomes very tedious
+-- to add/remove tests.
+-- This parser should be run after the proper parsing has been
+-- verified via event_parser_test_driver
+event_action_test_driver :: MoleculeMap -> Double -> [TestCase] 
+            -> Writer [TestResult] ()
+event_action_test_driver _   _    []     = return ()
+event_action_test_driver mol t (x:xs) =
+
+  let expr            = fst x
+      expectedTrigger = fst . snd $ x
+      expectedMols    = snd . snd $ x 
+  in
+
+    -- parse expression
+    case runParser parse_events testModelState "" expr of
+      Left er     -> tell [TestResult False expr "" (show er)]
+      Right event -> 
+        
+        -- make sure the trigger expression evaluated properly
+        let
+          actions    = evtActions event
+          outMols    = execute_actions actions (SymbolTable mol M.empty) t 
+          outTrigger = check_trigger mol t expectedTrigger event 
+        in
+          case outTrigger && ((molSymbols outMols) == expectedMols) of
+
+            False -> tell [TestResult False expr 
+              (show expectedTrigger ++ " => " ++ show expectedMols)
+              (show outTrigger ++ " => " ++ show (molSymbols outMols))]
+
+            True  -> tell [TestResult True expr "" ("good parse")]
+                       >> event_action_test_driver mol t xs
+
+
+
+-- | given an Event and a MoleculeMap check that the trigger
+-- evaluates to the expected value
+check_trigger :: MoleculeMap -> Double -> Bool -> Event
+                         -> Bool
+check_trigger molMap t expected (Event { evtTrigger = trigger }) =
+  computed == expected
+  
+  where
+    computed = compute_trigger (SymbolTable molMap M.empty) t trigger
+
diff --git a/test/check_event_parser/Makefile b/test/check_event_parser/Makefile
new file mode 100644
--- /dev/null
+++ b/test/check_event_parser/Makefile
@@ -0,0 +1,15 @@
+include ../test_helpers/make_common
+
+default: check
+
+
+.PHONY: check clean
+
+check:
+	ghc $(GHC_FLAGS_DEVEL) -i$(SRC_LOCATION) -i$(HELPER_LOCATION) \
+		--make EventParserTest.hs
+	@./EventParserTest
+
+
+clean:
+	rm -f *.hi *.o EventParserTest
diff --git a/test/check_irreversible/Makefile b/test/check_irreversible/Makefile
new file mode 100644
--- /dev/null
+++ b/test/check_irreversible/Makefile
@@ -0,0 +1,19 @@
+include ../test_helpers/make_common
+
+
+default: check
+
+.PHONY: check clean
+check:
+	@echo
+	@echo "Running irreversible reaction test 1"
+	@echo "This may take a while ...."
+	@./run_test.1.sh $(SIMGI_PATH)
+	@echo "Running irreversible reaction test 2"
+	@echo "This may take a while ...."
+	@./run_test.2.sh $(SIMGI_PATH)
+	@echo
+
+
+clean:
+	rm -f *.dat
diff --git a/test/check_irreversible/irreversible.1.sgl b/test/check_irreversible/irreversible.1.sgl
new file mode 100644
--- /dev/null
+++ b/test/check_irreversible/irreversible.1.sgl
@@ -0,0 +1,22 @@
+
+def parameters
+  time       = 1.0
+  systemVol  = 1.25e-19 
+  outputBuffer = 50000 
+  outputFreq = 100
+  outputFile = "irreversible.1.dat"
+end
+
+def molecules
+  a = 1000
+  b = 1000
+  c = 0
+end
+
+def reactions
+  a + b -> c  | 1e6 |
+end 
+
+def output
+  [a, b, c]
+end
diff --git a/test/check_irreversible/irreversible.2.sgl b/test/check_irreversible/irreversible.2.sgl
new file mode 100644
--- /dev/null
+++ b/test/check_irreversible/irreversible.2.sgl
@@ -0,0 +1,22 @@
+
+def parameters
+  time       = 1.0
+  systemVol  = 1.25e-19 
+  outputBuffer = 50000 
+  outputFreq = 100
+  outputFile = "irreversible.2.dat"
+end
+
+def molecules
+  a = 1000
+  b = 1000
+  c = 0
+end
+
+def reactions
+  b -> c  | 1e6 |
+end 
+
+def output
+  [a, b, c]
+end
diff --git a/test/check_irreversible/run_test.1.sh b/test/check_irreversible/run_test.1.sh
new file mode 100644
--- /dev/null
+++ b/test/check_irreversible/run_test.1.sh
@@ -0,0 +1,60 @@
+#!/bin/bash
+#
+# short wrapper for runing a number of reversible
+# simulations and then checking if we can reproduce
+# the correct average product concentration
+#
+
+# expected values of molecule counts
+a_expect=0.000000000000000
+b_expect=0.000000000000000
+c_expect=1000.000000000000000
+
+# 
+simgi_exe="${1}"
+
+# our global status
+status=0
+
+for ((counter=0; counter <= 100; counter++)); do
+
+  # provide a little "progressbar"
+  printf "."
+
+  # run and process
+  ${simgi_exe} irreversible.1.sgl >& /dev/null 
+
+  a=$(tail -n 1 irreversible.1.dat | gawk ' { print $3 }')
+  b=$(tail -n 1 irreversible.1.dat | gawk ' { print $4 }')
+  c=$(tail -n 1 irreversible.1.dat | gawk ' { print $5 }')
+
+  if [[ ${a_expect} != ${a} || ${b_expect} != ${b} \
+        || ${c_expect} != ${c} ]];
+  then
+    status=1
+  fi
+
+  # unlink
+  rm -f irreversible.1.dat || return 1
+  
+done
+
+# check if all the tests passed
+
+# brief output
+if [[ ${status} == 0 ]]; then
+  echo
+  echo
+  echo "Congratulations - the irreversible reaction test 1 passed!"
+  echo
+  echo
+else
+  echo
+  echo
+  echo "Error - the irreversible reaction test failed. Please check!"
+  echo
+  echo
+fi
+
+# return status of deviation check
+exit ${status}
diff --git a/test/check_irreversible/run_test.2.sh b/test/check_irreversible/run_test.2.sh
new file mode 100644
--- /dev/null
+++ b/test/check_irreversible/run_test.2.sh
@@ -0,0 +1,60 @@
+#!/bin/bash
+#
+# short wrapper for runing a number of reversible
+# simulations and then checking if we can reproduce
+# the correct average product concentration
+#
+
+# expected values of molecule counts
+a_expect=1000.000000000000000
+b_expect=0.000000000000000
+c_expect=1000.000000000000000
+
+# 
+simgi_exe="${1}"
+
+# our global status
+status=0
+
+for ((counter=0; counter <= 100; counter++)); do
+
+  # provide a little "progressbar"
+  printf "."
+
+  # run and process
+  ${simgi_exe} irreversible.2.sgl >& /dev/null 
+
+  a=$(tail -n 1 irreversible.2.dat | gawk ' { print $3 }')
+  b=$(tail -n 1 irreversible.2.dat | gawk ' { print $4 }')
+  c=$(tail -n 1 irreversible.2.dat | gawk ' { print $5 }')
+
+  if [[ ${a_expect} != ${a} || ${b_expect} != ${b} \
+        || ${c_expect} != ${c} ]];
+  then
+    status=1
+  fi
+
+  # unlink
+  rm -f irreversible.2.dat || return 1
+  
+done
+
+# check if all the tests passed
+
+# brief output
+if [[ ${status} == 0 ]]; then
+  echo
+  echo
+  echo "Congratulations - the irreversible reaction test 2 passed!"
+  echo
+  echo
+else
+  echo
+  echo
+  echo "Error - the irreversible reaction test failed. Please check!"
+  echo
+  echo
+fi
+
+# return status of deviation check
+exit ${status}
diff --git a/test/check_reaction_parser/Makefile b/test/check_reaction_parser/Makefile
new file mode 100644
--- /dev/null
+++ b/test/check_reaction_parser/Makefile
@@ -0,0 +1,15 @@
+include ../test_helpers/make_common
+
+default: check
+
+
+.PHONY: check clean
+
+check:
+	ghc $(GHC_FLAGS_DEVEL) -i$(SRC_LOCATION) -i$(HELPER_LOCATION) \
+		--make ReactionParserTest.hs
+	@./ReactionParserTest
+
+
+clean:
+	rm -f *.hi *.o ReactionParserTest
diff --git a/test/check_reaction_parser/ReactionParserTest.hs b/test/check_reaction_parser/ReactionParserTest.hs
new file mode 100644
--- /dev/null
+++ b/test/check_reaction_parser/ReactionParserTest.hs
@@ -0,0 +1,200 @@
+{-----------------------------------------------------------------
+ 
+  (c) 2009 Markus Dittrich 
+ 
+  This program is free software; you can redistribute it 
+  and/or modify it under the terms of the GNU General Public 
+  License Version 3 as published by the Free Software Foundation. 
+ 
+  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 Version 3 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, Inc., 59 Temple Place - Suite 330, 
+  Boston, MA 02111-1307, USA.
+
+--------------------------------------------------------------------}
+
+-- | this routine tests some aspects of the reaction block
+-- parsing routines
+module Main where
+
+
+-- imports 
+import Control.Monad.Writer
+import qualified Data.Map as M
+import Prelude
+import System.Exit
+
+-- local imports
+import GenericModel
+import PrettyPrint
+import InputParser
+import RpnData
+import TestHelpers
+import TokenParser
+
+
+-----------------------------------------------------------------
+-- 
+-- test data specifications
+--
+-----------------------------------------------------------------
+
+-- | expected values for parser events. The first
+-- entry in the tuple is of type Bool giving the expected
+-- value of the evaluated trigger expression. The second
+-- is a Map of expected molecule counts after the actions
+-- have been evaluated
+type ExpectedOutput = Reaction
+
+
+
+-- | set up the test data
+-- Format : (parse expression, expected result)
+type TestCase = (String, ExpectedOutput)
+
+
+
+-- | simple tests without access to local variables
+--- NOTE: for now we simply test if parsing succeeds
+simpleReactParseTests :: [TestCase]
+simpleReactParseTests = 
+  [ ("x + y -> z | 1e6 |",
+    Reaction (Constant 1e6) 
+             [("y",id),("x",id)]
+             [("x",-1),("y",-1),("z",1)])
+
+  , ("x + y + z -> w | 10.3 |",
+    Reaction (Constant 10.3)
+             [("z",id),("y",id),("x",id)]
+             [("x",-1),("y",-1),("z",-1),("w",1)])
+
+  , ("x + y + z -> w | {x*y} |",
+    Reaction (Function . RpnStack $ [Variable "x", Variable "y", 
+              BinFunc (*)])
+             [("z",id),("y",id),("x",id)]
+             [("x",-1),("y",-1),("z",-1),("w",1)])
+
+  , ("2x + 2y + z -> nil | { exp(-TIME) } |",
+    Reaction (Function . RpnStack $ [Variable "TIME", Number (-1.0)
+              , BinFunc (*), UnaFunc exp])
+             [("z",id),("y",\a -> 0.5 * a * (a-1))
+              ,("x",\a -> 0.5 * a * (a-1))]
+             [("x",-2),("y",-2),("z",-1)])
+
+  , ("2x + y + z -> 2x + 20y + 5z | {x+y+z } |",
+    Reaction (Function . RpnStack $ [Variable "x", Variable "y", 
+              BinFunc (+), Variable "z", BinFunc (+)])
+             [("z",id),("y",id),("x",\a -> 0.5 * a * (a-1))]
+             [("y",19),("z",4)])
+
+  , ("2x + y + z -> 2x + 20y + 5z | { x*   y*  z } |",
+    Reaction (Function . RpnStack $ [Variable "x", Variable "y", 
+              BinFunc (*), Variable "z", BinFunc (*)])
+             [("z",id),("y",id),("x",\a -> 0.5 * a * (a-1))]
+             [("y",19),("z",4)])
+
+  , ("2   x+ y+ z-> 2  x + 20 y + 5z|{x*y   *  z } |",
+    Reaction (Function . RpnStack $ [Variable "x", Variable "y", 
+              BinFunc (*), Variable "z", BinFunc (*)])
+             [("z",id),("y",id),("x",\a -> 0.5 * a * (a-1))]
+             [("y",19),("z",4)])
+
+  , ("x + y + z -> x + y + z |{ x-x+y-y+z-z }|",
+    Reaction (Function . RpnStack $ [Number 0.0])
+             [("z",id),("y",id),("x",id)]
+             [])
+
+  , ("nil -> x + y + z |{ x-x+y-y+z-z }|",
+    Reaction (Function . RpnStack $ [Number 0.0])
+             []
+             [("x",1),("y",1),("z",1)])
+
+  , ("nil -> nil |{ sqrt(2.0) }  |",
+    Reaction (Function . RpnStack $ [Number 2.0, UnaFunc sqrt])
+             []
+             [])
+  ]
+
+
+
+----------------------------------------------------------------
+-- Molecule maps and definitions for specific simulation times
+----------------------------------------------------------------
+
+-- | testmap containing a set of molecules and their 
+-- concentrations
+testMap_1 :: MoleculeMap
+testMap_1 = M.fromList [("x",1000),("y",2000),("z",100)]
+
+-- | a simulation time
+time_1 :: Double
+time_1 = 12.345
+
+
+
+
+----------------------------------------------------------------
+--
+--  main driver routines
+--
+----------------------------------------------------------------
+
+
+-- | main test driver
+main :: IO ()
+main = putStrLn "\n\n\nTesting Reaction Parser"
+
+  -- check event parser
+  >> (putStr $ color_string Cyan "\nReaction parse tests:\n")
+  >> let reactParseOut = execWriter $ react_parser_test_driver 
+           testMap_1 time_1 simpleReactParseTests
+     in
+  examine_output reactParseOut >>= \reactParseStatus ->
+
+
+  -- evaluate status and return
+  let status = reactParseStatus 
+  in
+    if status == True then
+      exitWith ExitSuccess
+    else
+      exitWith $ ExitFailure 1 
+
+
+
+-- the trigger and its actions action check out, i.e., we apply 
+-- them to our current state and see if we get the proper number 
+-- of molecules.
+-- NOTE: We do not thread the state through all tests but always
+-- start with our default, otherwise it becomes very tedious
+-- to add/remove tests.
+-- This parser should be run after the proper parsing has been
+-- verified via event_parser_test_driver
+react_parser_test_driver :: MoleculeMap -> Double -> [TestCase] 
+            -> Writer [TestResult] ()
+react_parser_test_driver _   _    []     = return ()
+react_parser_test_driver mol t (x:xs) =
+
+  let
+    parseString = fst x
+    expected    = snd x
+  in
+
+    -- parse expression
+    case runParser parse_reaction testModelState "" parseString of
+      Left er     -> tell [TestResult False parseString "" (show er)]
+      Right react -> 
+         
+         case (react == expected) of
+
+           False -> tell [TestResult False parseString "N/A" "N/A"]
+
+           True  -> tell [TestResult True parseString "" "good parse"]
+                      >> react_parser_test_driver mol t xs
+
+
diff --git a/test/check_reversible/Makefile b/test/check_reversible/Makefile
new file mode 100644
--- /dev/null
+++ b/test/check_reversible/Makefile
@@ -0,0 +1,19 @@
+include ../test_helpers/make_common
+
+
+default: check
+
+.PHONY: check clean
+check:
+	ghc $(GHC_FLAGS_DEVEL) -i$(SRC_LOCATION) --make check_deviation.hs
+	ghc $(GHC_FLAGS_DEVEL) -i$(SRC_LOCATION) --make average.hs
+	@echo
+	@echo "Running reversible reaction test."
+	@echo "This may take a while ...."
+	@./run_test.sh $(SIMGI_PATH)
+	@echo
+	
+
+clean:
+	rm -f *.hi *.o check_deviation average reversible.dat \
+		simgi.*
diff --git a/test/check_reversible/average.hs b/test/check_reversible/average.hs
new file mode 100644
--- /dev/null
+++ b/test/check_reversible/average.hs
@@ -0,0 +1,38 @@
+{-----------------------------------------------------------------
+ 
+  (c) 2009 Markus Dittrich 
+ 
+  This program is free software; you can redistribute it 
+  and/or modify it under the terms of the GNU General Public 
+  License Version 3 as published by the Free Software Foundation. 
+ 
+  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 Version 3 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, Inc., 59 Temple Place - Suite 330, 
+  Boston, MA 02111-1307, USA.
+
+--------------------------------------------------------------------}
+
+-- | short helper script for averaging a file containing a single
+-- column of doubles
+module Main where
+
+
+-- imports
+import Prelude
+
+
+main :: IO ()
+main = getContents
+       >>= \content -> 
+         let
+            items     = map (read) . lines $ content :: [Double]
+            theSum    = foldr (+) 0 items
+            average   = theSum/(fromIntegral . length $ items)
+         in
+           putStrLn (show average)
diff --git a/test/check_reversible/check_deviation.hs b/test/check_reversible/check_deviation.hs
new file mode 100644
--- /dev/null
+++ b/test/check_reversible/check_deviation.hs
@@ -0,0 +1,46 @@
+{-----------------------------------------------------------------
+ 
+  (c) 2009 Markus Dittrich 
+ 
+  This program is free software; you can redistribute it 
+  and/or modify it under the terms of the GNU General Public 
+  License Version 3 as published by the Free Software Foundation. 
+ 
+  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 Version 3 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, Inc., 59 Temple Place - Suite 330, 
+  Boston, MA 02111-1307, USA.
+
+--------------------------------------------------------------------}
+
+-- | short helper script for averaging a file containing a single
+-- column of doubles
+module Main where
+
+-- imports
+import Prelude
+import System
+import System.Environment
+
+
+main :: IO Int
+main = getArgs 
+       >>= \(name:value:tol:_) -> readFile name
+       >>= \content -> 
+         let
+            items     = map (read) . lines $ content :: [Double]
+            theSum    = foldr (+) 0 items
+            average   = theSum/(fromIntegral . length $ items)
+            expected  = read value :: Double
+            tolerance = read tol :: Double
+            deviation = abs((average - expected)/expected)
+         in
+           if ( deviation < tolerance )
+             then exitWith ExitSuccess
+             else exitWith $ ExitFailure 1
+           
diff --git a/test/check_reversible/reversible.sgl b/test/check_reversible/reversible.sgl
new file mode 100644
--- /dev/null
+++ b/test/check_reversible/reversible.sgl
@@ -0,0 +1,29 @@
+{-----------------------------------------------------
+
+  this is the input deck for the oregonator model 
+  (C) 2009 Markus Dittrich
+
+------------------------------------------------------}
+
+def parameters
+  time       = 1.0e-4
+  systemVol  = 1.25e-19
+  outputBuffer = 50000
+  outputFreq = 10
+  outputFile = "reversible.dat"
+end
+
+def molecules
+  a = 1000
+  b = 1000
+  c = 0
+end
+
+def reactions
+  a + b -> c  | {1e6*10} |
+  c -> a + b  | 1e5 |
+end
+
+def output
+  [a, b, c]
+end
diff --git a/test/check_reversible/run_test.sh b/test/check_reversible/run_test.sh
new file mode 100644
--- /dev/null
+++ b/test/check_reversible/run_test.sh
@@ -0,0 +1,53 @@
+#!/bin/bash
+#
+# short wrapper for runing a number of reversible
+# simulations and then checking if we can reproduce
+# the correct average product concentration
+#
+
+# 
+simgi_exe="${1}"
+
+# create a global tempfile to collect the data
+globalFile=$(mktemp simgi.XXXXXXXXXXXXX)
+
+for ((counter=0; counter <= 100; counter++)); do
+
+  # provide a little "progressbar"
+  printf "."
+
+  # run and process
+  ${simgi_exe} -s ${counter} reversible.sgl >& /dev/null 
+  tail -n 400 reversible.dat | gawk ' { print $5}' | ./average >> ${globalFile} || return 1
+
+  # unlink
+  rm -f reversible.dat || return 1
+  
+done
+
+# check if we're within the specs
+# the expected number of products is 430.643462709951
+# and we allow 0.5% tolerance
+./check_deviation ${globalFile} 430.643462709951 5e-3
+status=$?
+
+# brief output
+if [[ ${status} == 0 ]]; then
+  echo
+  echo
+  echo "Congratulations - the reversible reaction test passed!"
+  echo
+  echo
+else
+  echo
+  echo
+  echo "Error - the reversible reaction test failed. Please check!"
+  echo
+  echo
+fi
+
+# remove files
+rm -f ${globalFile}
+
+# return status of deviation check
+exit ${status}
diff --git a/test/check_rpnstack/Makefile b/test/check_rpnstack/Makefile
new file mode 100644
--- /dev/null
+++ b/test/check_rpnstack/Makefile
@@ -0,0 +1,15 @@
+include ../test_helpers/make_common
+
+
+default: check
+
+
+.PHONY: check clean
+check:
+	ghc $(GHC_FLAGS_DEVEL) -i$(SRC_LOCATION) -i$(HELPER_LOCATION) \
+		--make RpnStackTest.hs
+	@./RpnStackTest
+
+
+clean:
+	rm -f *.hi *.o RpnStackTest
diff --git a/test/check_rpnstack/RpnStackTest.hs b/test/check_rpnstack/RpnStackTest.hs
new file mode 100644
--- /dev/null
+++ b/test/check_rpnstack/RpnStackTest.hs
@@ -0,0 +1,182 @@
+{-----------------------------------------------------------------
+ 
+  (c) 2009 Markus Dittrich 
+ 
+  This program is free software; you can redistribute it 
+  and/or modify it under the terms of the GNU General Public 
+  License Version 3 as published by the Free Software Foundation. 
+ 
+  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 Version 3 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, Inc., 59 Temple Place - Suite 330, 
+  Boston, MA 02111-1307, USA.
+
+--------------------------------------------------------------------}
+
+-- | this routine tests the functionality in our RPNStack, i.e.,
+-- the infix to RPN parser as well as the compute engine
+module Main where
+
+
+-- imports 
+import Control.Monad.Writer
+import qualified Data.Map as M
+import Prelude
+import System.Exit
+
+-- local imports
+import ExtraFunctions
+import GenericModel
+import PrettyPrint
+import RpnParser
+import RpnCalc
+import TestHelpers
+import TokenParser
+
+
+-----------------------------------------------------------------
+-- 
+-- test data specifications
+--
+-----------------------------------------------------------------
+
+-- | set up the test data
+-- Format : (parse expression, expected result)
+type TestCase = (String, Double)
+
+
+
+-- | simple tests without access to local variables
+simpleTests :: [TestCase]
+simpleTests = 
+  [ ("3+4", 7.0) 
+  , ("3 +4", 7.0)
+  , ("3+  4", 7.0)
+  , ("3 +  4   ", 7.0)
+  , ("-3 + 3   ", 0.0)
+  , ("3*4+4", 16.0)
+  , ("3 * 4+4 ", 16.0)
+  , ("3+4*5", 23.0)
+  , ("(3+4)*5", 35.0)
+  , ("( 3 + 4    )*  5", 35.0)
+  , ("sqrt(2)^2", 2.0)
+  , ("(-2)^2", 4.0)
+  , ("exp(-1)", 0.36787944117144233)
+  , ("log(exp(3)) * 3^2", 27.0)
+  , ("log (  exp ( 3)   ) *3   ^2", 27.0)
+  , ("((((((((3*(3+(3*3) + 4) + 2) +3)-34) -4)+1))))", 16.0)
+  , ("(1*(1 * (1 *(((((3*(3+(3*3) + 4) + 2) +3)-34) -4)+1))))", 16.0)
+  , ("log(exp(2) - sqrt(2))", 1.7875577437560926)
+  , ("2 * 3.14 * sqrt(2)", 8.88126117170303786)
+  ]
+
+
+
+----------------------------------------------------------------
+-- tests with access to local variables and time
+----------------------------------------------------------------
+
+-- | variable tests
+variableTests :: [TestCase]
+variableTests =
+  [ ("3*x", 3000)
+  , ("x", 1000)
+  , ("sqrt(x)^2", 1000)
+  , ("x+y*z", 1000)
+  , ("x + y *   z", 1000)
+  , ("(x+y)*z", 0)
+  , ("(x  + y)* z", 0)
+  , ("exp(z)*TIME", 12.345)
+  , ("exp(z ) *  TIME  ", 12.345)
+  , ("-x * -y", 2.0e6)
+  , ("x*exp(-TIME)", 4.351456244655325e-3)
+  , ("x-x +x -x -y + y", 0.0)
+  , ("-TIME/TIME + TIME - TIME + 1.0", 0.0)
+  ]
+
+
+
+----------------------------------------------------------------
+-- Molecule maps and definitions for specific simulation times
+----------------------------------------------------------------
+
+-- | testmap containing a set of molecules and their 
+-- concentrations
+testMap_1 :: MoleculeMap
+testMap_1 = M.fromList [("x",1000),("y",2000),("z",0)]
+
+-- | a simulation time
+time_1 :: Double
+time_1 = 12.345
+
+
+
+
+----------------------------------------------------------------
+--
+--  main driver routines
+--
+----------------------------------------------------------------
+
+
+-- | main test driver
+main :: IO ()
+main = putStrLn "\n\n\nTesting RPN stack (parser/compute engine)"
+
+  -- run simple tests
+  >> (putStr $ color_string Cyan "\nSimple tests:\n")
+  >> let simpleOut = execWriter $ test_driver 
+           testMap_1 time_1 simpleTests
+     in
+  examine_output simpleOut >>= \simpleStatus ->
+
+
+  -- run variable tests
+  (putStr $ color_string Cyan "\n\nVariable tests:\n")
+  >> let varOut = execWriter $ test_driver testMap_1 time_1
+                               variableTests
+     in
+  examine_output varOut >>= \varStatus ->
+
+
+  -- evaluate status and return
+  let status = simpleStatus && varStatus in
+    if status == True then
+      exitWith ExitSuccess
+    else
+      exitWith $ ExitFailure 1 
+
+
+-- | driver for running a test routine that results in a
+-- successful evaluation of a test expression
+test_driver :: MoleculeMap -> Double -> [TestCase] 
+            -> Writer [TestResult] ()
+test_driver _ _ []          = return ()
+test_driver molMap t (x:xs) =
+
+  let expr     = fst x
+      expected = snd x
+  in
+
+    -- parse expression
+    case runParser parse_infix_to_rpn testModelState "" expr of
+      Left er -> tell [TestResult False expr (show expected) (show er)]
+      Right stack ->
+
+        -- evalute RPN stack
+        let result = rpn_compute (SymbolTable molMap M.empty) t stack in
+          examine_result expected result expr 
+          >> test_driver molMap t xs
+
+          where
+            examine_result target out anExpr = 
+              if is_equal target out
+                then 
+                  tell [TestResult True anExpr (show target) (show out)]
+                else 
+                  tell [TestResult False anExpr (show target) (show out)]
diff --git a/test/reversible_test/Makefile b/test/reversible_test/Makefile
deleted file mode 100644
--- a/test/reversible_test/Makefile
+++ /dev/null
@@ -1,19 +0,0 @@
-
-SRC_LOCATION=../../src
-SIMGI_PATH=${SRC_LOCATION}/simgi
-
-default: check
-
-.PHONY: check clean
-check:
-	ghc -i$(SRC_LOCATION) --make check_deviation.hs
-	ghc -i$(SRC_LOCATION) --make average.hs
-	@echo
-	@echo "Running reversible reaction test."
-	@echo "This may take a while ...."
-	@./run_test.sh ${SIMGI_PATH}
-	@echo
-	
-
-clean:
-	rm -f *.hi *.o check_deviation average
diff --git a/test/reversible_test/average.hs b/test/reversible_test/average.hs
deleted file mode 100644
--- a/test/reversible_test/average.hs
+++ /dev/null
@@ -1,34 +0,0 @@
-{-----------------------------------------------------------------
- 
-  (c) 2009 Markus Dittrich 
- 
-  This program is free software; you can redistribute it 
-  and/or modify it under the terms of the GNU General Public 
-  License Version 3 as published by the Free Software Foundation. 
- 
-  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 Version 3 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, Inc., 59 Temple Place - Suite 330, 
-  Boston, MA 02111-1307, USA.
-
---------------------------------------------------------------------}
-
--- | short helper script for averaging a file containing a single
--- column of doubles
-module Main where
-
-
-main :: IO ()
-main = getContents
-       >>= \content -> 
-         let
-            items     = map (read) . lines $ content :: [Double]
-            theSum    = foldr (+) 0 items
-            average   = theSum/(fromIntegral . length $ items)
-         in
-           putStrLn (show average)
diff --git a/test/reversible_test/check_deviation.hs b/test/reversible_test/check_deviation.hs
deleted file mode 100644
--- a/test/reversible_test/check_deviation.hs
+++ /dev/null
@@ -1,45 +0,0 @@
-{-----------------------------------------------------------------
- 
-  (c) 2009 Markus Dittrich 
- 
-  This program is free software; you can redistribute it 
-  and/or modify it under the terms of the GNU General Public 
-  License Version 3 as published by the Free Software Foundation. 
- 
-  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 Version 3 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, Inc., 59 Temple Place - Suite 330, 
-  Boston, MA 02111-1307, USA.
-
---------------------------------------------------------------------}
-
--- | short helper script for averaging a file containing a single
--- column of doubles
-module Main where
-
--- imports
-import System
-import System.Environment
-
-
-main :: IO Int
-main = getArgs 
-       >>= \(name:value:tol:_) -> readFile name
-       >>= \content -> 
-         let
-            items     = map (read) . lines $ content :: [Double]
-            theSum    = foldr (+) 0 items
-            average   = theSum/(fromIntegral . length $ items)
-            expected  = read value :: Double
-            tolerance = read tol :: Double
-            deviation = abs((average - expected)/expected)
-         in
-           if ( deviation < tolerance )
-             then exitWith ExitSuccess
-             else exitWith $ ExitFailure 1
-           
diff --git a/test/reversible_test/reversible.sgl b/test/reversible_test/reversible.sgl
deleted file mode 100644
--- a/test/reversible_test/reversible.sgl
+++ /dev/null
@@ -1,25 +0,0 @@
-{-----------------------------------------------------
-
-  this is the input deck for the oregonator model 
-  (C) 2009 Markus Dittrich
-
-------------------------------------------------------}
-
-def parameters
-  time       = 1.0e-4
-  systemVol  = 1.25e-19
-  outputIter = 50000
-  outputFreq = 10
-  outputFile = "reversible.dat"
-end
-
-def molecules
-  a 1000
-  b 1000
-  c 0
-end
-
-def reactions
-  a + b -> c  { 1e6*10 }
-  c -> a + b  { 1e5 }
-end 
diff --git a/test/reversible_test/run_test.sh b/test/reversible_test/run_test.sh
deleted file mode 100644
--- a/test/reversible_test/run_test.sh
+++ /dev/null
@@ -1,53 +0,0 @@
-#!/bin/bash
-#
-# short wrapper for runing a number of reversible
-# simulations and then checking if we can reproduce
-# the correct average product concentration
-#
-
-# 
-simgi_exe="${1}"
-
-# create a global tempfile to collect the data
-globalFile=$(mktemp simgi.XXXXXXXXXXXXX)
-
-for ((counter=0; counter <= 100; counter++)); do
-
-  # provide a little "progressbar"
-  printf "."
-
-  # run and process
-  ${simgi_exe} reversible.sgl >& /dev/null 
-  tail -n 400 reversible.dat | gawk ' { print $5}' | ./average >> ${globalFile} || return 1
-
-  # unlink
-  rm -f reversible.dat || return 1
-  
-done
-
-# check if we're within the specs
-# the expected number of products is 430.643462709951
-# and we allow 0.5% tolerance
-./check_deviation ${globalFile} 430.643462709951 5e-3
-status=$?
-
-# brief output
-if [[ ${status} == 0 ]]; then
-  echo
-  echo
-  echo "Congratulations - the reversible reaction test passed!"
-  echo
-  echo
-else
-  echo
-  echo
-  echo "Error - the reversible reaction test failed. Please check!"
-  echo
-  echo
-fi
-
-# remove files
-rm -f ${globalFile}
-
-# return status of deviation check
-exit ${status}
diff --git a/test/rpnstack_test/Makefile b/test/rpnstack_test/Makefile
deleted file mode 100644
--- a/test/rpnstack_test/Makefile
+++ /dev/null
@@ -1,15 +0,0 @@
-
-SRC_LOCATION=../../src
-
-
-default: check
-
-
-.PHONY: check clean
-check:
-	ghc -i$(SRC_LOCATION) --make RpnStackTest.hs
-	@./RpnStackTest
-
-
-clean:
-	rm -f *.hi *.o RpnStackTest
diff --git a/test/rpnstack_test/RpnStackTest.hs b/test/rpnstack_test/RpnStackTest.hs
deleted file mode 100644
--- a/test/rpnstack_test/RpnStackTest.hs
+++ /dev/null
@@ -1,227 +0,0 @@
-{-----------------------------------------------------------------
- 
-  (c) 2009 Markus Dittrich 
- 
-  This program is free software; you can redistribute it 
-  and/or modify it under the terms of the GNU General Public 
-  License Version 3 as published by the Free Software Foundation. 
- 
-  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 Version 3 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, Inc., 59 Temple Place - Suite 330, 
-  Boston, MA 02111-1307, USA.
-
---------------------------------------------------------------------}
-
--- | this routine tests the functionality in our RPNStack, i.e.,
--- the infix to RPN parser as well as the compute engine
-module Main where
-
-
--- imports 
-import Control.Monad.Writer
-import qualified Data.Map as M
-import Prelude
-import System.Exit
-
--- local imports
-import ExtraFunctions
-import GenericModel
-import PrettyPrint
-import RpnParser
-import RpnCalc
-import TokenParser
-
-
------------------------------------------------------------------
--- 
--- test data specifications
---
------------------------------------------------------------------
-
--- | set up the test data
--- Format : (parse expression, expected result)
-type TestCase = (String, Double)
-
-
-
--- | simple tests without access to local variables
-simpleTests :: [TestCase]
-simpleTests = 
-  [ ("3+4", 7.0) 
-  , ("3 +4", 7.0)
-  , ("  3+  4", 7.0)
-  , (" 3 +  4   ", 7.0)
-  , ("-3 + 3   ", 0.0)
-  , ("3*4+4", 16.0)
-  , ("3 * 4+4 ", 16.0)
-  , ("3+4*5", 23.0)
-  , ("(3+4)*5", 35.0)
-  , (" ( 3 + 4    )*  5", 35.0)
-  , ("sqrt(2)^2", 2.0)
-  , ("(-2)^2", 4.0)
-  , ("exp(-1)", 0.36787944117144233)
-  , ("log(exp(3)) * 3^2", 27.0)
-  , ("log (  exp ( 3)   ) *3   ^2", 27.0)
-  , ("((((((((3*(3+(3*3) + 4) + 2) +3)-34) -4)+1))))", 16.0)
-  , ("(1*(1 * (1 *(((((3*(3+(3*3) + 4) + 2) +3)-34) -4)+1))))", 16.0)
-  , ("log(exp(2) - sqrt(2))", 1.7875577437560926)
-  , ("2 * 3.14 * sqrt(2)", 8.88126117170303786)
-  ]
-
-
--- | variable tests
-variableTests :: [TestCase]
-variableTests =
-  [ ("3*x", 3000)
-  , ("sqrt(x)^2", 1000)
-  , ("x+y*z", 1000)
-  , ("x + y *   z", 1000)
-  , ("(x+y)*z", 0)
-  , ("(x  + y)* z", 0)
-  , ("exp(z)*TIME", 12.345)
-  , ("exp(z ) *  TIME  ", 12.345)
-  , ("-x * -y", 2.0e6)
-  , ("x*exp(-TIME)", 4.351456244655325e-3)
-  , ("x-x +x -x -y + y", 0.0)
-  , ("-TIME/TIME + TIME - TIME + 1.0", 0.0)
-  ]
-
-
-----------------------------------------------------------------
--- tests with access to local variables and time
-----------------------------------------------------------------
-
-
-
-----------------------------------------------------------------
--- Molecule maps and definitions for specific simulation times
-----------------------------------------------------------------
-
--- | testmap containing a set of molecules and their 
--- concentrations
-testMap_1 :: MoleculeMap
-testMap_1 = M.fromList [("x",1000),("y",2000),("z",0)]
-
--- | a simulation time
-time_1 :: Double
-time_1 = 12.345
-
-
-
-
-----------------------------------------------------------------
---
---  main driver routines
---
-----------------------------------------------------------------
-
-
--- | main test driver
-main :: IO ()
-main = putStrLn "\n\n\nTesting RPN stack (parser/compute engine)"
-
-  -- run simple tests
-  >> (putStr $ color_string Cyan "\nSimple tests:\n")
-  >> let simpleOut = execWriter $ test_driver 
-           testMap_1 time_1 simpleTests
-     in
-  examine_output simpleOut >>= \simpleStatus ->
-
-
-  -- run variable tests
-  (putStr $ color_string Cyan "\n\nVariable tests:\n")
-  >> let varOut = execWriter $ test_driver testMap_1 time_1
-                               variableTests
-     in
-  examine_output varOut >>= \varStatus ->
-
-
-  -- evaluate status and return
-  let status = simpleStatus && varStatus in
-    if status == True then
-      exitWith ExitSuccess
-    else
-      exitWith $ ExitFailure 1 
-
-
--- | examine the output of a test routine
--- | helper function for examining the output of a good test run
--- (i.e. one that should succeed), prints out the result for each 
--- test, collects the number of successes/failures and returns 
--- True in case all tests succeeded and False otherwise
-examine_output :: [TestResult] -> IO Bool
-examine_output = foldM examine_output_h True
-                 
-  where
-    examine_output_h :: Bool -> TestResult -> IO Bool
-    examine_output_h acc (TestResult status token target actual) = do
-      if status == True then do
-          putStr   $ color_string Blue "["
-          putStr   $ color_string White "OK"
-          putStr   $ color_string Blue  "] "
-          putStr   $ color_string Green " Successfully evaluated "
-          putStrLn $ color_string Yellow token
-          return $ acc && True
-        else do
-          putStr   $ color_string Blue "["
-          putStr   $ color_string Red "TROUBLE"
-          putStr   $ color_string Blue "] "
-          putStr   $ color_string Green " Failed to evaluate "
-          putStrLn $ color_string Yellow token
-          putStrLn $ color_string Green "\t\texpected : " 
-                       ++ (show target)
-          putStrLn $ color_string Green "\t\tgot      : " 
-                       ++ (show actual)
-          return False
-
-
--- | driver for running a test routine that results in a
--- successful evaluation of a test expression
-test_driver :: MoleculeMap -> Double -> [TestCase] 
-            -> Writer [TestResult] ()
-test_driver _ _ []          = return ()
-test_driver mol time (x:xs) =
-
-  let expr     = fst x
-      expected = snd x
-  in
-
-    -- parse expression
-    case runParser parse_infix_to_rpn initialModelState "" expr of
-      Left er -> tell [TestResult False expr (show expected) (show er)]
-      Right stack ->
-
-        -- evalute RPN stack
-        let result = rpn_compute mol time stack in
-          examine_result expected result expr 
-          >> test_driver mol time xs
-
-          where
-            examine_result target out expr = 
-              if is_equal target out
-                then 
-                  tell [TestResult True expr (show target) (show out)]
-                else 
-                  tell [TestResult False expr (show target) (show out)]
-          
-
--- | data structure for keeping track of our test results
--- which consist of a bool indicating success
--- or failure, the test token as well as the expected and
--- received result
-data TestResult = TestResult { status :: Bool
-                             , token  :: String
-                             , target :: String
-                             , actual :: String
-                             }
-
-
-defaultResult :: TestResult
-defaultResult = TestResult False "" "" ""
-
diff --git a/test/test_helpers/TestHelpers.hs b/test/test_helpers/TestHelpers.hs
new file mode 100644
--- /dev/null
+++ b/test/test_helpers/TestHelpers.hs
@@ -0,0 +1,110 @@
+{-----------------------------------------------------------------
+ 
+  (c) 2009 Markus Dittrich 
+ 
+  This program is free software; you can redistribute it 
+  and/or modify it under the terms of the GNU General Public 
+  License Version 3 as published by the Free Software Foundation. 
+ 
+  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 Version 3 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, Inc., 59 Temple Place - Suite 330, 
+  Boston, MA 02111-1307, USA.
+
+--------------------------------------------------------------------}
+
+-- | this module provides a few common routines used in our
+-- unit tests
+module TestHelpers ( examine_output
+                   , testModelState
+                   , TestResult(TestResult)
+                   , defaultResult
+                   ) where
+
+
+-- imports 
+import Control.Monad
+import qualified Data.Map as M
+import Prelude
+
+-- local imports
+import GenericModel
+import PrettyPrint
+
+
+
+-- | examine the output of a test routine
+-- | helper function for examining the output of a good test run
+-- (i.e. one that should succeed), prints out the result for each 
+-- test, collects the number of successes/failures and returns 
+-- True in case all tests succeeded and False otherwise
+examine_output :: [TestResult] -> IO Bool
+examine_output = foldM examine_output_h True
+                 
+  where
+    examine_output_h :: Bool -> TestResult -> IO Bool
+    examine_output_h acc (TestResult state tok targ act) = do
+      if state == True then do
+          putStr   $ color_string Blue "["
+          putStr   $ color_string White "OK"
+          putStr   $ color_string Blue  "] "
+          putStr   $ color_string Green " Successfully evaluated "
+          putStrLn $ color_string Yellow tok
+          return $ acc && True
+        else do
+          putStr   $ color_string Blue "["
+          putStr   $ color_string Red "TROUBLE"
+          putStr   $ color_string Blue "] "
+          putStr   $ color_string Green " Failed to evaluate "
+          putStrLn $ color_string Yellow tok
+          putStrLn $ color_string Green "\t\texpected : " 
+                       ++ (show targ)
+          putStrLn $ color_string Green "\t\tgot      : " 
+                       ++ (show act)
+          return False
+
+
+
+-- | data structure for keeping track of our test results
+-- which consist of a bool indicating success
+-- or failure, the test token as well as the expected and
+-- received result
+data TestResult = TestResult { status :: Bool
+                             , token  :: String
+                             , target :: String
+                             , actual :: String
+                             }
+
+
+defaultResult :: TestResult
+defaultResult = TestResult { status = False 
+                           , token  = ""
+                           , target = ""
+                           , actual = ""
+                           }
+
+
+-- | initial model state we use for our unit tests
+-- we use a negative system volume so we don't have
+-- to deal with unit conversion when testing the
+-- reaction parser
+testModelState :: ModelState
+testModelState = ModelState { molCount         = M.empty
+                            , rates            = []
+                            , reactions        = []
+                            , randNums         = []
+                            , events           = []
+                            , systemVol        = -1.0 -- negative!
+                            , currentTime      = 0.0
+                            , currentIter      = 0
+                            , maxTime          = 0.0
+                            , outputBufferSize = 10000
+                            , outputFreq       = 1000
+                            , outputCache      = []
+                            , outfileName      = ""
+                            }
diff --git a/test/test_helpers/make_common b/test/test_helpers/make_common
new file mode 100644
--- /dev/null
+++ b/test/test_helpers/make_common
@@ -0,0 +1,7 @@
+# some common variables all Makefiles in the test directory use
+
+GHC_FLAGS_DEVEL = -O -Wall -fwarn-simple-patterns -fwarn-tabs -fwarn-incomplete-record-updates -fwarn-monomorphism-restriction -fwarn-implicit-prelude -fno-warn-orphans
+
+SRC_LOCATION=../../src
+HELPER_LOCATION=../test_helpers
+SIMGI_PATH=$(SRC_LOCATION)/simgi