diff --git a/Simulation/Aivika/Dynamics.hs b/Simulation/Aivika/Dynamics.hs
--- a/Simulation/Aivika/Dynamics.hs
+++ b/Simulation/Aivika/Dynamics.hs
@@ -1,1628 +1,35 @@
 
--- Copyright (c) 2009, 2010, 2011 David Sorokin <david.sorokin@gmail.com>
--- 
--- All rights reserved.
--- 
--- Redistribution and use in source and binary forms, with or without
--- modification, are permitted provided that the following conditions
--- are met:
--- 
--- 1. Redistributions of source code must retain the above copyright
---    notice, this list of conditions and the following disclaimer.
--- 
--- 2. Redistributions in binary form must reproduce the above copyright
---    notice, this list of conditions and the following disclaimer in the
---    documentation and/or other materials provided with the distribution.
--- 
--- 3. Neither the name of the author nor the names of his contributors
---    may be used to endorse or promote products derived from this software
---    without specific prior written permission.
--- 
--- THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND
--- ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
--- IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
--- ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE
--- FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
--- DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
--- OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
--- HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
--- LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
--- OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
--- SUCH DAMAGE.
-
-{-# LANGUAGE FlexibleContexts, FlexibleInstances, UndecidableInstances #-}
-
--- |
--- Module     : Simulation.Aivika.Dynamics
--- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>
--- License    : BSD3
--- Maintainer : David Sorokin <david.sorokin@gmail.com>
--- Stability  : experimental
--- Tested with: GHC 6.12.1
---
--- Aivika is a multi-paradigm simulation library. It allows us to integrate 
--- a system of ordinary differential equations. Also it can be applied to
--- the Discrete Event Simulation. It supports the event-oriented, 
--- process-oriented and activity-oriented paradigms. Aivika also supports 
--- the Agent-based Modeling. Finally, it can be applied to System Dynamics.
---
-module Simulation.Aivika.Dynamics 
-       (-- * Dynamics
-        Dynamics,
-        DynamicsTrans(..),
-        Specs(..),
-        Method(..),
-        runDynamics1,
-        runDynamics,
-        runDynamicsIO,
-        -- ** Time parameters
-        starttime,
-        stoptime,
-        dt,
-        time,
-        -- ** Maximum and Minimum
-        maxD,
-        minD,
-        -- ** Integrals
-        Integ,
-        newInteg,
-        integInit,
-        integValue,
-        integDiff,
-        -- ** Table Functions
-        lookupD,
-        lookupStepwiseD,
-        -- ** Interpolation
-        initD,
-        discrete,
-        interpolate,
-        -- ** Memoization and Sequential Calculations
-        Memo,
-        UMemo,
-        memo,
-        umemo,
-        memo0,
-        umemo0,
-        -- ** Utility
-        once,
-        -- * Event Queue
-        DynamicsQueue,
-        newQueue,
-        enqueueDC,
-        enqueueD,
-        runQueue,
-        -- * References
-        DynamicsRef,
-        newRef,
-        refQueue,
-        readRef,
-        writeRef,
-        writeRef',
-        modifyRef,
-        modifyRef',
-        -- * Discontinuous Processes
-        DynamicsPID,
-        DynamicsProc,
-        newPID,
-        pidQueue,
-        holdProcD,
-        holdProc,
-        passivateProc,
-        procPassive,
-        reactivateProc,
-        procPID,
-        runProc,
-        -- * Resources
-        DynamicsResource,
-        newResource,
-        resourceQueue,
-        resourceInitCount,
-        resourceCount,
-        requestResource,
-        releaseResource,
-        -- * Agent-based Modeling
-        Agent,
-        AgentState,
-        newAgent,
-        newState,
-        newSubstate,
-        agentQueue,
-        agentState,
-        activateState,
-        initState,
-        stateAgent,
-        stateParent,
-        addTimeoutD,
-        addTimeout,
-        addTimerD,
-        addTimer,
-        stateActivation,
-        stateDeactivation) where
-
-import Data.Array
-import Data.Array.IO
-import Data.IORef
-import Control.Monad
-import Control.Monad.Trans
-
-import qualified Simulation.Aivika.Queue as Q
-import qualified Simulation.Aivika.PriorityQueue as PQ
-
---
--- The Dynamics Monad
---
--- A value of the Dynamics monad represents an abstract dynamic 
--- process, i.e. a time varying polymorphic function. This is 
--- a key point of the Aivika simulation library.
---
-
--- | A value in the 'Dynamics' monad represents a dynamic process, i.e.
--- a polymorphic time varying function.
-newtype Dynamics a = Dynamics (Parameters -> IO a)
-
--- | It defines the simulation time appended with additional information.
-data Parameters = Parameters { parSpecs :: Specs,    -- ^ the simulation specs
-                               parTime :: Double,    -- ^ the current time
-                               parIteration :: Int,  -- ^ the current iteration
-                               parPhase :: Int }     -- ^ the current phase
-
--- | It defines the simulation specs.
-data Specs = Specs { spcStartTime :: Double,    -- ^ the start time
-                     spcStopTime :: Double,     -- ^ the stop time
-                     spcDT :: Double,           -- ^ the integration time step
-                     spcMethod :: Method        -- ^ the integration method
-                   } deriving (Eq, Ord, Show)
-
--- | It defines the integration method.
-data Method = Euler          -- ^ Euler's method
-            | RungeKutta2    -- ^ the 2nd order Runge-Kutta method
-            | RungeKutta4    -- ^ the 4th order Runge-Kutta method
-            deriving (Eq, Ord, Show)
-
-iterations :: Specs -> [Int]
-iterations sc = [i1 .. i2] where
-  i1 = 0
-  i2 = round ((spcStopTime sc - 
-               spcStartTime sc) / spcDT sc)
-
-iterationBnds :: Specs -> (Int, Int)
-iterationBnds sc = (0, round ((spcStopTime sc - 
-                               spcStartTime sc) / spcDT sc))
-
-iterationLoBnd :: Specs -> Int
-iterationLoBnd sc = 0
-
-iterationHiBnd :: Specs -> Int
-iterationHiBnd sc = round ((spcStopTime sc - 
-                            spcStartTime sc) / spcDT sc)
-
-phases :: Specs -> [Int]
-phases sc = 
-  case spcMethod sc of
-    Euler -> [0]
-    RungeKutta2 -> [0, 1]
-    RungeKutta4 -> [0, 1, 2, 3]
-
-phaseBnds :: Specs -> (Int, Int)
-phaseBnds sc = 
-  case spcMethod sc of
-    Euler -> (0, 0)
-    RungeKutta2 -> (0, 1)
-    RungeKutta4 -> (0, 3)
-
-phaseLoBnd :: Specs -> Int
-phaseLoBnd sc = 0
-                  
-phaseHiBnd :: Specs -> Int
-phaseHiBnd sc = 
-  case spcMethod sc of
-    Euler -> 0
-    RungeKutta2 -> 1
-    RungeKutta4 -> 3
-
-basicTime :: Specs -> Int -> Int -> Double
-basicTime sc n ph =
-  if ph < 0 then 
-    error "Incorrect phase: basicTime"
-  else
-    spcStartTime sc + n' * spcDT sc + delta (spcMethod sc) ph 
-      where n' = fromInteger (toInteger n)
-            delta Euler       0 = 0
-            delta RungeKutta2 0 = 0
-            delta RungeKutta2 1 = spcDT sc
-            delta RungeKutta4 0 = 0
-            delta RungeKutta4 1 = spcDT sc / 2
-            delta RungeKutta4 2 = spcDT sc / 2
-            delta RungeKutta4 3 = spcDT sc
-
-neighborhood :: Specs -> Double -> Double -> Bool
-neighborhood sc t t' = 
-  abs (t - t') <= spcDT sc / 1.0e6
-
-instance Monad Dynamics where
-  return  = returnD
-  m >>= k = bindD m k
-
-returnD :: a -> Dynamics a
-returnD a = Dynamics (\ps -> return a)
-
-bindD :: Dynamics a -> (a -> Dynamics b) -> Dynamics b
-bindD (Dynamics m) k = 
-  Dynamics $ \ps -> 
-  do a <- m ps
-     let Dynamics m' = k a
-     m' ps
-
-subrunDynamics1 :: Dynamics a -> Specs -> IO a
-subrunDynamics1 (Dynamics m) sc =
-  do let n = iterationHiBnd sc
-         t = basicTime sc n 0
-     m Parameters { parSpecs = sc,
-                    parTime = t,
-                    parIteration = n,
-                    parPhase = 0 }
-
-subrunDynamics :: Dynamics a -> Specs -> [IO a]
-subrunDynamics (Dynamics m) sc =
-  do let (nl, nu) = iterationBnds sc
-         parameterise n = Parameters { parSpecs = sc,
-                                       parTime = basicTime sc n 0,
-                                       parIteration = n,
-                                       parPhase = 0 }
-     map (m . parameterise) [nl .. nu]
-
--- | Run the simulation and return the result in the last 
--- time point using the specified simulation specs.
-runDynamics1 :: Dynamics (Dynamics a) -> Specs -> IO a
-runDynamics1 (Dynamics m) sc = 
-  do d <- m Parameters { parSpecs = sc,
-                         parTime = spcStartTime sc,
-                         parIteration = 0,
-                         parPhase = 0 }
-     subrunDynamics1 d sc
-
--- | Run the simulation and return the results in all 
--- integration time points using the specified simulation specs.
-runDynamics :: Dynamics (Dynamics a) -> Specs -> IO [a]
-runDynamics (Dynamics m) sc = 
-  do d <- m Parameters { parSpecs = sc,
-                         parTime = spcStartTime sc,
-                         parIteration = 0,
-                         parPhase = 0 }
-     sequence $ subrunDynamics d sc
-
--- | Run the simulation and return the results in all 
--- integration time points using the specified simulation specs.
-runDynamicsIO :: Dynamics (Dynamics a) -> Specs -> IO [IO a]
-runDynamicsIO (Dynamics m) sc =
-  do d <- m Parameters { parSpecs = sc,
-                         parTime = spcStartTime sc,
-                         parIteration = 0,
-                         parPhase = 0 }
-     return $ subrunDynamics d sc
-
-instance Functor Dynamics where
-  fmap f (Dynamics m) = 
-    Dynamics $ \ps -> do { a <- m ps; return $ f a }
-
-instance Eq (Dynamics a) where
-  x == y = error "Can't compare dynamics." 
-
-instance Show (Dynamics a) where
-  showsPrec _ x = showString "<< Dynamics >>"
-
-liftMD :: (a -> b) -> Dynamics a -> Dynamics b
-liftMD f (Dynamics x) =
-  Dynamics $ \ps -> do { a <- x ps; return $ f a }
-
-liftM2D :: (a -> b -> c) -> Dynamics a -> Dynamics b -> Dynamics c
-liftM2D f (Dynamics x) (Dynamics y) =
-  Dynamics $ \ps -> do { a <- x ps; b <- y ps; return $ f a b }
-
-instance (Num a) => Num (Dynamics a) where
-  x + y = liftM2D (+) x y
-  x - y = liftM2D (-) x y
-  x * y = liftM2D (*) x y
-  negate = liftMD negate
-  abs = liftMD abs
-  signum = liftMD signum
-  fromInteger i = return $ fromInteger i
-
-instance (Fractional a) => Fractional (Dynamics a) where
-  x / y = liftM2D (/) x y
-  recip = liftMD recip
-  fromRational t = return $ fromRational t
-
-instance (Floating a) => Floating (Dynamics a) where
-  pi = return pi
-  exp = liftMD exp
-  log = liftMD log
-  sqrt = liftMD sqrt
-  x ** y = liftM2D (**) x y
-  sin = liftMD sin
-  cos = liftMD cos
-  tan = liftMD tan
-  asin = liftMD asin
-  acos = liftMD acos
-  atan = liftMD atan
-  sinh = liftMD sinh
-  cosh = liftMD cosh
-  tanh = liftMD tanh
-  asinh = liftMD asinh
-  acosh = liftMD acosh
-  atanh = liftMD atanh
-
-instance MonadIO Dynamics where
-  liftIO m = Dynamics $ const m
-
--- | The 'DynamicsTrans' class defines a type which the 'Dynamics' 
--- computation can be lifted to.
-class DynamicsTrans m where
-  -- | Lift the computation.
-  liftD :: Dynamics a -> m a
-
---
--- Integration Parameters and Time
---
-
--- | Return the start simulation time.
-starttime :: Dynamics Double
-starttime = Dynamics $ return . spcStartTime . parSpecs
-
--- | Return the stop simulation time.
-stoptime :: Dynamics Double
-stoptime = Dynamics $ return . spcStopTime . parSpecs
-
--- | Return the integration time step.
-dt :: Dynamics Double
-dt = Dynamics $ return . spcDT . parSpecs
-
--- | Return the current simulation time.
-time :: Dynamics Double
-time = Dynamics $ return . parTime 
-
---
--- Maximum and Minimum
---
-
--- | Return the maximum.
-maxD :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics a
-maxD = liftM2D max
-
--- | Return the minimum.
-minD :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics a
-minD = liftM2D min
-           
---
--- System Dynamics
---
-
--- | The 'Integ' type represents an integral.
-data Integ = Integ { integInit     :: Dynamics Double,   -- ^ The initial value.
-                     integExternal :: IORef (Dynamics Double),
-                     integInternal :: IORef (Dynamics Double) }
-
--- | Create a new integral with the specified initial value.
-newInteg :: Dynamics Double -> Dynamics Integ
-newInteg i = 
-  do r1 <- liftIO $ newIORef $ initD i 
-     r2 <- liftIO $ newIORef $ initD i 
-     let integ = Integ { integInit     = i, 
-                         integExternal = r1,
-                         integInternal = r2 }
-         z = Dynamics $ \ps -> 
-           do (Dynamics m) <- readIORef (integInternal integ)
-              m ps
-     y <- umemo interpolate z
-     liftIO $ writeIORef (integExternal integ) y
-     return integ
-
--- | Return the integral's value.
-integValue :: Integ -> Dynamics Double
-integValue integ = 
-  Dynamics $ \ps ->
-  do (Dynamics m) <- readIORef (integExternal integ)
-     m ps
-
--- | Set the derivative for the integral.
-integDiff :: Integ -> Dynamics Double -> Dynamics ()
-integDiff integ diff =
-  do let z = Dynamics $ \ps ->
-           do y <- readIORef (integExternal integ)
-              let i = integInit integ
-              case spcMethod (parSpecs ps) of
-                Euler -> integEuler diff i y ps
-                RungeKutta2 -> integRK2 diff i y ps
-                RungeKutta4 -> integRK4 diff i y ps
-     liftIO $ writeIORef (integInternal integ) z
-
-integEuler :: Dynamics Double
-             -> Dynamics Double 
-             -> Dynamics Double 
-             -> Parameters -> IO Double
-integEuler (Dynamics f) (Dynamics i) (Dynamics y) ps = 
-  case parIteration ps of
-    0 -> 
-      i ps
-    n -> do 
-      let sc  = parSpecs ps
-          ty  = basicTime sc (n - 1) 0
-          psy = ps { parTime = ty, parIteration = n - 1, parPhase = 0 }
-      a <- y psy
-      b <- f psy
-      let !v = a + spcDT (parSpecs ps) * b
-      return v
-
-integRK2 :: Dynamics Double
-           -> Dynamics Double
-           -> Dynamics Double
-           -> Parameters -> IO Double
-integRK2 (Dynamics f) (Dynamics i) (Dynamics y) ps =
-  case parPhase ps of
-    0 -> case parIteration ps of
-      0 ->
-        i ps
-      n -> do
-        let sc = parSpecs ps
-            ty = basicTime sc (n - 1) 0
-            t1 = ty
-            t2 = basicTime sc (n - 1) 1
-            psy = ps { parTime = ty, parIteration = n - 1, parPhase = 0 }
-            ps1 = psy
-            ps2 = ps { parTime = t2, parIteration = n - 1, parPhase = 1 }
-        vy <- y psy
-        k1 <- f ps1
-        k2 <- f ps2
-        let !v = vy + spcDT sc / 2.0 * (k1 + k2)
-        return v
-    1 -> do
-      let sc = parSpecs ps
-          n  = parIteration ps
-          ty = basicTime sc n 0
-          t1 = ty
-          psy = ps { parTime = ty, parIteration = n, parPhase = 0 }
-          ps1 = psy
-      vy <- y psy
-      k1 <- f ps1
-      let !v = vy + spcDT sc * k1
-      return v
-    _ -> 
-      error "Incorrect phase: integ"
-
-integRK4 :: Dynamics Double
-           -> Dynamics Double
-           -> Dynamics Double
-           -> Parameters -> IO Double
-integRK4 (Dynamics f) (Dynamics i) (Dynamics y) ps =
-  case parPhase ps of
-    0 -> case parIteration ps of
-      0 -> 
-        i ps
-      n -> do
-        let sc = parSpecs ps
-            ty = basicTime sc (n - 1) 0
-            t1 = ty
-            t2 = basicTime sc (n - 1) 1
-            t3 = basicTime sc (n - 1) 2
-            t4 = basicTime sc (n - 1) 3
-            psy = ps { parTime = ty, parIteration = n - 1, parPhase = 0 }
-            ps1 = psy
-            ps2 = ps { parTime = t2, parIteration = n - 1, parPhase = 1 }
-            ps3 = ps { parTime = t3, parIteration = n - 1, parPhase = 2 }
-            ps4 = ps { parTime = t4, parIteration = n - 1, parPhase = 3 }
-        vy <- y psy
-        k1 <- f ps1
-        k2 <- f ps2
-        k3 <- f ps3
-        k4 <- f ps4
-        let !v = vy + spcDT sc / 6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
-        return v
-    1 -> do
-      let sc = parSpecs ps
-          n  = parIteration ps
-          ty = basicTime sc n 0
-          t1 = ty
-          psy = ps { parTime = ty, parIteration = n, parPhase = 0 }
-          ps1 = psy
-      vy <- y psy
-      k1 <- f ps1
-      let !v = vy + spcDT sc / 2.0 * k1
-      return v
-    2 -> do
-      let sc = parSpecs ps
-          n  = parIteration ps
-          ty = basicTime sc n 0
-          t2 = basicTime sc n 1
-          psy = ps { parTime = ty, parIteration = n, parPhase = 0 }
-          ps2 = ps { parTime = t2, parIteration = n, parPhase = 1 }
-      vy <- y psy
-      k2 <- f ps2
-      let !v = vy + spcDT sc / 2.0 * k2
-      return v
-    3 -> do
-      let sc = parSpecs ps
-          n  = parIteration ps
-          ty = basicTime sc n 0
-          t3 = basicTime sc n 2
-          psy = ps { parTime = ty, parIteration = n, parPhase = 0 }
-          ps3 = ps { parTime = t3, parIteration = n, parPhase = 2 }
-      vy <- y psy
-      k3 <- f ps3
-      let !v = vy + spcDT sc * k3
-      return v
-    _ -> 
-      error "Incorrect phase: integ"
-
--- smoothI :: Dynamics Double -> Dynamics Double -> Dynamics Double 
---           -> Dynamics Double
--- smoothI x t i = y where
---   y = integ ((x - y) / t) i
-
--- smooth :: Dynamics Double -> Dynamics Double -> Dynamics Double
--- smooth x t = smoothI x t x
-
--- smooth3I :: Dynamics Double -> Dynamics Double -> Dynamics Double 
---            -> Dynamics Double
--- smooth3I x t i = y where
---   y  = integ ((s1 - y) / t') i
---   s1 = integ ((s0 - s1) / t') i
---   s0 = integ ((x - s0) / t') i
---   t' = t / 3.0
-
--- smooth3 :: Dynamics Double -> Dynamics Double -> Dynamics Double
--- smooth3 x t = smooth3I x t x
-
--- smoothNI :: Dynamics Double -> Dynamics Double -> Int -> Dynamics Double 
---            -> Dynamics Double
--- smoothNI x t n i = s ! n where
---   s   = array (1, n) [(k, f k) | k <- [1 .. n]]
---   f 0 = integ ((x - s ! 0) / t') i
---   f k = integ ((s ! (k - 1) - s ! k) / t') i
---   t'  = t / fromIntegral n
-
--- smoothN :: Dynamics Double -> Dynamics Double -> Int -> Dynamics Double
--- smoothN x t n = smoothNI x t n x
-
--- delay1I :: Dynamics Double -> Dynamics Double -> Dynamics Double 
---           -> Dynamics Double
--- delay1I x t i = y where
---   y = integ (x - y) (i * t) / t
-
--- delay1 :: Dynamics Double -> Dynamics Double -> Dynamics Double
--- delay1 x t = delay1I x t x
-
--- delay3I :: Dynamics Double -> Dynamics Double -> Dynamics Double 
---           -> Dynamics Double
--- delay3I x t i = y where
---   y  = integ (s1 - y) (i * t') / t'
---   s1 = integ (s0 - s1) (i * t') / t'
---   s0 = integ (x - s0) (i * t') / t'
---   t' = t / 3.0
-
--- delay3 :: Dynamics Double -> Dynamics Double -> Dynamics Double
--- delay3 x t = delay3I x t x
-
--- delayNI :: Dynamics Double -> Dynamics Double -> Int -> Dynamics Double 
---           -> Dynamics Double
--- delayNI x t n i = s ! n where
---   s   = array (1, n) [(k, f k) | k <- [1 .. n]]
---   f 0 = integ (x - s ! 0) (i * t') / t'
---   f k = integ (s ! (k - 1) - s ! k) (i * t') / t'
---   t'  = t / fromIntegral n
-
--- delayN :: Dynamics Double -> Dynamics Double -> Int -> Dynamics Double
--- delayN x t n = delayNI x t n x
-
--- forecast :: Dynamics Double -> Dynamics Double -> Dynamics Double 
---            -> Dynamics Double
--- forecast x at hz =
---   x * (1.0 + (x / smooth x at - 1.0) / at * hz)
-
--- trend :: Dynamics Double -> Dynamics Double -> Dynamics Double 
---         -> Dynamics Double
--- trend x at i =
---   (x / smoothI x at (x / (1.0 + i * at)) - 1.0) / at
-
---
--- Table Functions
---
-
--- | Lookup @x@ in a table of pairs @(x, y)@ using linear interpolation.
-lookupD :: Dynamics Double -> Array Int (Double, Double) -> Dynamics Double
-lookupD (Dynamics m) tbl =
-  Dynamics (\ps -> do a <- m ps; return $ find first last a) where
-    (first, last) = bounds tbl
-    find left right x =
-      if left > right then
-        error "Incorrect index: table"
-      else
-        let index = (left + 1 + right) `div` 2
-            x1    = fst $ tbl ! index
-        in if x1 <= x then 
-             let y | index < right = find index right x
-                   | right == last  = snd $ tbl ! right
-                   | otherwise     = 
-                     let x2 = fst $ tbl ! (index + 1)
-                         y1 = snd $ tbl ! index
-                         y2 = snd $ tbl ! (index + 1)
-                     in y1 + (y2 - y1) * (x - x1) / (x2 - x1) 
-             in y
-           else
-             let y | left < index = find left (index - 1) x
-                   | left == first = snd $ tbl ! left
-                   | otherwise    = error "Incorrect index: table"
-             in y
-
--- | Lookup @x@ in a table of pairs @(x, y)@ using stepwise function.
-lookupStepwiseD :: Dynamics Double -> Array Int (Double, Double)
-                  -> Dynamics Double
-lookupStepwiseD (Dynamics m) tbl =
-  Dynamics (\ps -> do a <- m ps; return $ find first last a) where
-    (first, last) = bounds tbl
-    find left right x =
-      if left > right then
-        error "Incorrect index: table"
-      else
-        let index = (left + 1 + right) `div` 2
-            x1    = fst $ tbl ! index
-        in if x1 <= x then 
-             let y | index < right = find index right x
-                   | right == last  = snd $ tbl ! right
-                   | otherwise     = snd $ tbl ! right
-             in y
-           else
-             let y | left < index = find left (index - 1) x
-                   | left == first = snd $ tbl ! left
-                   | otherwise    = error "Incorrect index: table"
-             in y
-
--- --
--- -- Discrete Functions
--- --
-    
--- delayTrans :: Dynamics a -> Dynamics Double -> Dynamics a 
---               -> (Dynamics a -> Dynamics a) -> Dynamics a
--- delayTrans (Dynamics x) (Dynamics d) (Dynamics i) tr = tr $ Dynamics r 
---   where
---     r ps = do 
---       let t  = parTime ps
---           sc = parSpecs ps
---           n  = parIteration ps
---       a <- d ps
---       let t' = (t - a) - spcStartTime sc
---           n' = fromInteger $ toInteger $ floor $ t' / spcDT sc
---           y | n' < 0    = i $ ps { parTime = spcStartTime sc,
---                                    parIteration = 0, 
---                                    parPhase = 0 }
---             | n' < n    = x $ ps { parTime = t',
---                                    parIteration = n',
---                                    parPhase = -1 }
---             | n' > n    = error "Cannot return the future data: delay"
---             | otherwise = error "Cannot return the current data: delay"
---       y    
-
--- delay :: (Memo a) => Dynamics a -> Dynamics Double -> Dynamics a
--- delay x d = delayTrans x d x $ memo0 discrete
-
--- delay' :: (UMemo a) => Dynamics a -> Dynamics Double -> Dynamics a
--- delay' x d = delayTrans x d x $ memo0' discrete
-
--- delayI :: (Memo a) => Dynamics a -> Dynamics Double -> Dynamics a -> Dynamics a
--- delayI x d i = delayTrans x d i $ memo0 discrete         
-
--- delayI' :: (UMemo a) => Dynamics a -> Dynamics Double -> Dynamics a -> Dynamics a
--- delayI' x d i = delayTrans x d i $ memo0' discrete         
-
---
--- Interpolation and Initial Value
---
--- These functions complement the memoization, possibly except for 
--- the initial function which can be also useful to get an initial 
--- value of any dynamic process. See comments to the Memoization 
--- section.
---      
-
--- | Return the initial value.
-initD :: Dynamics a -> Dynamics a
-initD (Dynamics m) =
-  Dynamics $ \ps ->
-  if parIteration ps == 0 && parPhase ps == 0 then
-    m ps
-  else
-    let sc = parSpecs ps
-    in m $ ps { parTime = basicTime sc 0 0,
-                parIteration = 0,
-                parPhase = 0 } 
-
--- | Discretize the computation in the integration time points.
-discrete :: Dynamics a -> Dynamics a
-discrete (Dynamics m) =
-  Dynamics $ \ps ->
-  let ph = parPhase ps
-      r | ph == 0    = m ps
-        | ph > 0    = let sc = parSpecs ps
-                          n  = parIteration ps
-                      in m $ ps { parTime = basicTime sc n 0,
-                                  parPhase = 0 }
-        | otherwise = let sc = parSpecs ps
-                          t  = parTime ps
-                          n  = parIteration ps
-                          t' = spcStartTime sc + fromIntegral (n + 1) * spcDT sc
-                          n' = if neighborhood sc t t' then n + 1 else n
-                      in m $ ps { parTime = basicTime sc n' 0,
-                                  parIteration = n',
-                                  parPhase = 0 }
-  in r
-
--- | Interpolate the computation based on the integration time points only.
-interpolate :: Dynamics Double -> Dynamics Double
-interpolate (Dynamics m) = 
-  Dynamics $ \ps -> 
-  if parPhase ps >= 0 then 
-    m ps
-  else 
-    let sc = parSpecs ps
-        t  = parTime ps
-        x  = (t - spcStartTime sc) / spcDT sc
-        n1 = max (floor x) (iterationLoBnd sc)
-        n2 = min (ceiling x) (iterationHiBnd sc)
-        t1 = basicTime sc n1 0
-        t2 = basicTime sc n2 0
-        z1 = m $ ps { parTime = t1, 
-                      parIteration = n1, 
-                      parPhase = 0 }
-        z2 = m $ ps { parTime = t2,
-                      parIteration = n2,
-                      parPhase = 0 }
-        r | t == t1   = z1
-          | t == t2   = z2
-          | otherwise = 
-            do y1 <- z1
-               y2 <- z2
-               return $ y1 + (y2 - y1) * (t - t1) / (t2 - t1)
-    in r
-
--- --
--- -- Memoization
--- --
--- -- The memoization creates such processes, which values are 
--- -- defined and then stored in the cache for the points of
--- -- integration. You should use some kind of interpolation 
--- -- like the interpolate function to process all other time 
--- -- points that don't coincide with the integration points:
--- --
--- --   x = memo interpolate y    -- a linear interpolation
--- --   x = memo discrete y       -- a discrete process
--- --
-
--- | The 'Memo' class specifies a type for which an array can be created.
-class (MArray IOArray e IO) => Memo e where
-  newMemoArray_ :: Ix i => (i, i) -> IO (IOArray i e)
-
--- | The 'UMemo' class specifies a type for which an unboxed array exists.
-class (MArray IOUArray e IO) => UMemo e where
-  newMemoUArray_ :: Ix i => (i, i) -> IO (IOUArray i e)
-
-instance Memo e where
-  newMemoArray_ = newArray_
-    
-instance (MArray IOUArray e IO) => UMemo e where
-  newMemoUArray_ = newArray_
-
--- | Memoize and order the computation in the integration time points using 
--- the specified interpolation and being aware of the Runge-Kutta method.
-memo :: Memo e => (Dynamics e -> Dynamics e) -> Dynamics e 
-       -> Dynamics (Dynamics e)
-memo tr (Dynamics m) = 
-  Dynamics $ \ps ->
-  do let sc = parSpecs ps
-         (phl, phu) = phaseBnds sc
-         (nl, nu)   = iterationBnds sc
-     arr   <- newMemoArray_ ((phl, nl), (phu, nu))
-     nref  <- newIORef 0
-     phref <- newIORef 0
-     let r ps = 
-           do let sc  = parSpecs ps
-                  n   = parIteration ps
-                  ph  = parPhase ps
-                  phu = phaseHiBnd sc 
-                  loop n' ph' = 
-                    if (n' > n) || ((n' == n) && (ph' > ph)) 
-                    then 
-                      readArray arr (ph, n)
-                    else 
-                      let ps' = ps { parIteration = n', parPhase = ph',
-                                     parTime = basicTime sc n' ph' }
-                      in do a <- m ps'
-                            a `seq` writeArray arr (ph', n') a
-                            if ph' >= phu 
-                              then do writeIORef phref 0
-                                      writeIORef nref (n' + 1)
-                                      loop (n' + 1) 0
-                              else do writeIORef phref (ph' + 1)
-                                      loop n' (ph' + 1)
-              n'  <- readIORef nref
-              ph' <- readIORef phref
-              loop n' ph'
-     return $ tr $ Dynamics r
-
--- | Memoize and order the computation in the integration time points using 
--- the specified interpolation and being aware of the Runge-Kutta method.
-umemo :: UMemo e => (Dynamics e -> Dynamics e) -> Dynamics e 
-        -> Dynamics (Dynamics e)
-umemo tr (Dynamics m) = 
-  Dynamics $ \ps ->
-  do let sc = parSpecs ps
-         (phl, phu) = phaseBnds sc
-         (nl, nu)   = iterationBnds sc
-     arr   <- newMemoUArray_ ((phl, nl), (phu, nu))
-     nref  <- newIORef 0
-     phref <- newIORef 0
-     let r ps =
-           do let sc  = parSpecs ps
-                  n   = parIteration ps
-                  ph  = parPhase ps
-                  phu = phaseHiBnd sc 
-                  loop n' ph' = 
-                    if (n' > n) || ((n' == n) && (ph' > ph)) 
-                    then 
-                      readArray arr (ph, n)
-                    else 
-                      let ps' = ps { parIteration = n', 
-                                     parPhase = ph',
-                                     parTime = basicTime sc n' ph' }
-                      in do a <- m ps'
-                            a `seq` writeArray arr (ph', n') a
-                            if ph' >= phu 
-                              then do writeIORef phref 0
-                                      writeIORef nref (n' + 1)
-                                      loop (n' + 1) 0
-                              else do writeIORef phref (ph' + 1)
-                                      loop n' (ph' + 1)
-              n'  <- readIORef nref
-              ph' <- readIORef phref
-              loop n' ph'
-     return $ tr $ Dynamics r
-
--- | Memoize and order the computation in the integration time points using 
--- the specified interpolation and without knowledge of the Runge-Kutta method.
-memo0 :: Memo e => (Dynamics e -> Dynamics e) -> Dynamics e 
-        -> Dynamics (Dynamics e)
-memo0 tr (Dynamics m) = 
-  Dynamics $ \ps ->
-  do let sc   = parSpecs ps
-         bnds = iterationBnds sc
-     arr   <- newMemoArray_ bnds
-     nref  <- newIORef 0
-     let r ps =
-           do let sc = parSpecs ps
-                  n  = parIteration ps
-                  loop n' = 
-                    if n' > n
-                    then 
-                      readArray arr n
-                    else 
-                      let ps' = ps { parIteration = n', parPhase = 0,
-                                     parTime = basicTime sc n' 0 }
-                      in do a <- m ps'
-                            a `seq` writeArray arr n' a
-                            writeIORef nref (n' + 1)
-                            loop (n' + 1)
-              n' <- readIORef nref
-              loop n'
-     return $ tr $ Dynamics r
-
--- | Memoize and order the computation in the integration time points using 
--- the specified interpolation and without knowledge of the Runge-Kutta method.
-umemo0 :: UMemo e => (Dynamics e -> Dynamics e) -> Dynamics e 
-         -> Dynamics (Dynamics e)
-umemo0 tr (Dynamics m) = 
-  Dynamics $ \ps ->
-  do let sc   = parSpecs ps
-         bnds = iterationBnds sc
-     arr   <- newMemoUArray_ bnds
-     nref  <- newIORef 0
-     let r ps =
-           do let sc = parSpecs ps
-                  n  = parIteration ps
-                  loop n' = 
-                    if n' > n
-                    then 
-                      readArray arr n
-                    else 
-                      let ps' = ps { parIteration = n', parPhase = 0,
-                                     parTime = basicTime sc n' 0 }
-                      in do a <- m ps'
-                            a `seq` writeArray arr n' a
-                            writeIORef nref (n' + 1)
-                            loop (n' + 1)
-              n' <- readIORef nref
-              loop n'
-     return $ tr $ Dynamics r
-
---
--- Once
---
-
--- | Call the computation only once.
-once :: Dynamics a -> Dynamics (Dynamics a)
-once (Dynamics m) =
-  Dynamics $ \ps ->
-  do x <- newIORef Nothing
-     let r ps =
-           do a <- readIORef x
-              case a of
-                Just b -> 
-                  return b
-                Nothing ->
-                  do b <- m ps
-                     writeIORef x $ Just b
-                     return $! b
-     return $ Dynamics r
-
---
--- The DynamicsCont Monad
---
--- It looks somewhere like the ContT monad transformer parameterized by 
--- the Dynamics monad, although this analogy is not strong. The main 
--- idea is to represent the continuation as a dynamic process varying 
--- in time.
---
-
-newtype DynamicsCont a = DynamicsCont (Dynamics (a -> IO ()) -> Dynamics ())
-
-instance Monad DynamicsCont where
-  return  = returnDC
-  m >>= k = bindDC m k
-
-returnDC :: a -> DynamicsCont a
-returnDC a = 
-  DynamicsCont $ \(Dynamics c) -> 
-  Dynamics $ \ps -> 
-  do cont' <- c ps
-     cont' a
-                          
-bindDC :: DynamicsCont a -> (a -> DynamicsCont b) -> DynamicsCont b
-bindDC (DynamicsCont m) k =
-  DynamicsCont $ \c ->
-  m $ Dynamics $ \ps -> 
-  let cont' a = let (DynamicsCont m') = k a
-                    (Dynamics u) = m' c
-                in u ps
-  in return cont'
-
-runCont :: DynamicsCont a -> IO (a -> IO ()) -> Dynamics ()
-runCont (DynamicsCont m) f = m $ Dynamics $ const f
-
-instance DynamicsTrans DynamicsCont where
-  liftD (Dynamics m) =
-    DynamicsCont $ \(Dynamics c) ->
-    Dynamics $ \ps ->
-    do cont' <- c ps
-       a <- m ps
-       cont' a
-
-instance Functor DynamicsCont where
-  fmap = liftM
-
-instance MonadIO DynamicsCont where
-  liftIO m =
-    DynamicsCont $ \(Dynamics c) ->
-    Dynamics $ \ps ->
-    do cont' <- c ps
-       a <- m
-       cont' a
-
---
--- The Event Queue (The Dynamics Queue)
---
--- The most exciting thing is that any event is just some value in 
--- the Dynamics monad, i.e. a computation, or saying differently, 
--- a dynamic process that has a single purpose to perform some 
--- side effect. To pass the message, we actually use a closure.
---
-
--- | The 'DynamicsQueue' type represents the event queue.
-data DynamicsQueue = DynamicsQueue { 
-  queuePQ   :: PQ.PriorityQueue (Dynamics (() -> IO ())),
-  queueBusy :: IORef Bool,
-  queueTime :: IORef Double, 
-  queueRun  :: Dynamics () }
-
--- | Create a new event queue.
-newQueue :: Dynamics DynamicsQueue
-newQueue = 
-  Dynamics $ \ps ->
-  do let sc = parSpecs ps
-     f <- newIORef False
-     t <- newIORef $ spcStartTime sc
-     let cont () = return ()
-     pq <- PQ.newQueue $ return cont
-     let q = DynamicsQueue { queuePQ   = pq,
-                             queueBusy = f,
-                             queueTime = t, 
-                             queueRun  = subrunQueue q }
-     return q
-             
--- | Enqueue the event which must be actuated at the specified time.
-enqueueDC :: DynamicsQueue -> Dynamics Double -> Dynamics (() -> IO ()) 
-            -> Dynamics ()
-enqueueDC q (Dynamics t) c = Dynamics r where
-  r ps =
-    do t' <- t ps
-       let pq = queuePQ q
-       PQ.enqueue pq t' c
-    
--- | Enqueue the event which must be actuated at the specified time.
-enqueueD :: DynamicsQueue -> Dynamics Double -> Dynamics () -> Dynamics ()
-enqueueD q t (Dynamics m) = enqueueDC q t (Dynamics c) where
-  c ps = let f () = m ps in return f
-    
-subrunQueue :: DynamicsQueue -> Dynamics ()
-subrunQueue q = Dynamics r where
-  r ps =
-    do let f = queueBusy q
-       f' <- readIORef f
-       unless f' $
-         do writeIORef f True
-            call q ps
-            writeIORef f False
-  call q ps =
-    do let pq = queuePQ q
-       f <- PQ.queueNull pq
-       unless f $
-         do (t2, Dynamics c2) <- PQ.queueFront pq
-            let t = queueTime q
-            t' <- readIORef t
-            when (t2 < t') $ 
-              error "The time value is too small: subrunQueue"
-            when (t2 <= parTime ps) $
-              do writeIORef t t2
-                 PQ.dequeue pq
-                 let sc  = parSpecs ps
-                     t0  = spcStartTime sc
-                     dt  = spcDT sc
-                     n2  = fromInteger $ toInteger $ floor ((t2 - t0) / dt)
-                 k <- c2 $ ps { parTime = t2,
-                               parIteration = n2,
-                               parPhase = -1 }
-                 k ()    -- raise the event
-                 call q ps
-
--- | Run the event queue processing its events.
-runQueue :: DynamicsQueue -> Dynamics ()
-runQueue = queueRun
-
---
--- DynamicsPID and DynamicsProc
---
--- A value in the DynamicsProc monad represents a control process that can be 
--- suspended and resumed at any time. It behaves like a dynamic process too. 
--- Any value in the Dynamics monad can be lifted to the DynamicsProc monad. 
--- Moreover, a value in the DynamicsProc monad can be run in the Dynamics monad.
---
--- A value of the DynamicsPID type is just an identifier of such a process.
---
-
--- Public functions:
---
---   pidQueue
---   holdProcD
---   holdProc
---   passivateProc
---   procPassive
---   reactivateProc
---   runProc
---   procPID
---   newPID
-
--- | Represents a process handler, its PID.
-data DynamicsPID = 
-  DynamicsPID { pidQueue   :: DynamicsQueue,  -- ^ Return the bound event queue.
-                pidStarted :: IORef Bool,
-                pidCont    :: IORef (Maybe (Dynamics (() -> IO ()))) }
-
--- | Specifies a discontinuous process that can be suspended at any time
--- and then resumed later.
-newtype DynamicsProc a = DynamicsProc (DynamicsPID -> DynamicsCont a)
-
--- | Hold the process for the specified time period.
-holdProcD :: Dynamics Double -> DynamicsProc ()
-holdProcD t =
-  DynamicsProc $ \pid ->
-  DynamicsCont $ \c ->
-  enqueueDC (pidQueue pid) (t + time) c
-
--- | Hold the process for the specified time period.
-holdProc :: Double -> DynamicsProc ()
-holdProc t = holdProcD $ return t
-
--- | Passivate the process.
-passivateProc :: DynamicsProc ()
-passivateProc =
-  DynamicsProc $ \pid ->
-  DynamicsCont $ \c ->
-  Dynamics $ \ps ->
-  do let x = pidCont pid
-     a <- readIORef x
-     case a of
-       Nothing -> writeIORef x $ Just c
-       Just _  -> error "Cannot passivate the process twice: passivate"
-
--- | Test whether the process with the specified PID is passivated.
-procPassive :: DynamicsPID -> DynamicsProc Bool
-procPassive pid =
-  DynamicsProc $ \_ ->
-  DynamicsCont $ \(Dynamics c) ->
-  Dynamics $ \ps ->
-  do cont' <- c ps 
-     let x = pidCont pid
-     a <- readIORef x
-     case a of
-       Nothing -> cont' False
-       Just _  -> cont' True
-
--- | Reactivate a process with the specified PID.
-reactivateProc :: DynamicsPID -> DynamicsProc ()
-reactivateProc pid =
-  DynamicsProc $ \pid' ->
-  DynamicsCont $ \c@(Dynamics cont) ->
-  Dynamics $ \ps ->
-  do let x = pidCont pid
-     a <- readIORef x
-     case a of
-       Nothing ->
-         do cont' <- cont ps
-            cont' ()
-       Just (Dynamics cont2) ->
-         do writeIORef x Nothing
-            let Dynamics m = enqueueDC (pidQueue pid') time c
-            m ps
-            cont2' <- cont2 ps
-            cont2' ()
-
--- | Start the process with the specified PID at the desired time.
-runProc :: DynamicsProc () -> DynamicsPID -> Dynamics Double -> Dynamics ()
-runProc (DynamicsProc p) pid t =
-  runCont m r
-    where m = do y <- liftIO $ readIORef (pidStarted pid)
-                 if y 
-                   then error $
-                        "A process with such PID " ++
-                        "has been started already: runProc"
-                   else liftIO $ writeIORef (pidStarted pid) True
-                 DynamicsCont $ \c -> enqueueDC (pidQueue pid) t c
-                 p pid
-          r = let f () = return () in return f
-
--- | Return the current process PID.
-procPID :: DynamicsProc DynamicsPID
-procPID = DynamicsProc $ \pid -> return pid
-
--- | Create a new process PID.
-newPID :: DynamicsQueue -> Dynamics DynamicsPID
-newPID q =
-  do x <- liftIO $ newIORef Nothing
-     y <- liftIO $ newIORef False
-     return DynamicsPID { pidQueue   = q,
-                          pidStarted = y,
-                          pidCont    = x }
-
-instance Eq DynamicsPID where
-  x == y = pidCont x == pidCont y    -- for the references are unique
-
-instance Monad DynamicsProc where
-  return  = returnDP
-  m >>= k = bindDP m k
-
-returnDP :: a -> DynamicsProc a
-returnDP a = DynamicsProc (\pid -> return a)
-
-bindDP :: DynamicsProc a -> (a -> DynamicsProc b) -> DynamicsProc b
-bindDP (DynamicsProc m) k = 
-  DynamicsProc $ \pid -> 
-  do a <- m pid
-     let DynamicsProc m' = k a
-     m' pid
-
-instance Functor DynamicsProc where
-  fmap = liftM
-
-instance DynamicsTrans DynamicsProc where
-  liftD m = DynamicsProc $ \pid -> liftD m
-  
-instance MonadIO DynamicsProc where
-  liftIO m = DynamicsProc $ \pid -> liftIO m
-
---
--- DynamicsResource
---
-  
--- Public functions:  
---
---   resourceQueue
---   resourceInitCount
---   resourceCount
---   requestResource
---   releaseResource
---   newResource
-  
--- | Represents a limited resource.
-data DynamicsResource = 
-  DynamicsResource { resourceQueue     :: DynamicsQueue,  
-                     -- ^ Return the bound event queue.
-                     resourceInitCount :: Int,
-                     -- ^ Return the initial count of the resource.
-                     resourceCountRef  :: IORef Int, 
-                     resourceWaitQueue :: Q.Queue (Dynamics (() -> IO ()))}
-
-instance Eq DynamicsResource where
-  x == y = resourceCountRef x == resourceCountRef y  -- unique references
-
--- | Create a new resource with the specified initial count.
-newResource :: DynamicsQueue -> Int -> Dynamics DynamicsResource
-newResource q initCount =
-  Dynamics $ \ps ->
-  do countRef  <- newIORef initCount
-     waitQueue <- Q.newQueue
-     return DynamicsResource { resourceQueue     = q,
-                               resourceInitCount = initCount,
-                               resourceCountRef  = countRef,
-                               resourceWaitQueue = waitQueue }
-
--- | Return the current count of the resource.
-resourceCount :: DynamicsResource -> DynamicsProc Int
-resourceCount r =
-  DynamicsProc $ \_ ->
-  DynamicsCont $ \(Dynamics c) ->
-  Dynamics $ \ps ->
-  do cont' <- c ps 
-     a <- readIORef (resourceCountRef r)
-     cont' a
-
--- | Request for the resource decreasing its count in case of success,
--- otherwise suspending the discontinuous process until some other 
--- process releases the resource.
-requestResource :: DynamicsResource -> DynamicsProc ()
-requestResource r =
-  DynamicsProc $ \_ ->
-  DynamicsCont $ \c@(Dynamics cont) ->
-  Dynamics $ \ps ->
-  do a <- readIORef (resourceCountRef r)
-     if a == 0 
-       then Q.enqueue (resourceWaitQueue r) c
-       else do let a' = a - 1
-               a' `seq` writeIORef (resourceCountRef r) a'
-               cont' <- cont ps
-               cont' ()
-
--- | Release the resource increasing its count and resuming one of the
--- previously suspended processes as possible.
-releaseResource :: DynamicsResource -> DynamicsProc ()
-releaseResource r =
-  DynamicsProc $ \_ ->
-  DynamicsCont $ \(Dynamics c) ->
-  Dynamics $ \ps ->
-  do a <- readIORef (resourceCountRef r)
-     let a' = a + 1
-     when (a' > resourceInitCount r) $
-       error $
-       "The resource count cannot be greater than " ++
-       "its initial value: releaseResource."
-     f <- Q.queueNull (resourceWaitQueue r)
-     if f 
-       then a' `seq` writeIORef (resourceCountRef r) a'
-       else do c2 <- Q.queueFront (resourceWaitQueue r)
-               Q.dequeue (resourceWaitQueue r)
-               let Dynamics m = enqueueDC (resourceQueue r) time c2
-               m ps
-     cont' <- c ps
-     cont' ()
-
---
--- DynamicsRef
---
-
--- | The 'DynamicsRef' type represents a mutable variable similar to 
--- the 'IORef' variable but only bound to some event queue, which makes 
--- the variable coordinated with that queue.
-data DynamicsRef a = 
-  DynamicsRef { refQueue  :: DynamicsQueue,    -- ^ Return the bound event queue.
-                refRunner :: Dynamics (),
-                refValue  :: IORef a }
-
--- | Create a new reference bound to the specified event queue.
-newRef :: DynamicsQueue -> a -> Dynamics (DynamicsRef a)
-newRef q a =
-  do x <- liftIO $ newIORef a
-     return DynamicsRef { refQueue  = q,
-                          refRunner = runQueue q,
-                          refValue  = x }
-     
--- | Read the value of a reference, forcing the bound event queue to raise 
--- the events in case of need.
-readRef :: DynamicsRef a -> Dynamics a
-readRef r = Dynamics $ \ps -> 
-  do let Dynamics m = refRunner r
-     m ps
-     readIORef (refValue r)
-
--- | Write a new value into the reference.
-writeRef :: DynamicsRef a -> a -> Dynamics ()
-writeRef r a = Dynamics $ \ps -> 
-  do writeIORef (refValue r) a
-     let Dynamics m = refRunner r 
-     m ps
-
--- | Mutate the contents of the reference, forcing the bound event queue to
--- raise all pending events in case of need.
-modifyRef :: DynamicsRef a -> (a -> a) -> Dynamics ()
-modifyRef r f = Dynamics $ \ps -> 
-  do let Dynamics m = refRunner r 
-     m ps
-     modifyIORef (refValue r) f
-
--- | A strict version of the 'writeRef' function.
-writeRef' :: DynamicsRef a -> a -> Dynamics ()
-writeRef' r a = a `seq` writeRef r a
-
--- | A strict version of the 'modifyRef' function.
-modifyRef' :: DynamicsRef a -> (a -> a) -> Dynamics ()
-modifyRef' r f = Dynamics $ \ps ->
-  do let Dynamics m = refRunner r
-     m ps
-     a <- readIORef (refValue r)
-     let b = f a
-     b `seq` writeIORef (refValue r) b
-
---
--- Agent-based Modeling
---
-
--- Public functions:
---
---   agentQueue
---   agentState
---   activateState
---   initState
---   stateAgent     
---   stateParent
---   addTimeoutD
---   addTimeout
---   addTimerD
---   addTimer
---   stateActivation
---   stateDeactivation
---   newState     
---   newSubstate
---   newAgent
-
--- | Represents an agent.
-data Agent = Agent { agentQueue :: DynamicsQueue,
-                     -- ^ Return the bound event queue.
-                     agentModeRef :: IORef AgentMode,
-                     agentStateRef :: IORef (Maybe AgentState) }
-
--- | Represents the agent state.
-data AgentState = AgentState { stateAgent :: Agent,
-                               -- ^ Return the corresponded agent.
-                               stateParent :: Maybe AgentState,
-                               -- ^ Return the parent state or 'Nothing'.
-                               stateActivateRef :: IORef (Dynamics ()),
-                               stateDeactivateRef :: IORef (Dynamics ()), 
-                               stateVersionRef :: IORef Int }
-                  
-data AgentMode = CreationMode
-               | InitialMode
-               | TransientMode
-               | ProcessingMode
-                      
-instance Eq Agent where
-  x == y = agentStateRef x == agentStateRef y      -- unique references
-  
-instance Eq AgentState where
-  x == y = stateVersionRef x == stateVersionRef y  -- unique references
-
-findPath :: AgentState -> AgentState -> ([AgentState], [AgentState])
-findPath source target = 
-  if stateAgent source == stateAgent target 
-  then
-    partitionPath path1 path2
-  else
-    error "Different agents: findPath."
-      where
-        path1 = fullPath source []
-        path2 = fullPath target []
-        fullPath st acc =
-          case stateParent st of
-            Nothing  -> st : acc
-            Just st' -> fullPath st' (st : acc)
-        partitionPath path1 path2 =
-          case (path1, path2) of
-            (h1 : t1, [h2]) | h1 == h2 -> 
-              (reverse path1, path2)
-            (h1 : t1, h2 : t2) | h1 == h2 -> 
-              partitionPath t1 t2
-            _ -> 
-              (reverse path1, path2)
-            
-traversePath :: AgentState -> AgentState -> Dynamics ()
-traversePath source target =
-  let (path1, path2) = findPath source target
-      agent = stateAgent source
-      activate st ps   =
-        do Dynamics m <- readIORef (stateActivateRef st)
-           m ps
-      deactivate st ps =
-        do Dynamics m <- readIORef (stateDeactivateRef st)
-           m ps
-  in Dynamics $ \ps ->
-       do writeIORef (agentModeRef agent) TransientMode
-          forM_ path1 $ \st ->
-            do writeIORef (agentStateRef agent) (Just st)
-               deactivate st ps
-               -- it makes all timeout and timer handlers obsolete
-               modifyIORef (stateVersionRef st) (1 +)
-          forM_ path2 $ \st ->
-            do when (st == target) $
-                 writeIORef (agentModeRef agent) InitialMode
-               writeIORef (agentStateRef agent) (Just st)
-               activate st ps
-               when (st == target) $
-                 writeIORef (agentModeRef agent) ProcessingMode
-
--- | Add to the state a timeout handler that will be actuated 
--- in the specified time period, while the state remains active.
-addTimeoutD :: AgentState -> Dynamics Double -> Dynamics () -> Dynamics ()
-addTimeoutD st t (Dynamics action) =
-  Dynamics $ \ps ->
-  do v <- readIORef (stateVersionRef st)
-     let m1 = Dynamics $ \ps ->
-           do v' <- readIORef (stateVersionRef st)
-              when (v == v') $ action ps
-         q = agentQueue (stateAgent st)
-         Dynamics m2 = enqueueD q (t + time) m1
-     m2 ps
-
--- | Add to the state a timer handler that will be actuated
--- in the specified time period and then repeated again many times,
--- while the state remains active.
-addTimerD :: AgentState -> Dynamics Double -> Dynamics () -> Dynamics ()
-addTimerD st t (Dynamics action) =
-  Dynamics $ \ps ->
-  do v <- readIORef (stateVersionRef st)
-     let m1 = Dynamics $ \ps ->
-           do v' <- readIORef (stateVersionRef st)
-              when (v == v') $ do { m2 ps; action ps }
-         q = agentQueue (stateAgent st)
-         Dynamics m2 = enqueueD q (t + time) m1
-     m2 ps
-
--- | Add to the state a timeout handler that will be actuated 
--- in the specified time period, while the state remains active.
-addTimeout :: AgentState -> Double -> Dynamics () -> Dynamics ()
-addTimeout st t = addTimeoutD st (return t)
-
--- | Add to the state a timer handler that will be actuated
--- in the specified time period and then repeated again many times,
--- while the state remains active.
-addTimer :: AgentState -> Double -> Dynamics () -> Dynamics ()
-addTimer st t = addTimerD st (return t)
-
--- | Create a new state.
-newState :: Agent -> Dynamics AgentState
-newState agent =
-  Dynamics $ \ps ->
-  do aref <- newIORef $ return ()
-     dref <- newIORef $ return ()
-     vref <- newIORef 0
-     return AgentState { stateAgent = agent,
-                         stateParent = Nothing,
-                         stateActivateRef = aref,
-                         stateDeactivateRef = dref,
-                         stateVersionRef = vref }
-
--- | Create a child state.
-newSubstate :: AgentState -> Dynamics AgentState
-newSubstate parent =
-  Dynamics $ \ps ->
-  do let agent = stateAgent parent 
-     aref <- newIORef $ return ()
-     dref <- newIORef $ return ()
-     vref <- newIORef 0
-     return AgentState { stateAgent = agent,
-                         stateParent = Just parent,
-                         stateActivateRef= aref,
-                         stateDeactivateRef = dref,
-                         stateVersionRef = vref }
-
--- | Create an agent bound with the specified event queue.
-newAgent :: DynamicsQueue -> Dynamics Agent
-newAgent queue =
-  Dynamics $ \ps ->
-  do modeRef    <- newIORef CreationMode
-     stateRef   <- newIORef Nothing
-     return Agent { agentQueue = queue,
-                    agentModeRef = modeRef,
-                    agentStateRef = stateRef }
-
--- | Return the selected downmost active state.
-agentState :: Agent -> Dynamics (Maybe AgentState)
-agentState agent =
-  Dynamics $ \ps -> 
-  do let Dynamics m = queueRun $ agentQueue agent 
-     m ps    -- ensure that the agent state is actual
-     readIORef (agentStateRef agent)
-                   
--- | Select the next downmost active state.       
-activateState :: AgentState -> Dynamics ()
-activateState st =
-  Dynamics $ \ps ->
-  do let agent = stateAgent st
-         Dynamics m = queueRun $ agentQueue agent 
-     m ps    -- ensure that the agent state is actual
-     mode <- readIORef (agentModeRef agent)
-     case mode of
-       CreationMode ->
-         case stateParent st of
-           Just _ ->
-             error $ 
-             "To run the agent for the first time, an initial state " ++
-             "must be top-level: activateState."
-           Nothing ->
-             do writeIORef (agentModeRef agent) InitialMode
-                writeIORef (agentStateRef agent) (Just st)
-                Dynamics m <- readIORef (stateActivateRef st)
-                m ps
-                writeIORef (agentModeRef agent) ProcessingMode
-       InitialMode ->
-         error $ 
-         "Use the initState function during " ++
-         "the state activation: activateState."
-       TransientMode ->
-         error $
-         "Use the initState function during " ++
-         "the state activation: activateState."
-       ProcessingMode ->
-         do Just st0 <- readIORef (agentStateRef agent)
-            let Dynamics m = traversePath st0 st
-            m ps
-              
--- | Activate the child state during the direct activation of 
--- the parent state. This call is ignored in other cases.
-initState :: AgentState -> Dynamics ()
-initState st =
-  Dynamics $ \ps ->
-  do let agent = stateAgent st
-         Dynamics m = queueRun $ agentQueue agent 
-     m ps    -- ensure that the agent state is actual
-     mode <- readIORef (agentModeRef agent)
-     case mode of
-       CreationMode ->
-         error $
-         "To run the agent for the fist time, use " ++
-         "the activateState function: initState."
-       InitialMode ->
-         do Just st0 <- readIORef (agentStateRef agent)
-            let Dynamics m = traversePath st0 st
-            m ps
-       TransientMode -> 
-         return ()
-       ProcessingMode ->
-         error $
-         "Use the activateState function everywhere outside " ++
-         "the state activation: initState."
-
--- | Set the activation computation for the specified state.
-stateActivation :: AgentState -> Dynamics () -> Dynamics ()
-stateActivation st action =
-  Dynamics $ \ps ->
-  writeIORef (stateActivateRef st) action
-  
--- | Set the deactivation computation for the specified state.
-stateDeactivation :: AgentState -> Dynamics () -> Dynamics ()
-stateDeactivation st action =
-  Dynamics $ \ps ->
-  writeIORef (stateDeactivateRef st) action
-  
-          
+-- |
+-- Module     : Simulation.Aivika.Dynamics
+-- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.0.3
+--
+-- The module defines the 'Dynamics' monad representing an abstract dynamic 
+-- process, i.e. a time varying polymorphic function. 
+-- 
+-- This is a key point of the Aivika simulation library. With help of this monad 
+-- we can simulate the system of ordinary differential equations (ODEs) of 
+-- System Dynamics, define the tasks of Discrete Event Simulation (DES) supporting 
+-- different paradigms. Also we can use the Agent-based Modeling. Thus, 
+-- we can create hybrid simulation models.
+--
+module Simulation.Aivika.Dynamics 
+       (Dynamics,
+        Specs(..),
+        Method(..),
+        runDynamics1,
+        runDynamics1_,
+        runDynamics,
+        runDynamics_,
+        runDynamicsIO,
+        runDynamicsSeries1,
+        runDynamicsSeries1_,
+        runDynamicsSeries,
+        runDynamicsSeries_,
+        printDynamics1,
+        printDynamics) where
+
+import Simulation.Aivika.Dynamics.Internal.Dynamics
diff --git a/Simulation/Aivika/Dynamics/Agent.hs b/Simulation/Aivika/Dynamics/Agent.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Dynamics/Agent.hs
@@ -0,0 +1,260 @@
+
+-- |
+-- Module     : Simulation.Aivika.Dynamics.Agent
+-- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.0.3
+--
+-- This module introduces an agent-based modeling.
+--
+
+module Simulation.Aivika.Dynamics.Agent
+       (Agent,
+        AgentState,
+        newAgent,
+        newState,
+        newSubstate,
+        agentQueue,
+        agentState,
+        activateState,
+        initState,
+        stateAgent,
+        stateParent,
+        addTimeout,
+        addTimer,
+        stateActivation,
+        stateDeactivation) where
+
+import Data.IORef
+import Control.Monad
+
+import Simulation.Aivika.Dynamics.Internal.Dynamics
+import Simulation.Aivika.Dynamics.EventQueue
+
+--
+-- Agent-based Modeling
+--
+
+-- | Represents an agent.
+data Agent = Agent { agentQueue :: EventQueue,
+                     -- ^ Return the bound event queue.
+                     agentModeRef :: IORef AgentMode,
+                     agentStateRef :: IORef (Maybe AgentState) }
+
+-- | Represents the agent state.
+data AgentState = AgentState { stateAgent :: Agent,
+                               -- ^ Return the corresponded agent.
+                               stateParent :: Maybe AgentState,
+                               -- ^ Return the parent state or 'Nothing'.
+                               stateActivateRef :: IORef (Dynamics ()),
+                               stateDeactivateRef :: IORef (Dynamics ()), 
+                               stateVersionRef :: IORef Int }
+                  
+data AgentMode = CreationMode
+               | InitialMode
+               | TransientMode
+               | ProcessingMode
+                      
+instance Eq Agent where
+  x == y = agentStateRef x == agentStateRef y      -- unique references
+  
+instance Eq AgentState where
+  x == y = stateVersionRef x == stateVersionRef y  -- unique references
+
+findPath :: AgentState -> AgentState -> ([AgentState], [AgentState])
+findPath source target = 
+  if stateAgent source == stateAgent target 
+  then
+    partitionPath path1 path2
+  else
+    error "Different agents: findPath."
+      where
+        path1 = fullPath source []
+        path2 = fullPath target []
+        fullPath st acc =
+          case stateParent st of
+            Nothing  -> st : acc
+            Just st' -> fullPath st' (st : acc)
+        partitionPath path1 path2 =
+          case (path1, path2) of
+            (h1 : t1, [h2]) | h1 == h2 -> 
+              (reverse path1, path2)
+            (h1 : t1, h2 : t2) | h1 == h2 -> 
+              partitionPath t1 t2
+            _ -> 
+              (reverse path1, path2)
+            
+traversePath :: AgentState -> AgentState -> Dynamics ()
+traversePath source target =
+  let (path1, path2) = findPath source target
+      agent = stateAgent source
+      activate st p =
+        do Dynamics m <- readIORef (stateActivateRef st)
+           m p
+      deactivate st p =
+        do Dynamics m <- readIORef (stateDeactivateRef st)
+           m p
+  in Dynamics $ \p ->
+       do writeIORef (agentModeRef agent) TransientMode
+          forM_ path1 $ \st ->
+            do writeIORef (agentStateRef agent) (Just st)
+               deactivate st p
+               -- it makes all timeout and timer handlers obsolete
+               modifyIORef (stateVersionRef st) (1 +)
+          forM_ path2 $ \st ->
+            do when (st == target) $
+                 writeIORef (agentModeRef agent) InitialMode
+               writeIORef (agentStateRef agent) (Just st)
+               activate st p
+               when (st == target) $
+                 writeIORef (agentModeRef agent) ProcessingMode
+
+-- | Add to the state a timeout handler that will be actuated 
+-- in the specified time period, while the state remains active.
+addTimeout :: AgentState -> Double -> Dynamics () -> Dynamics ()
+addTimeout st dt (Dynamics action) =
+  Dynamics $ \p ->
+  do v <- readIORef (stateVersionRef st)
+     let m1 = Dynamics $ \p ->
+           do v' <- readIORef (stateVersionRef st)
+              when (v == v') $ action p
+         q = agentQueue (stateAgent st)
+         Dynamics m2 = enqueue q (pointTime p + dt) m1
+     m2 p
+
+-- | Add to the state a timer handler that will be actuated
+-- in the specified time period and then repeated again many times,
+-- while the state remains active.
+addTimer :: AgentState -> Dynamics Double -> Dynamics () -> Dynamics ()
+addTimer st (Dynamics dt) (Dynamics action) =
+  Dynamics $ \p ->
+  do v <- readIORef (stateVersionRef st)
+     let m1 = Dynamics $ \p ->
+           do v' <- readIORef (stateVersionRef st)
+              when (v == v') $ do { m2 p; action p }
+         q = agentQueue (stateAgent st)
+         Dynamics m2 = 
+           Dynamics $ \p ->
+           do dt' <- dt p
+              let Dynamics m3 = enqueue q (pointTime p + dt') m1
+              m3 p
+     m2 p
+
+-- | Create a new state.
+newState :: Agent -> Dynamics AgentState
+newState agent =
+  Dynamics $ \p ->
+  do aref <- newIORef $ return ()
+     dref <- newIORef $ return ()
+     vref <- newIORef 0
+     return AgentState { stateAgent = agent,
+                         stateParent = Nothing,
+                         stateActivateRef = aref,
+                         stateDeactivateRef = dref,
+                         stateVersionRef = vref }
+
+-- | Create a child state.
+newSubstate :: AgentState -> Dynamics AgentState
+newSubstate parent =
+  Dynamics $ \p ->
+  do let agent = stateAgent parent 
+     aref <- newIORef $ return ()
+     dref <- newIORef $ return ()
+     vref <- newIORef 0
+     return AgentState { stateAgent = agent,
+                         stateParent = Just parent,
+                         stateActivateRef= aref,
+                         stateDeactivateRef = dref,
+                         stateVersionRef = vref }
+
+-- | Create an agent bound with the specified event queue.
+newAgent :: EventQueue -> Dynamics Agent
+newAgent queue =
+  Dynamics $ \p ->
+  do modeRef    <- newIORef CreationMode
+     stateRef   <- newIORef Nothing
+     return Agent { agentQueue = queue,
+                    agentModeRef = modeRef,
+                    agentStateRef = stateRef }
+
+-- | Return the selected downmost active state.
+agentState :: Agent -> Dynamics (Maybe AgentState)
+agentState agent =
+  Dynamics $ \p -> 
+  do let Dynamics m = queueRun $ agentQueue agent 
+     m p    -- ensure that the agent state is actual
+     readIORef (agentStateRef agent)
+                   
+-- | Select the next downmost active state.       
+activateState :: AgentState -> Dynamics ()
+activateState st =
+  Dynamics $ \p ->
+  do let agent = stateAgent st
+         Dynamics m = queueRun $ agentQueue agent 
+     m p    -- ensure that the agent state is actual
+     mode <- readIORef (agentModeRef agent)
+     case mode of
+       CreationMode ->
+         case stateParent st of
+           Just _ ->
+             error $ 
+             "To run the agent for the first time, an initial state " ++
+             "must be top-level: activateState."
+           Nothing ->
+             do writeIORef (agentModeRef agent) InitialMode
+                writeIORef (agentStateRef agent) (Just st)
+                Dynamics m <- readIORef (stateActivateRef st)
+                m p
+                writeIORef (agentModeRef agent) ProcessingMode
+       InitialMode ->
+         error $ 
+         "Use the initState function during " ++
+         "the state activation: activateState."
+       TransientMode ->
+         error $
+         "Use the initState function during " ++
+         "the state activation: activateState."
+       ProcessingMode ->
+         do Just st0 <- readIORef (agentStateRef agent)
+            let Dynamics m = traversePath st0 st
+            m p
+              
+-- | Activate the child state during the direct activation of 
+-- the parent state. This call is ignored in other cases.
+initState :: AgentState -> Dynamics ()
+initState st =
+  Dynamics $ \p ->
+  do let agent = stateAgent st
+         Dynamics m = queueRun $ agentQueue agent 
+     m p    -- ensure that the agent state is actual
+     mode <- readIORef (agentModeRef agent)
+     case mode of
+       CreationMode ->
+         error $
+         "To run the agent for the fist time, use " ++
+         "the activateState function: initState."
+       InitialMode ->
+         do Just st0 <- readIORef (agentStateRef agent)
+            let Dynamics m = traversePath st0 st
+            m p
+       TransientMode -> 
+         return ()
+       ProcessingMode ->
+         error $
+         "Use the activateState function everywhere outside " ++
+         "the state activation: initState."
+
+-- | Set the activation computation for the specified state.
+stateActivation :: AgentState -> Dynamics () -> Dynamics ()
+stateActivation st action =
+  Dynamics $ \p ->
+  writeIORef (stateActivateRef st) action
+  
+-- | Set the deactivation computation for the specified state.
+stateDeactivation :: AgentState -> Dynamics () -> Dynamics ()
+stateDeactivation st action =
+  Dynamics $ \p ->
+  writeIORef (stateDeactivateRef st) action
+  
diff --git a/Simulation/Aivika/Dynamics/Base.hs b/Simulation/Aivika/Dynamics/Base.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Dynamics/Base.hs
@@ -0,0 +1,40 @@
+
+-- |
+-- Module     : Simulation.Aivika.Dynamics.Base
+-- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.0.3
+--
+-- This module defines basic functions for the 'Dynamics' monad.
+--
+
+module Simulation.Aivika.Dynamics.Base
+       (-- * Time Parameters
+        starttime,
+        stoptime,
+        dt,
+        time,
+        -- * Interpolation and Initial Value
+        initD,
+        discrete,
+        interpolate,
+        -- * Memoization
+        memo,
+        umemo,
+        memo0,
+        umemo0,
+        -- * Iterating
+        iterateD,
+        -- * Fold
+        foldD1,
+        foldD,
+        -- * Norming
+        divideD) where
+
+import Simulation.Aivika.Dynamics.Internal.Dynamics
+import Simulation.Aivika.Dynamics.Internal.Time
+import Simulation.Aivika.Dynamics.Internal.Interpolate
+import Simulation.Aivika.Dynamics.Internal.Memo
+import Simulation.Aivika.Dynamics.Internal.Fold
diff --git a/Simulation/Aivika/Dynamics/Cont.hs b/Simulation/Aivika/Dynamics/Cont.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Dynamics/Cont.hs
@@ -0,0 +1,20 @@
+
+-- |
+-- Module     : Simulation.Aivika.Dynamics.Cont
+-- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.0.3
+--
+-- The 'Cont' monad looks somewhere like the standard ContT monad transformer 
+-- parameterized by the 'Dynamics' monad, although this analogy is not strong. 
+-- The main idea is to represent the continuation as a dynamic process varying 
+-- in time.
+--
+module Simulation.Aivika.Dynamics.Cont
+       (Cont,
+        runCont) where
+
+import Simulation.Aivika.Dynamics.Internal.Dynamics
+import Simulation.Aivika.Dynamics.Internal.Cont
diff --git a/Simulation/Aivika/Dynamics/EventQueue.hs b/Simulation/Aivika/Dynamics/EventQueue.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Dynamics/EventQueue.hs
@@ -0,0 +1,89 @@
+
+-- |
+-- Module     : Simulation.Aivika.Dynamics.EventQueue
+-- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.0.3
+--
+-- The module introduces the event queue. Any event is the Dynamics computation,
+-- or, saying differently, a dynamic process that has a single purpose 
+-- to perform some side effect at the desired time. To pass the message, 
+-- we actually use a closure.
+--
+module Simulation.Aivika.Dynamics.EventQueue
+       (EventQueue,
+        newQueue,
+        enqueueCont,
+        enqueue,
+        queueRun) where
+
+import Data.IORef
+import Control.Monad
+
+import Simulation.Aivika.Dynamics.Internal.Dynamics
+import qualified Simulation.Aivika.PriorityQueue as PQ
+
+-- | The 'EventQueue' type represents the event queue.
+data EventQueue = EventQueue { 
+  queuePQ   :: PQ.PriorityQueue (Dynamics (() -> IO ())),
+  queueRun  :: Dynamics (),   -- ^ Run the event queue processing its events
+  queueBusy :: IORef Bool,
+  queueTime :: IORef Double }
+
+-- | Create a new event queue.
+newQueue :: Dynamics EventQueue
+newQueue = 
+  Dynamics $ \p ->
+  do let sc = pointSpecs p
+     f <- newIORef False
+     t <- newIORef $ spcStartTime sc
+     pq <- PQ.newQueue
+     let q = EventQueue { queuePQ   = pq,
+                          queueRun  = runQueue q,
+                          queueBusy = f,
+                          queueTime = t }
+     return q
+             
+-- | Enqueue the event which must be actuated at the specified time.
+enqueueCont :: EventQueue -> Double -> Dynamics (() -> IO ()) -> Dynamics ()
+enqueueCont q t c = Dynamics r where
+  r p = let pq = queuePQ q in PQ.enqueue pq t c
+    
+-- | Enqueue the event which must be actuated at the specified time.
+enqueue :: EventQueue -> Double -> Dynamics () -> Dynamics ()
+enqueue q t (Dynamics m) = enqueueCont q t (Dynamics c) where
+  c p = let f () = m p in return f
+    
+-- | Run the event queue processing its events.
+runQueue :: EventQueue -> Dynamics ()
+runQueue q = Dynamics r where
+  r p =
+    do let f = queueBusy q
+       f' <- readIORef f
+       unless f' $
+         do writeIORef f True
+            call q p
+            writeIORef f False
+  call q p =
+    do let pq = queuePQ q
+       f <- PQ.queueNull pq
+       unless f $
+         do (t2, Dynamics c2) <- PQ.queueFront pq
+            let t = queueTime q
+            t' <- readIORef t
+            when (t2 < t') $ 
+              error "The time value is too small: subrunQueue"
+            when (t2 <= pointTime p) $
+              do writeIORef t t2
+                 PQ.dequeue pq
+                 let sc  = pointSpecs p
+                     t0  = spcStartTime sc
+                     dt  = spcDT sc
+                     n2  = fromInteger $ toInteger $ floor ((t2 - t0) / dt)
+                 k <- c2 $ p { pointTime = t2,
+                               pointIteration = n2,
+                               pointPhase = -1 }
+                 k ()    -- raise the event
+                 call q p
diff --git a/Simulation/Aivika/Dynamics/Internal/Cont.hs b/Simulation/Aivika/Dynamics/Internal/Cont.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Dynamics/Internal/Cont.hs
@@ -0,0 +1,84 @@
+
+-- |
+-- Module     : Simulation.Aivika.Dynamics.Internal.Cont
+-- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.0.3
+--
+-- The 'Cont' monad looks somewhere like the standard ContT monad transformer 
+-- parameterized by the 'Dynamics' monad, although this analogy is not strong. 
+-- The main idea is to represent the continuation as a dynamic process varying 
+-- in time.
+--
+module Simulation.Aivika.Dynamics.Internal.Cont
+       (Cont(..),
+        runCont) where
+
+import Control.Monad
+import Control.Monad.Trans
+
+import Simulation.Aivika.Dynamics.Internal.Dynamics
+import Simulation.Aivika.Dynamics.Lift
+
+-- | The 'Cont' type is similar to the standard Cont monad but only
+-- the continuation is represented as a dynamic process varying in time.
+newtype Cont a = Cont (Dynamics (a -> IO ()) -> Dynamics ())
+
+instance Monad Cont where
+  return  = returnC
+  m >>= k = bindC m k
+
+instance Lift Cont where
+  liftD = liftC
+
+instance Functor Cont where
+  fmap = liftM
+
+instance MonadIO Cont where
+  liftIO = liftIOC 
+
+returnC :: a -> Cont a
+{-# INLINE returnC #-}
+returnC a = 
+  Cont $ \(Dynamics c) -> 
+  Dynamics $ \p -> 
+  do cont' <- c p
+     cont' a
+                          
+bindC :: Cont a -> (a -> Cont b) -> Cont b
+{-# INLINE bindC #-}
+bindC (Cont m) k =
+  Cont $ \c ->
+  m $ Dynamics $ \p -> 
+  let cont' a = let (Cont m') = k a
+                    (Dynamics u) = m' c
+                in u p
+  in return cont'
+
+-- | Run the 'Cont' computation.
+runCont :: Cont a -> IO (a -> IO ()) -> Dynamics ()
+{-# INLINE runCont #-}
+runCont (Cont m) f = m $ Dynamics $ const f
+
+-- | Lift the 'Dynamics' computation.
+liftC :: Dynamics a -> Cont a
+{-# INLINE liftC #-}
+liftC (Dynamics m) =
+  Cont $ \(Dynamics c) ->
+  Dynamics $ \p ->
+  do cont' <- c p
+     a <- m p
+     cont' a
+     
+-- | Lift the IO computation.
+liftIOC :: IO a -> Cont a
+{-# INLINE liftIOC #-}
+liftIOC m =
+  Cont $ \(Dynamics c) ->
+  Dynamics $ \p ->
+  do cont' <- c p
+     a <- m
+     cont' a
+  
diff --git a/Simulation/Aivika/Dynamics/Internal/Dynamics.hs b/Simulation/Aivika/Dynamics/Internal/Dynamics.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Dynamics/Internal/Dynamics.hs
@@ -0,0 +1,389 @@
+
+-- |
+-- Module     : Simulation.Aivika.Dynamics.Internal.Dynamics
+-- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.0.3
+--
+-- The module defines the 'Dynamics' monad representing an abstract dynamic 
+-- process, i.e. a time varying polymorphic function. 
+-- 
+-- This is a key point of the Aivika simulation library. With help of this monad 
+-- we can simulate the system of ordinary differential equations (ODEs) of 
+-- System Dynamics, define the tasks of Discrete Event Simulation (DES) supporting 
+-- different paradigms. Also we can use the Agent-based Modeling. Thus, 
+-- we can create hybrid simulation models.
+--
+module Simulation.Aivika.Dynamics.Internal.Dynamics
+       (-- * Dynamics
+        Dynamics(..),
+        Point(..),
+        Specs(..),
+        Method(..),
+        Run(..),
+        runDynamics1,
+        runDynamics1_,
+        runDynamics,
+        runDynamics_,
+        runDynamicsIO,
+        runDynamicsSeries1,
+        runDynamicsSeries1_,
+        runDynamicsSeries,
+        runDynamicsSeries_,
+        printDynamics1,
+        printDynamics,
+        -- * Utilities
+        basicTime,
+        iterationBnds,
+        iterationHiBnd,
+        iterationLoBnd,
+        phaseBnds,
+        phaseHiBnd,
+        phaseLoBnd) where
+
+import Control.Monad
+import Control.Monad.Trans
+
+--
+-- The Dynamics Monad
+--
+-- A value of the Dynamics monad represents an abstract dynamic 
+-- process, i.e. a time varying polymorphic function. This is 
+-- a key point of the Aivika simulation library.
+--
+
+-- | A value in the 'Dynamics' monad represents a dynamic process, i.e.
+-- a polymorphic time varying function.
+newtype Dynamics a = Dynamics (Point -> IO a)
+
+-- | It defines the simulation point appended with the additional information.
+data Point = Point { pointSpecs :: Specs,    -- ^ the simulation specs
+                     pointRun :: Run,        -- ^ the simulation run
+                     pointTime :: Double,    -- ^ the current time
+                     pointIteration :: Int,  -- ^ the current iteration
+                     pointPhase :: Int       -- ^ the current phase
+                   } deriving (Eq, Ord, Show)
+
+-- | It defines the simulation specs.
+data Specs = Specs { spcStartTime :: Double,    -- ^ the start time
+                     spcStopTime :: Double,     -- ^ the stop time
+                     spcDT :: Double,           -- ^ the integration time step
+                     spcMethod :: Method        -- ^ the integration method
+                   } deriving (Eq, Ord, Show)
+
+-- | It defines the integration method.
+data Method = Euler          -- ^ Euler's method
+            | RungeKutta2    -- ^ the 2nd order Runge-Kutta method
+            | RungeKutta4    -- ^ the 4th order Runge-Kutta method
+            deriving (Eq, Ord, Show)
+
+-- | It defined the simulation run as part of some experiment.
+data Run = Run { runIndex :: Int,    -- ^ the current simulation run
+                 runCount :: Int     -- ^ the total number of runs in this experiment
+               } deriving (Eq, Ord, Show)
+           
+-- | Returns the iterations starting from zero.
+iterations :: Specs -> [Int]
+iterations sc = [i1 .. i2] where
+  i1 = 0
+  i2 = round ((spcStopTime sc - 
+               spcStartTime sc) / spcDT sc)
+
+-- | Returns the first and last iterations.
+iterationBnds :: Specs -> (Int, Int)
+iterationBnds sc = (0, round ((spcStopTime sc - 
+                               spcStartTime sc) / spcDT sc))
+
+-- | Returns the first iteration, i.e. zero.
+iterationLoBnd :: Specs -> Int
+iterationLoBnd sc = 0
+
+-- | Returns the last iteration.
+iterationHiBnd :: Specs -> Int
+iterationHiBnd sc = round ((spcStopTime sc - 
+                            spcStartTime sc) / spcDT sc)
+
+-- | Returns the phases for the specified simulation specs starting from zero.
+phases :: Specs -> [Int]
+phases sc = 
+  case spcMethod sc of
+    Euler -> [0]
+    RungeKutta2 -> [0, 1]
+    RungeKutta4 -> [0, 1, 2, 3]
+
+-- | Returns the first and last phases.
+phaseBnds :: Specs -> (Int, Int)
+phaseBnds sc = 
+  case spcMethod sc of
+    Euler -> (0, 0)
+    RungeKutta2 -> (0, 1)
+    RungeKutta4 -> (0, 3)
+
+-- | Returns the first phase, i.e. zero.
+phaseLoBnd :: Specs -> Int
+phaseLoBnd sc = 0
+                  
+-- | Returns the last phase, 1 for Euler's method, 2 for RK2 and 4 for RK4.
+phaseHiBnd :: Specs -> Int
+phaseHiBnd sc = 
+  case spcMethod sc of
+    Euler -> 0
+    RungeKutta2 -> 1
+    RungeKutta4 -> 3
+
+-- | Returns a simulation time for the integration point specified by 
+-- the specs, iteration and phase.
+basicTime :: Specs -> Int -> Int -> Double
+basicTime sc n ph =
+  if ph < 0 then 
+    error "Incorrect phase: basicTime"
+  else
+    spcStartTime sc + n' * spcDT sc + delta (spcMethod sc) ph 
+      where n' = fromInteger (toInteger n)
+            delta Euler       0 = 0
+            delta RungeKutta2 0 = 0
+            delta RungeKutta2 1 = spcDT sc
+            delta RungeKutta4 0 = 0
+            delta RungeKutta4 1 = spcDT sc / 2
+            delta RungeKutta4 2 = spcDT sc / 2
+            delta RungeKutta4 3 = spcDT sc
+
+instance Monad Dynamics where
+  return  = returnD
+  m >>= k = bindD m k
+
+returnD :: a -> Dynamics a
+returnD a = Dynamics (\p -> return a)
+
+bindD :: Dynamics a -> (a -> Dynamics b) -> Dynamics b
+bindD (Dynamics m) k = 
+  Dynamics $ \p -> 
+  do a <- m p
+     let Dynamics m' = k a
+     m' p
+
+subrunDynamics1 :: Dynamics a -> Specs -> Run -> IO a
+subrunDynamics1 (Dynamics m) sc r =
+  do let n = iterationHiBnd sc
+         t = basicTime sc n 0
+     m Point { pointSpecs = sc,
+               pointRun = r,
+               pointTime = t,
+               pointIteration = n,
+               pointPhase = 0 }
+
+subrunDynamics1_ :: Dynamics a -> Specs -> Run -> IO ()
+subrunDynamics1_ (Dynamics m) sc r =
+  do let n = iterationHiBnd sc
+         t = basicTime sc n 0
+     m Point { pointSpecs = sc,
+               pointRun = r,
+               pointTime = t,
+               pointIteration = n,
+               pointPhase = 0 }
+     return ()
+
+subrunDynamics :: Dynamics a -> Specs -> Run -> [IO a]
+subrunDynamics (Dynamics m) sc r =
+  do let (nl, nu) = iterationBnds sc
+         point n = Point { pointSpecs = sc,
+                           pointRun = r,
+                           pointTime = basicTime sc n 0,
+                           pointIteration = n,
+                           pointPhase = 0 }
+     map (m . point) [nl .. nu]
+
+subrunDynamics_ :: Dynamics a -> Specs -> Run -> IO ()
+subrunDynamics_ (Dynamics m) sc r =
+  do let (nl, nu) = iterationBnds sc
+         point n = Point { pointSpecs = sc,
+                           pointRun = r,
+                           pointTime = basicTime sc n 0,
+                           pointIteration = n,
+                           pointPhase = 0 }
+     mapM_ (m . point) [nl .. nu]
+
+-- | Run the simulation and return the result in the last 
+-- time point using the specified simulation specs.
+runDynamics1 :: Dynamics (Dynamics a) -> Specs -> IO a
+runDynamics1 (Dynamics m) sc = 
+  do let r = Run { runIndex = 1, runCount = 1 }
+     d <- m Point { pointSpecs = sc,
+                    pointRun = r,
+                    pointTime = spcStartTime sc,
+                    pointIteration = 0,
+                    pointPhase = 0 }
+     subrunDynamics1 d sc r
+
+-- | Run the simulation and return the result in the last 
+-- time point using the specified simulation specs.
+runDynamics1_ :: Dynamics (Dynamics a) -> Specs -> IO ()
+runDynamics1_ (Dynamics m) sc = 
+  do let r = Run { runIndex = 1, runCount = 1 }
+     d <- m Point { pointSpecs = sc,
+                    pointRun = r,
+                    pointTime = spcStartTime sc,
+                    pointIteration = 0,
+                    pointPhase = 0 }
+     subrunDynamics1_ d sc r
+
+-- | Run the simulation and return the results in all 
+-- integration time points using the specified simulation specs.
+runDynamics :: Dynamics (Dynamics a) -> Specs -> IO [a]
+runDynamics (Dynamics m) sc = 
+  do let r = Run { runIndex = 1, runCount = 1 }
+     d <- m Point { pointSpecs = sc,
+                    pointRun = r,
+                    pointTime = spcStartTime sc,
+                    pointIteration = 0,
+                    pointPhase = 0 }
+     sequence $ subrunDynamics d sc r
+
+-- | Run the simulation and return the results in all 
+-- integration time points using the specified simulation specs.
+runDynamics_ :: Dynamics (Dynamics a) -> Specs -> IO ()
+runDynamics_ (Dynamics m) sc = 
+  do let r = Run { runIndex = 1, runCount = 1 }
+     d <- m Point { pointSpecs = sc,
+                    pointRun = r,
+                    pointTime = spcStartTime sc,
+                    pointIteration = 0,
+                    pointPhase = 0 }
+     sequence_ $ subrunDynamics d sc r
+
+-- | Run the simulation and return the results in all 
+-- integration time points using the specified simulation specs.
+runDynamicsIO :: Dynamics (Dynamics a) -> Specs -> IO [IO a]
+runDynamicsIO (Dynamics m) sc =
+  do let r = Run { runIndex = 1, runCount = 1 }
+     d <- m Point { pointSpecs = sc,
+                    pointRun = r,
+                    pointTime = spcStartTime sc,
+                    pointIteration = 0,
+                    pointPhase = 0 }
+     return $ subrunDynamics d sc r
+
+-- | Run an experiment consisting of the given number of simulations, where each 
+-- model is created and then requested in the last integration time point using 
+-- the specified specs.
+runDynamicsSeries1_ :: Dynamics (Dynamics a) -> Specs -> Int -> [IO ()]
+runDynamicsSeries1_ (Dynamics m) sc runs = map f [1 .. runs]
+  where f i =
+          do let r = Run { runIndex = i, runCount = runs }
+             d <- m Point { pointSpecs = sc,
+                            pointRun = r,
+                            pointTime = spcStartTime sc,
+                            pointIteration = 0,
+                            pointPhase = 0 }
+             subrunDynamics1_ d sc r
+
+-- | Run an experiment consisting of the given number of simulations, where each 
+-- model is created and then requested sequentially in all integration time points 
+-- using the specified specs.
+runDynamicsSeries_ :: Dynamics (Dynamics a) -> Specs -> Int -> [IO ()]
+runDynamicsSeries_ (Dynamics m) sc runs = map f [1 .. runs]
+  where f i =
+          do let r = Run { runIndex = i, runCount = runs }
+             d <- m Point { pointSpecs = sc,
+                            pointRun = r,
+                            pointTime = spcStartTime sc,
+                            pointIteration = 0,
+                            pointPhase = 0 }
+             subrunDynamics_ d sc r
+
+-- | Run an experiment consisting of the given number of simulations, where each 
+-- model is created and then requested in the last integration time point using 
+-- the specified specs.
+runDynamicsSeries1 :: Dynamics (Dynamics a) -> Specs -> Int -> [IO a]
+runDynamicsSeries1 (Dynamics m) sc runs = map f [1 .. runs]
+  where f i =
+          do let r = Run { runIndex = i, runCount = runs }
+             d <- m Point { pointSpecs = sc,
+                            pointRun = r,
+                            pointTime = spcStartTime sc,
+                            pointIteration = 0,
+                            pointPhase = 0 }
+             subrunDynamics1 d sc r
+
+-- | Run an experiment consisting of the given number of simulations, where each 
+-- model is created and then requested sequentially in all integration time points 
+-- using the specified specs.
+runDynamicsSeries :: Dynamics (Dynamics a) -> Specs -> Int -> [IO [a]]
+runDynamicsSeries (Dynamics m) sc runs = map f [1 .. runs]
+  where f i =
+          do let r = Run { runIndex = i, runCount = runs }
+             d <- m Point { pointSpecs = sc,
+                            pointRun = r,
+                            pointTime = spcStartTime sc,
+                            pointIteration = 0,
+                            pointPhase = 0 }
+             sequence $ subrunDynamics d sc r
+
+-- | Run the simulation and print the result in the last 
+-- time point using the specified simulation specs.
+printDynamics1 :: (Show a) => Dynamics (Dynamics a) -> Specs -> IO ()
+printDynamics1 m sc = runDynamics1 m sc >>= print
+
+-- | Run the simulation and print lazily the results in all
+-- integration time points using the specified simulation specs.
+printDynamics :: (Show a) => Dynamics (Dynamics a) -> Specs -> IO ()
+printDynamics m sc = runDynamicsIO m sc >>= loop
+  where loop [] = return ()
+        loop (x : xs) = do { a <- x; print a; loop xs }
+
+instance Functor Dynamics where
+  fmap = liftMD
+
+instance Eq (Dynamics a) where
+  x == y = error "Can't compare dynamics." 
+
+instance Show (Dynamics a) where
+  showsPrec _ x = showString "<< Dynamics >>"
+
+liftMD :: (a -> b) -> Dynamics a -> Dynamics b
+{-# INLINE liftMD #-}
+liftMD f (Dynamics x) =
+  Dynamics $ \p -> do { a <- x p; return $ f a }
+
+liftM2D :: (a -> b -> c) -> Dynamics a -> Dynamics b -> Dynamics c
+{-# INLINE liftM2D #-}
+liftM2D f (Dynamics x) (Dynamics y) =
+  Dynamics $ \p -> do { a <- x p; b <- y p; return $ f a b }
+
+instance (Num a) => Num (Dynamics a) where
+  x + y = liftM2D (+) x y
+  x - y = liftM2D (-) x y
+  x * y = liftM2D (*) x y
+  negate = liftMD negate
+  abs = liftMD abs
+  signum = liftMD signum
+  fromInteger i = return $ fromInteger i
+
+instance (Fractional a) => Fractional (Dynamics a) where
+  x / y = liftM2D (/) x y
+  recip = liftMD recip
+  fromRational t = return $ fromRational t
+
+instance (Floating a) => Floating (Dynamics a) where
+  pi = return pi
+  exp = liftMD exp
+  log = liftMD log
+  sqrt = liftMD sqrt
+  x ** y = liftM2D (**) x y
+  sin = liftMD sin
+  cos = liftMD cos
+  tan = liftMD tan
+  asin = liftMD asin
+  acos = liftMD acos
+  atan = liftMD atan
+  sinh = liftMD sinh
+  cosh = liftMD cosh
+  tanh = liftMD tanh
+  asinh = liftMD asinh
+  acosh = liftMD acosh
+  atanh = liftMD atanh
+
+instance MonadIO Dynamics where
+  liftIO m = Dynamics $ const m
diff --git a/Simulation/Aivika/Dynamics/Internal/Fold.hs b/Simulation/Aivika/Dynamics/Internal/Fold.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Dynamics/Internal/Fold.hs
@@ -0,0 +1,91 @@
+
+-- |
+-- Module     : Simulation.Aivika.Dynamics.Internal.Fold
+-- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.0.3
+--
+-- This module defines the fold functions that allows traversing the values of
+-- any dynamic process in the integration time points.
+--
+module Simulation.Aivika.Dynamics.Internal.Fold
+       (foldD1,
+        foldD,
+        divideD) where
+
+import Data.IORef
+import Control.Monad
+import Control.Monad.Trans
+
+import Simulation.Aivika.Dynamics.Internal.Dynamics
+import Simulation.Aivika.Dynamics.Internal.Interpolate
+import Simulation.Aivika.Dynamics.Internal.Memo
+
+--
+-- Fold
+--
+
+-- | Like the standard 'foldl1' function but applied to values in 
+-- the integration time points. The accumulator values are transformed
+-- according to the first argument, which should be either function 
+-- 'memo0' or 'umemo0'.
+foldD1 :: (Dynamics a -> Dynamics (Dynamics a))
+         -> (a -> a -> a) 
+         -> Dynamics a 
+         -> Dynamics (Dynamics a)
+foldD1 tr f (Dynamics m) =
+  do r <- liftIO $ newIORef m
+     let z = Dynamics $ \p ->
+           case pointIteration p of
+             0 -> 
+               m p
+             n -> do 
+               let sc = pointSpecs p
+                   ty = basicTime sc (n - 1) 0
+                   py = p { pointTime = ty, pointIteration = n - 1, pointPhase = 0 }
+               y <- readIORef r
+               s <- y py
+               x <- m p
+               return $! f s x
+     y@(Dynamics m) <- tr z
+     liftIO $ writeIORef r m
+     return y
+
+-- | Like the standard 'foldl' function but applied to values in 
+-- the integration time points. The accumulator values are transformed
+-- according to the first argument, which should be either function
+-- 'memo0' or 'umemo0'.
+foldD :: (Dynamics a -> Dynamics (Dynamics a))
+        -> (a -> b -> a) 
+        -> a
+        -> Dynamics b 
+        -> Dynamics (Dynamics a)
+foldD tr f acc (Dynamics m) =
+  do r <- liftIO $ newIORef $ const $ return acc
+     let z = Dynamics $ \p ->
+           case pointIteration p of
+             0 -> do
+               x <- m p
+               return $! f acc x
+             n -> do 
+               let sc = pointSpecs p
+                   ty = basicTime sc (n - 1) 0
+                   py = p { pointTime = ty, pointIteration = n - 1, pointPhase = 0 }
+               y <- readIORef r
+               s <- y py
+               x <- m p
+               return $! f s x
+     y@(Dynamics m) <- tr z
+     liftIO $ writeIORef r m
+     return y
+
+-- | Divide the values in integration time points by the number of
+-- the current iteration. It can be useful for statistic functions in
+-- combination with the fold.
+divideD :: Dynamics Double -> Dynamics Double
+divideD (Dynamics m) = 
+  discrete $ Dynamics $ \p ->
+  do a <- m p
+     return $ a / fromInteger (toInteger (pointIteration p + 1))
diff --git a/Simulation/Aivika/Dynamics/Internal/Interpolate.hs b/Simulation/Aivika/Dynamics/Internal/Interpolate.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Dynamics/Internal/Interpolate.hs
@@ -0,0 +1,62 @@
+
+-- |
+-- Module     : Simulation.Aivika.Dynamics.Internal.Interpolate
+-- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.0.3
+--
+-- This module defines interpolation functions.
+-- These functions complement the memoization, possibly except for 
+-- the 'initD' function which is useful to get an initial 
+-- value of any dynamic process.
+--
+
+module Simulation.Aivika.Dynamics.Internal.Interpolate
+       (initD,
+        discrete,
+        interpolate) where
+
+import Simulation.Aivika.Dynamics.Internal.Dynamics
+
+-- | Return the initial value.
+initD :: Dynamics a -> Dynamics a
+{-# INLINE initD #-}
+initD (Dynamics m) =
+  Dynamics $ \p ->
+  if pointIteration p == 0 && pointPhase p == 0 then
+    m p
+  else
+    let sc = pointSpecs p
+    in m $ p { pointTime = basicTime sc 0 0,
+               pointIteration = 0,
+               pointPhase = 0 } 
+
+-- | Discretize the computation in the integration time points.
+discrete :: Dynamics a -> Dynamics a
+{-# INLINE discrete #-}
+discrete (Dynamics m) =
+  Dynamics $ \p ->
+  if pointPhase p == 0 then
+    m p
+  else
+    let sc = pointSpecs p
+        n  = pointIteration p
+    in m $ p { pointTime = basicTime sc n 0,
+               pointPhase = 0 }
+
+-- | Interpolate the computation based on the integration time points only.
+-- Unlike the 'discrete' function it knows about the intermediate time points 
+-- that are used in the Runge-Kutta method.
+interpolate :: Dynamics a -> Dynamics a
+{-# INLINE interpolate #-}
+interpolate (Dynamics m) = 
+  Dynamics $ \p -> 
+  if pointPhase p >= 0 then 
+    m p
+  else 
+    let sc = pointSpecs p
+        n  = pointIteration p
+    in m $ p { pointTime = basicTime sc n 0,
+               pointPhase = 0 }
diff --git a/Simulation/Aivika/Dynamics/Internal/Memo.hs b/Simulation/Aivika/Dynamics/Internal/Memo.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Dynamics/Internal/Memo.hs
@@ -0,0 +1,195 @@
+
+{-# LANGUAGE FlexibleContexts #-}
+
+-- |
+-- Module     : Simulation.Aivika.Dynamics.Internal.Memo
+-- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.0.3
+--
+-- This module defines memo functions. The memoization creates such dynamic processes, 
+-- which values are cached in the integration time points. Then these values are 
+-- interpolated in all other time points.
+--
+
+module Simulation.Aivika.Dynamics.Internal.Memo
+       (memo,
+        umemo,
+        memo0,
+        umemo0,
+        iterateD) where
+
+import Data.Array
+import Data.Array.IO
+import Data.IORef
+import Control.Monad
+
+import Simulation.Aivika.Dynamics.Internal.Dynamics
+import Simulation.Aivika.Dynamics.Internal.Interpolate
+
+newMemoArray_ :: Ix i => (i, i) -> IO (IOArray i e)
+newMemoArray_ = newArray_
+
+newMemoUArray_ :: (MArray IOUArray e IO, Ix i) => (i, i) -> IO (IOUArray i e)
+newMemoUArray_ = newArray_
+
+-- | Memoize and order the computation in the integration time points using 
+-- the interpolation that knows of the Runge-Kutta method.
+memo :: Dynamics e -> Dynamics (Dynamics e)
+{-# INLINE memo #-}
+memo (Dynamics m) = 
+  Dynamics $ \p ->
+  do let sc = pointSpecs p
+         (phl, phu) = phaseBnds sc
+         (nl, nu)   = iterationBnds sc
+     arr   <- newMemoArray_ ((phl, nl), (phu, nu))
+     nref  <- newIORef 0
+     phref <- newIORef 0
+     let r p = 
+           do let sc  = pointSpecs p
+                  n   = pointIteration p
+                  ph  = pointPhase p
+                  phu = phaseHiBnd sc 
+                  loop n' ph' = 
+                    if (n' > n) || ((n' == n) && (ph' > ph)) 
+                    then 
+                      readArray arr (ph, n)
+                    else 
+                      let p' = p { pointIteration = n', pointPhase = ph',
+                                   pointTime = basicTime sc n' ph' }
+                      in do a <- m p'
+                            a `seq` writeArray arr (ph', n') a
+                            if ph' >= phu 
+                              then do writeIORef phref 0
+                                      writeIORef nref (n' + 1)
+                                      loop (n' + 1) 0
+                              else do writeIORef phref (ph' + 1)
+                                      loop n' (ph' + 1)
+              n'  <- readIORef nref
+              ph' <- readIORef phref
+              loop n' ph'
+     return $ interpolate $ Dynamics r
+
+-- | This is a more efficient version the 'memo' function which uses 
+-- an unboxed array to store the values.
+umemo :: (MArray IOUArray e IO) => Dynamics e -> Dynamics (Dynamics e)
+{-# INLINE umemo #-}
+umemo (Dynamics m) = 
+  Dynamics $ \p ->
+  do let sc = pointSpecs p
+         (phl, phu) = phaseBnds sc
+         (nl, nu)   = iterationBnds sc
+     arr   <- newMemoUArray_ ((phl, nl), (phu, nu))
+     nref  <- newIORef 0
+     phref <- newIORef 0
+     let r p =
+           do let sc  = pointSpecs p
+                  n   = pointIteration p
+                  ph  = pointPhase p
+                  phu = phaseHiBnd sc 
+                  loop n' ph' = 
+                    if (n' > n) || ((n' == n) && (ph' > ph)) 
+                    then 
+                      readArray arr (ph, n)
+                    else 
+                      let p' = p { pointIteration = n', 
+                                   pointPhase = ph',
+                                   pointTime = basicTime sc n' ph' }
+                      in do a <- m p'
+                            a `seq` writeArray arr (ph', n') a
+                            if ph' >= phu 
+                              then do writeIORef phref 0
+                                      writeIORef nref (n' + 1)
+                                      loop (n' + 1) 0
+                              else do writeIORef phref (ph' + 1)
+                                      loop n' (ph' + 1)
+              n'  <- readIORef nref
+              ph' <- readIORef phref
+              loop n' ph'
+     return $ interpolate $ Dynamics r
+
+-- | Memoize and order the computation in the integration time points using 
+-- the 'discrete' interpolation. It consumes less memory than the 'memo'
+-- function but it is not aware of the Runge-Kutta method. There is a subtle
+-- difference when we request for values in the intermediate time points
+-- that are used by this method to integrate. In general case you should 
+-- prefer the 'memo0' function above 'memo'.
+memo0 :: Dynamics e -> Dynamics (Dynamics e)
+{-# INLINE memo0 #-}
+memo0 (Dynamics m) = 
+  Dynamics $ \p ->
+  do let sc   = pointSpecs p
+         bnds = iterationBnds sc
+     arr  <- newMemoArray_ bnds
+     nref <- newIORef 0
+     let r p =
+           do let sc = pointSpecs p
+                  n  = pointIteration p
+                  loop n' = 
+                    if n' > n
+                    then 
+                      readArray arr n
+                    else 
+                      let p' = p { pointIteration = n', pointPhase = 0,
+                                   pointTime = basicTime sc n' 0 }
+                      in do a <- m p'
+                            a `seq` writeArray arr n' a
+                            writeIORef nref (n' + 1)
+                            loop (n' + 1)
+              n' <- readIORef nref
+              loop n'
+     return $ discrete $ Dynamics r
+
+-- | This is a more efficient version the 'memo0' function which uses 
+-- an unboxed array to store the values.
+umemo0 :: (MArray IOUArray e IO) => Dynamics e -> Dynamics (Dynamics e)
+{-# INLINE umemo0 #-}
+umemo0 (Dynamics m) = 
+  Dynamics $ \p ->
+  do let sc   = pointSpecs p
+         bnds = iterationBnds sc
+     arr  <- newMemoUArray_ bnds
+     nref <- newIORef 0
+     let r p =
+           do let sc = pointSpecs p
+                  n  = pointIteration p
+                  loop n' = 
+                    if n' > n
+                    then 
+                      readArray arr n
+                    else 
+                      let p' = p { pointIteration = n', pointPhase = 0,
+                                   pointTime = basicTime sc n' 0 }
+                      in do a <- m p'
+                            a `seq` writeArray arr n' a
+                            writeIORef nref (n' + 1)
+                            loop (n' + 1)
+              n' <- readIORef nref
+              loop n'
+     return $ discrete $ Dynamics r
+
+-- | Iterate sequentially the dynamic process with side effects in 
+-- the integration time points. It is equivalent to a call of the
+-- 'memo0' function but significantly more efficient, for the array 
+-- is not created.
+iterateD :: Dynamics () -> Dynamics (Dynamics ())
+{-# INLINE iterateD #-}
+iterateD (Dynamics m) = 
+  Dynamics $ \p ->
+  do let sc = pointSpecs p
+     nref <- newIORef 0
+     let r p =
+           do let sc = pointSpecs p
+                  n  = pointIteration p
+                  loop n' = 
+                    unless (n' > n) $
+                    let p' = p { pointIteration = n', pointPhase = 0,
+                                 pointTime = basicTime sc n' 0 }
+                    in do a <- m p'
+                          a `seq` writeIORef nref (n' + 1)
+                          loop (n' + 1)
+              n' <- readIORef nref
+              loop n'
+     return $ discrete $ Dynamics r
diff --git a/Simulation/Aivika/Dynamics/Internal/Process.hs b/Simulation/Aivika/Dynamics/Internal/Process.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Dynamics/Internal/Process.hs
@@ -0,0 +1,162 @@
+
+-- |
+-- Module     : Simulation.Aivika.Dynamics.Internal.Process
+-- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.0.3
+--
+-- A value in the 'Process' monad represents a discontinuous process that 
+-- can suspend and resume at any time. It behaves like a dynamic process too. 
+-- Any value in the 'Dynamics' monad can be lifted to the Process monad. 
+-- Moreover, a value in the Process monad can be run in the Dynamics monad.
+--
+-- A value of the 'ProcessID' type is just an identifier of such a process.
+--
+module Simulation.Aivika.Dynamics.Internal.Process
+       (ProcessID,
+        Process(..),
+        processQueue,
+        newProcessID,
+        holdProcess,
+        passivateProcess,
+        processPassive,
+        reactivateProcess,
+        processID,
+        runProcess) where
+
+import Data.IORef
+import Control.Monad
+import Control.Monad.Trans
+
+import Simulation.Aivika.Dynamics.Internal.Dynamics
+import Simulation.Aivika.Dynamics.Internal.Cont
+import Simulation.Aivika.Dynamics.Lift
+import Simulation.Aivika.Dynamics.EventQueue
+
+-- | Represents a process identificator.
+data ProcessID = 
+  ProcessID { processQueue   :: EventQueue,  -- ^ Return the event queue.
+              processStarted :: IORef Bool,
+              processCont    :: IORef (Maybe (Dynamics (() -> IO ()))) }
+
+-- | Specifies a discontinuous process that can suspend at any time
+-- and then resume later.
+newtype Process a = Process (ProcessID -> Cont a)
+
+-- | Hold the process for the specified time period.
+holdProcess :: Double -> Process ()
+holdProcess dt =
+  Process $ \pid ->
+  Cont $ \c ->
+  Dynamics $ \ps ->
+  do let Dynamics m = enqueueCont (processQueue pid) (pointTime ps + dt) c
+     m ps
+
+-- | Passivate the process.
+passivateProcess :: Process ()
+passivateProcess =
+  Process $ \pid ->
+  Cont $ \c ->
+  Dynamics $ \p ->
+  do let x = processCont pid
+     a <- readIORef x
+     case a of
+       Nothing -> writeIORef x $ Just c
+       Just _  -> error "Cannot passivate the process twice: passivate"
+
+-- | Test whether the process with the specified ID is passivated.
+processPassive :: ProcessID -> Process Bool
+processPassive pid =
+  Process $ \_ ->
+  Cont $ \(Dynamics c) ->
+  Dynamics $ \p ->
+  do cont' <- c p
+     let x = processCont pid
+     a <- readIORef x
+     case a of
+       Nothing -> cont' False
+       Just _  -> cont' True
+
+-- | Reactivate a process with the specified ID.
+reactivateProcess :: ProcessID -> Process ()
+reactivateProcess pid =
+  Process $ \pid' ->
+  Cont $ \c@(Dynamics cont) ->
+  Dynamics $ \p ->
+  do let x = processCont pid
+     a <- readIORef x
+     case a of
+       Nothing ->
+         do cont' <- cont p
+            cont' ()
+       Just (Dynamics cont2) ->
+         do writeIORef x Nothing
+            let Dynamics m = enqueueCont (processQueue pid') (pointTime p) c
+            m p
+            cont2' <- cont2 p
+            cont2' ()
+
+-- | Start the process with the specified ID at the desired time.
+runProcess :: Process () -> ProcessID -> Double -> Dynamics ()
+runProcess (Process p) pid t =
+  runCont m r
+    where m = do y <- liftIO $ readIORef (processStarted pid)
+                 if y 
+                   then error $
+                        "A process with such ID " ++
+                        "has been started already: runProc"
+                   else liftIO $ writeIORef (processStarted pid) True
+                 Cont $ \c -> enqueueCont (processQueue pid) t c
+                 p pid
+          r = let f () = return () in return f
+
+-- | Return the current process ID.
+processID :: Process ProcessID
+processID = Process $ \pid -> return pid
+
+-- | Create a new process ID.
+newProcessID :: EventQueue -> Dynamics ProcessID
+newProcessID q =
+  do x <- liftIO $ newIORef Nothing
+     y <- liftIO $ newIORef False
+     return ProcessID { processQueue   = q,
+                        processStarted = y,
+                        processCont    = x }
+
+instance Eq ProcessID where
+  x == y = processCont x == processCont y    -- for the references are unique
+
+instance Monad Process where
+  return  = returnP
+  m >>= k = bindP m k
+
+instance Functor Process where
+  fmap = liftM
+
+instance Lift Process where
+  liftD = liftP
+  
+instance MonadIO Process where
+  liftIO = liftIOP
+  
+returnP :: a -> Process a
+{-# INLINE returnP #-}
+returnP a = Process (\pid -> return a)
+
+bindP :: Process a -> (a -> Process b) -> Process b
+{-# INLINE bindP #-}
+bindP (Process m) k = 
+  Process $ \pid -> 
+  do a <- m pid
+     let Process m' = k a
+     m' pid
+
+liftP :: Dynamics a -> Process a
+{-# INLINE liftP #-}
+liftP m = Process $ \pid -> liftD m
+
+liftIOP :: IO a -> Process a
+{-# INLINE liftIOP #-}
+liftIOP m = Process $ \pid -> liftIO m
diff --git a/Simulation/Aivika/Dynamics/Internal/Time.hs b/Simulation/Aivika/Dynamics/Internal/Time.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Dynamics/Internal/Time.hs
@@ -0,0 +1,35 @@
+
+-- |
+-- Module     : Simulation.Aivika.Dynamics.Internal.Time
+-- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.0.3
+--
+-- This module defines the time parameters.
+--
+
+module Simulation.Aivika.Dynamics.Internal.Time
+       (starttime,
+        stoptime,
+        dt,
+        time) where
+
+import Simulation.Aivika.Dynamics.Internal.Dynamics
+
+-- | Return the start simulation time.
+starttime :: Dynamics Double
+starttime = Dynamics $ return . spcStartTime . pointSpecs
+
+-- | Return the stop simulation time.
+stoptime :: Dynamics Double
+stoptime = Dynamics $ return . spcStopTime . pointSpecs
+
+-- | Return the integration time step.
+dt :: Dynamics Double
+dt = Dynamics $ return . spcDT . pointSpecs
+
+-- | Return the current simulation time.
+time :: Dynamics Double
+time = Dynamics $ return . pointTime 
diff --git a/Simulation/Aivika/Dynamics/Lift.hs b/Simulation/Aivika/Dynamics/Lift.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Dynamics/Lift.hs
@@ -0,0 +1,21 @@
+
+-- |
+-- Module     : Simulation.Aivika.Dynamics.Lift
+-- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.0.3
+--
+-- This module defines the 'liftD' function that allows embedding
+-- the 'Dynamics' computation.
+--
+module Simulation.Aivika.Dynamics.Lift (Lift(..)) where
+
+import Simulation.Aivika.Dynamics
+
+-- | The 'Lift' class defines a type which the 'Dynamics' 
+-- computation can be lifted to.
+class Lift m where
+  -- | Lift the computation.
+  liftD :: Dynamics a -> m a
diff --git a/Simulation/Aivika/Dynamics/Parameter.hs b/Simulation/Aivika/Dynamics/Parameter.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Dynamics/Parameter.hs
@@ -0,0 +1,58 @@
+
+-- |
+-- Module     : Simulation.Aivika.Dynamics.Parameter
+-- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.0.3
+--
+-- This module defines the parameters of simulation experiments.
+--
+
+module Simulation.Aivika.Dynamics.Parameter
+       (newParameter,
+        newTableParameter,
+        newIndexedParameter) where
+
+import Data.Array
+import Data.IORef
+import qualified Data.Map as M
+import Control.Concurrent.MVar
+
+import Simulation.Aivika.Dynamics.Internal.Dynamics
+
+-- | Create a thread-safe parameter that returns always the same value during the simulation run, 
+-- where the value is recalculated for each new run.
+newParameter :: IO a -> IO (Dynamics a)
+newParameter a = newIndexedParameter $ \_ -> a
+
+-- | Create a thread-safe parameter that returns always the same value during the simulation run,
+-- where the value is taken consequently from the specified table based on the number of the 
+-- current run starting from zero. After all values from the table are used, it takes the first 
+-- value of the table, then the second one and so on.
+newTableParameter :: Array Int a -> IO (Dynamics a)
+newTableParameter t = newIndexedParameter (\i -> return $ t ! (((i - i1) `mod` n) + i1))
+  where (i1, i2) = bounds t
+        n = i2 - i1 + 1
+
+-- | Create a thread-safe parameter that returns always the same value during the simulation run, 
+-- where the value depends on the number of this run starting from zero.
+newIndexedParameter :: (Int -> IO a) -> IO (Dynamics a)
+newIndexedParameter f = 
+  do lock <- newMVar ()
+     dict <- newIORef M.empty
+     return $ Dynamics $ \p ->
+       do let i = runIndex $ pointRun p
+          m <- readIORef dict
+          if M.member i m
+            then do let Just v = M.lookup i m
+                    return v
+            else withMVar lock $ 
+                 \() -> do { m <- readIORef dict;
+                            if M.member i m
+                            then do let Just v = M.lookup i m
+                                    return v
+                            else do v <- f i
+                                    writeIORef dict $ M.insert i v m
+                                    return v }
diff --git a/Simulation/Aivika/Dynamics/Process.hs b/Simulation/Aivika/Dynamics/Process.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Dynamics/Process.hs
@@ -0,0 +1,30 @@
+
+-- |
+-- Module     : Simulation.Aivika.Dynamics.Process
+-- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.0.3
+--
+-- A value in the 'Process' monad represents a discontinuous process that 
+-- can suspend and resume at any time. It behaves like a dynamic process too. 
+-- Any value in the 'Dynamics' monad can be lifted to the Process monad. 
+-- Moreover, a value in the Process monad can be run in the Dynamics monad.
+--
+-- A value of the 'ProcessID' type is just an identifier of such a process.
+--
+module Simulation.Aivika.Dynamics.Process
+       (ProcessID,
+        Process,
+        processQueue,
+        newProcessID,
+        holdProcess,
+        passivateProcess,
+        processPassive,
+        reactivateProcess,
+        processID,
+        runProcess) where
+
+import Simulation.Aivika.Dynamics.Internal.Dynamics
+import Simulation.Aivika.Dynamics.Internal.Process
diff --git a/Simulation/Aivika/Dynamics/Random.hs b/Simulation/Aivika/Dynamics/Random.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Dynamics/Random.hs
@@ -0,0 +1,73 @@
+
+-- |
+-- Module     : Simulation.Aivika.Dynamics.Random
+-- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.0.3
+--
+-- Below are defined random functions that mostly return discrete processes. 
+-- Literally, it means that the values are initially defined in integration 
+-- time points and then they are passed to the 'discrete' function.
+--
+
+module Simulation.Aivika.Dynamics.Random 
+       (newRandom, newNormal, normalGen) where
+
+import Random
+import Data.IORef
+import Control.Monad.Trans
+
+import Simulation.Aivika.Dynamics
+import Simulation.Aivika.Dynamics.Base
+
+-- | Return the uniform random numbers between 0.0 and 1.0 in
+-- the integration time points.
+newRandom :: Dynamics (Dynamics Double)
+newRandom =
+  memo0 $ liftIO $ getStdRandom random
+     
+-- | Return the normal random numbers with mean 0.0 and variance 1.0 in
+-- the integration time points.
+newNormal :: Dynamics (Dynamics Double)
+newNormal =
+  do g <- liftIO normalGen
+     memo0 $ liftIO g
+
+-- | Normal random number generator.
+normalGen :: IO (IO Double)
+normalGen =
+  do nextRef <- newIORef 0.0
+     flagRef <- newIORef False
+     xi1Ref  <- newIORef 0.0
+     xi2Ref  <- newIORef 0.0
+     psiRef  <- newIORef 0.0
+     let loop =
+           do psi <- readIORef psiRef
+              if (psi >= 1.0) || (psi == 0.0)
+                then do g1 <- getStdRandom random
+                        g2 <- getStdRandom random
+                        let xi1 = 2.0 * g1 - 1.0
+                            xi2 = 2.0 * g2 - 1.0
+                            psi = xi1 * xi1 + xi2 * xi2
+                        writeIORef xi1Ref xi1
+                        writeIORef xi2Ref xi2
+                        writeIORef psiRef psi
+                        loop
+                else writeIORef psiRef $ sqrt (- 2.0 * log psi / psi)
+     return $
+       do flag <- readIORef flagRef
+          if flag
+            then do writeIORef flagRef False
+                    readIORef nextRef
+            else do writeIORef xi1Ref 0.0
+                    writeIORef xi2Ref 0.0
+                    writeIORef psiRef 0.0
+                    loop
+                    xi1 <- readIORef xi1Ref
+                    xi2 <- readIORef xi2Ref
+                    psi <- readIORef psiRef
+                    writeIORef flagRef True
+                    writeIORef nextRef $ xi2 * psi
+                    return $ xi1 * psi
diff --git a/Simulation/Aivika/Dynamics/Ref.hs b/Simulation/Aivika/Dynamics/Ref.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Dynamics/Ref.hs
@@ -0,0 +1,63 @@
+
+-- |
+-- Module     : Simulation.Aivika.Dynamics.Ref
+-- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.0.3
+--
+-- This module defines an updatable reference that depends on the event queue.
+--
+module Simulation.Aivika.Dynamics.Ref
+       (Ref,
+        newRef,
+        refQueue,
+        readRef,
+        writeRef,
+        modifyRef) where
+
+import Data.IORef
+import Control.Monad.Trans
+
+import Simulation.Aivika.Dynamics.Internal.Dynamics
+import Simulation.Aivika.Dynamics.EventQueue
+
+-- | The 'Ref' type represents a mutable variable similar to the 'IORef' variable 
+-- but only bound to some event queue, which makes the variable coordinated 
+-- with that queue.
+data Ref a = 
+  Ref { refQueue :: EventQueue,  -- ^ Return the bound event queue.
+        refRun   :: Dynamics (),
+        refValue :: IORef a }
+
+-- | Create a new reference bound to the specified event queue.
+newRef :: EventQueue -> a -> Dynamics (Ref a)
+newRef q a =
+  do x <- liftIO $ newIORef a
+     return Ref { refQueue = q,
+                  refRun   = queueRun q,
+                  refValue = x }
+     
+-- | Read the value of a reference, forcing the bound event queue to raise 
+-- the events in case of need.
+readRef :: Ref a -> Dynamics a
+readRef r = Dynamics $ \p -> 
+  do let Dynamics m = refRun r
+     m p
+     readIORef (refValue r)
+
+-- | Write a new value into the reference.
+writeRef :: Ref a -> a -> Dynamics ()
+writeRef r a = Dynamics $ \p -> 
+  a `seq` writeIORef (refValue r) a
+
+-- | Mutate the contents of the reference, forcing the bound event queue to
+-- raise all pending events in case of need.
+modifyRef :: Ref a -> (a -> a) -> Dynamics ()
+modifyRef r f = Dynamics $ \p -> 
+  do let Dynamics m = refRun r
+     m p
+     a <- readIORef (refValue r)
+     let b = f a
+     b `seq` writeIORef (refValue r) b
diff --git a/Simulation/Aivika/Dynamics/Resource.hs b/Simulation/Aivika/Dynamics/Resource.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Dynamics/Resource.hs
@@ -0,0 +1,101 @@
+
+-- |
+-- Module     : Simulation.Aivika.Dynamics.Resource
+-- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.0.3
+--
+-- This module defines a limited resource which can be acquired and 
+-- then released by the discontinuous process 'DynamicProc'.
+--
+module Simulation.Aivika.Dynamics.Resource
+       (Resource,
+        newResource,
+        resourceQueue,
+        resourceInitCount,
+        resourceCount,
+        requestResource,
+        releaseResource) where
+
+import Data.IORef
+import Control.Monad
+
+import Simulation.Aivika.Dynamics.Internal.Dynamics
+import Simulation.Aivika.Dynamics.Internal.Cont
+import Simulation.Aivika.Dynamics.Internal.Process
+import Simulation.Aivika.Dynamics.EventQueue
+import qualified Simulation.Aivika.Queue as Q
+
+-- | Represents a limited resource.
+data Resource = 
+  Resource { resourceQueue     :: EventQueue,  
+             -- ^ Return the bound event queue.
+             resourceInitCount :: Int,
+             -- ^ Return the initial count of the resource.
+             resourceCountRef  :: IORef Int, 
+             resourceWaitQueue :: Q.Queue (Dynamics (() -> IO ()))}
+
+instance Eq Resource where
+  x == y = resourceCountRef x == resourceCountRef y  -- unique references
+
+-- | Create a new resource with the specified initial count.
+newResource :: EventQueue -> Int -> Dynamics Resource
+newResource q initCount =
+  Dynamics $ \p ->
+  do countRef  <- newIORef initCount
+     waitQueue <- Q.newQueue
+     return Resource { resourceQueue     = q,
+                       resourceInitCount = initCount,
+                       resourceCountRef  = countRef,
+                       resourceWaitQueue = waitQueue }
+
+-- | Return the current count of the resource.
+resourceCount :: Resource -> Process Int
+resourceCount r =
+  Process $ \_ ->
+  Cont $ \(Dynamics c) ->
+  Dynamics $ \p ->
+  do cont' <- c p 
+     a <- readIORef (resourceCountRef r)
+     cont' a
+
+-- | Request for the resource decreasing its count in case of success,
+-- otherwise suspending the discontinuous process until some other 
+-- process releases the resource.
+requestResource :: Resource -> Process ()
+requestResource r =
+  Process $ \_ ->
+  Cont $ \c@(Dynamics cont) ->
+  Dynamics $ \p ->
+  do a <- readIORef (resourceCountRef r)
+     if a == 0 
+       then Q.enqueue (resourceWaitQueue r) c
+       else do let a' = a - 1
+               a' `seq` writeIORef (resourceCountRef r) a'
+               cont' <- cont p
+               cont' ()
+
+-- | Release the resource increasing its count and resuming one of the
+-- previously suspended processes as possible.
+releaseResource :: Resource -> Process ()
+releaseResource r =
+  Process $ \_ ->
+  Cont $ \(Dynamics c) ->
+  Dynamics $ \p ->
+  do a <- readIORef (resourceCountRef r)
+     let a' = a + 1
+     when (a' > resourceInitCount r) $
+       error $
+       "The resource count cannot be greater than " ++
+       "its initial value: releaseResource."
+     f <- Q.queueNull (resourceWaitQueue r)
+     if f 
+       then a' `seq` writeIORef (resourceCountRef r) a'
+       else do c2 <- Q.queueFront (resourceWaitQueue r)
+               Q.dequeue (resourceWaitQueue r)
+               let Dynamics m = enqueueCont (resourceQueue r) (pointTime p) c2
+               m p
+     cont' <- c p
+     cont' ()
diff --git a/Simulation/Aivika/Dynamics/SystemDynamics.hs b/Simulation/Aivika/Dynamics/SystemDynamics.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Dynamics/SystemDynamics.hs
@@ -0,0 +1,441 @@
+
+{-# LANGUAGE FlexibleContexts #-}
+
+-- |
+-- Module     : Simulation.Aivika.Dynamics.SystemDynamics
+-- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.0.3
+--
+-- This module defines integrals and other functions of System Dynamics.
+--
+
+module Simulation.Aivika.Dynamics.SystemDynamics
+       (-- * Maximum and Minimum
+        maxD,
+        minD,
+        -- * Integrals
+        Integ,
+        newInteg,
+        integInit,
+        integValue,
+        integDiff,
+        -- * Integral Functions
+        integ,
+        -- * Difference Equations
+        Sum,
+        newSum,
+        sumInit,
+        sumValue,
+        sumDiff,
+        -- * Table Functions
+        lookupD,
+        lookupStepwiseD) where
+
+import Data.Array
+import Data.Array.IO
+import Data.IORef
+import Control.Monad
+import Control.Monad.Trans
+
+import Simulation.Aivika.Dynamics.Internal.Dynamics
+import Simulation.Aivika.Dynamics.Base
+
+--
+-- Maximum and Minimum
+--
+
+-- | Return the maximum.
+maxD :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics a
+maxD = liftM2 max
+
+-- | Return the minimum.
+minD :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics a
+minD = liftM2 min
+
+--
+-- Integrals
+--
+
+-- | The 'Integ' type represents an integral.
+data Integ = Integ { integInit     :: Dynamics Double,   -- ^ The initial value.
+                     integExternal :: IORef (Dynamics Double),
+                     integInternal :: IORef (Dynamics Double) }
+
+-- | Create a new integral with the specified initial value.
+newInteg :: Dynamics Double -> Dynamics Integ
+newInteg i = 
+  do r1 <- liftIO $ newIORef $ initD i 
+     r2 <- liftIO $ newIORef $ initD i 
+     let integ = Integ { integInit     = i, 
+                         integExternal = r1,
+                         integInternal = r2 }
+         z = Dynamics $ \p -> 
+           do (Dynamics m) <- readIORef (integInternal integ)
+              m p
+     y <- umemo z
+     liftIO $ writeIORef (integExternal integ) y
+     return integ
+
+-- | Return the integral's value.
+integValue :: Integ -> Dynamics Double
+integValue integ = 
+  Dynamics $ \p ->
+  do (Dynamics m) <- readIORef (integExternal integ)
+     m p
+
+-- | Set the derivative for the integral.
+integDiff :: Integ -> Dynamics Double -> Dynamics ()
+integDiff integ diff =
+  do let z = Dynamics $ \p ->
+           do y <- readIORef (integExternal integ)
+              let i = integInit integ
+              case spcMethod (pointSpecs p) of
+                Euler -> integEuler diff i y p
+                RungeKutta2 -> integRK2 diff i y p
+                RungeKutta4 -> integRK4 diff i y p
+     liftIO $ writeIORef (integInternal integ) z
+
+integEuler :: Dynamics Double
+             -> Dynamics Double 
+             -> Dynamics Double 
+             -> Point -> IO Double
+integEuler (Dynamics f) (Dynamics i) (Dynamics y) p = 
+  case pointIteration p of
+    0 -> 
+      i p
+    n -> do 
+      let sc = pointSpecs p
+          ty = basicTime sc (n - 1) 0
+          py = p { pointTime = ty, pointIteration = n - 1, pointPhase = 0 }
+      a <- y py
+      b <- f py
+      let !v = a + spcDT (pointSpecs p) * b
+      return v
+
+integRK2 :: Dynamics Double
+           -> Dynamics Double
+           -> Dynamics Double
+           -> Point -> IO Double
+integRK2 (Dynamics f) (Dynamics i) (Dynamics y) p =
+  case pointPhase p of
+    0 -> case pointIteration p of
+      0 ->
+        i p
+      n -> do
+        let sc = pointSpecs p
+            ty = basicTime sc (n - 1) 0
+            t1 = ty
+            t2 = basicTime sc (n - 1) 1
+            py = p { pointTime = ty, pointIteration = n - 1, pointPhase = 0 }
+            p1 = py
+            p2 = p { pointTime = t2, pointIteration = n - 1, pointPhase = 1 }
+        vy <- y py
+        k1 <- f p1
+        k2 <- f p2
+        let !v = vy + spcDT sc / 2.0 * (k1 + k2)
+        return v
+    1 -> do
+      let sc = pointSpecs p
+          n  = pointIteration p
+          ty = basicTime sc n 0
+          t1 = ty
+          py = p { pointTime = ty, pointIteration = n, pointPhase = 0 }
+          p1 = py
+      vy <- y py
+      k1 <- f p1
+      let !v = vy + spcDT sc * k1
+      return v
+    _ -> 
+      error "Incorrect phase: integRK2"
+
+integRK4 :: Dynamics Double
+           -> Dynamics Double
+           -> Dynamics Double
+           -> Point -> IO Double
+integRK4 (Dynamics f) (Dynamics i) (Dynamics y) p =
+  case pointPhase p of
+    0 -> case pointIteration p of
+      0 -> 
+        i p
+      n -> do
+        let sc = pointSpecs p
+            ty = basicTime sc (n - 1) 0
+            t1 = ty
+            t2 = basicTime sc (n - 1) 1
+            t3 = basicTime sc (n - 1) 2
+            t4 = basicTime sc (n - 1) 3
+            py = p { pointTime = ty, pointIteration = n - 1, pointPhase = 0 }
+            p1 = py
+            p2 = p { pointTime = t2, pointIteration = n - 1, pointPhase = 1 }
+            p3 = p { pointTime = t3, pointIteration = n - 1, pointPhase = 2 }
+            p4 = p { pointTime = t4, pointIteration = n - 1, pointPhase = 3 }
+        vy <- y py
+        k1 <- f p1
+        k2 <- f p2
+        k3 <- f p3
+        k4 <- f p4
+        let !v = vy + spcDT sc / 6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
+        return v
+    1 -> do
+      let sc = pointSpecs p
+          n  = pointIteration p
+          ty = basicTime sc n 0
+          t1 = ty
+          py = p { pointTime = ty, pointIteration = n, pointPhase = 0 }
+          p1 = py
+      vy <- y py
+      k1 <- f p1
+      let !v = vy + spcDT sc / 2.0 * k1
+      return v
+    2 -> do
+      let sc = pointSpecs p
+          n  = pointIteration p
+          ty = basicTime sc n 0
+          t2 = basicTime sc n 1
+          py = p { pointTime = ty, pointIteration = n, pointPhase = 0 }
+          p2 = p { pointTime = t2, pointIteration = n, pointPhase = 1 }
+      vy <- y py
+      k2 <- f p2
+      let !v = vy + spcDT sc / 2.0 * k2
+      return v
+    3 -> do
+      let sc = pointSpecs p
+          n  = pointIteration p
+          ty = basicTime sc n 0
+          t3 = basicTime sc n 2
+          py = p { pointTime = ty, pointIteration = n, pointPhase = 0 }
+          p3 = p { pointTime = t3, pointIteration = n, pointPhase = 2 }
+      vy <- y py
+      k3 <- f p3
+      let !v = vy + spcDT sc * k3
+      return v
+    _ -> 
+      error "Incorrect phase: integRK4"
+
+-- smoothI :: Dynamics Double -> Dynamics Double -> Dynamics Double 
+--           -> Dynamics Double
+-- smoothI x t i = y where
+--   y = integ ((x - y) / t) i
+
+-- smooth :: Dynamics Double -> Dynamics Double -> Dynamics Double
+-- smooth x t = smoothI x t x
+
+-- smooth3I :: Dynamics Double -> Dynamics Double -> Dynamics Double 
+--            -> Dynamics Double
+-- smooth3I x t i = y where
+--   y  = integ ((s1 - y) / t') i
+--   s1 = integ ((s0 - s1) / t') i
+--   s0 = integ ((x - s0) / t') i
+--   t' = t / 3.0
+
+-- smooth3 :: Dynamics Double -> Dynamics Double -> Dynamics Double
+-- smooth3 x t = smooth3I x t x
+
+-- smoothNI :: Dynamics Double -> Dynamics Double -> Int -> Dynamics Double 
+--            -> Dynamics Double
+-- smoothNI x t n i = s ! n where
+--   s   = array (1, n) [(k, f k) | k <- [1 .. n]]
+--   f 0 = integ ((x - s ! 0) / t') i
+--   f k = integ ((s ! (k - 1) - s ! k) / t') i
+--   t'  = t / fromIntegral n
+
+-- smoothN :: Dynamics Double -> Dynamics Double -> Int -> Dynamics Double
+-- smoothN x t n = smoothNI x t n x
+
+-- delay1I :: Dynamics Double -> Dynamics Double -> Dynamics Double 
+--           -> Dynamics Double
+-- delay1I x t i = y where
+--   y = integ (x - y) (i * t) / t
+
+-- delay1 :: Dynamics Double -> Dynamics Double -> Dynamics Double
+-- delay1 x t = delay1I x t x
+
+-- delay3I :: Dynamics Double -> Dynamics Double -> Dynamics Double 
+--           -> Dynamics Double
+-- delay3I x t i = y where
+--   y  = integ (s1 - y) (i * t') / t'
+--   s1 = integ (s0 - s1) (i * t') / t'
+--   s0 = integ (x - s0) (i * t') / t'
+--   t' = t / 3.0
+
+-- delay3 :: Dynamics Double -> Dynamics Double -> Dynamics Double
+-- delay3 x t = delay3I x t x
+
+-- delayNI :: Dynamics Double -> Dynamics Double -> Int -> Dynamics Double 
+--           -> Dynamics Double
+-- delayNI x t n i = s ! n where
+--   s   = array (1, n) [(k, f k) | k <- [1 .. n]]
+--   f 0 = integ (x - s ! 0) (i * t') / t'
+--   f k = integ (s ! (k - 1) - s ! k) (i * t') / t'
+--   t'  = t / fromIntegral n
+
+-- delayN :: Dynamics Double -> Dynamics Double -> Int -> Dynamics Double
+-- delayN x t n = delayNI x t n x
+
+-- forecast :: Dynamics Double -> Dynamics Double -> Dynamics Double 
+--            -> Dynamics Double
+-- forecast x at hz =
+--   x * (1.0 + (x / smooth x at - 1.0) / at * hz)
+
+-- trend :: Dynamics Double -> Dynamics Double -> Dynamics Double 
+--         -> Dynamics Double
+-- trend x at i =
+--   (x / smoothI x at (x / (1.0 + i * at)) - 1.0) / at
+
+--
+-- Integral Functions
+--
+
+-- | Return an integral with the specified derivative and initial value.
+-- If you want to create a loopback then you should use the 'Integ' type 
+-- directly. The 'integ' function is just a wrapper that uses this type.
+integ :: Dynamics Double -> Dynamics Double -> Dynamics (Dynamics Double)
+integ diff i =
+  do x <- newInteg i
+     integDiff x diff
+     return $ integValue x
+
+--
+-- Difference Equations
+--
+
+-- | The 'Sum' type represents a sum defined by some difference equation.
+data Sum a = Sum { sumInit     :: Dynamics a,   -- ^ The initial value.
+                   sumExternal :: IORef (Dynamics a),
+                   sumInternal :: IORef (Dynamics a) }
+
+-- | Create a new sum with the specified initial value.
+newSum :: (MArray IOUArray a IO, Num a) => Dynamics a -> Dynamics (Sum a)
+newSum i =   
+  do r1 <- liftIO $ newIORef $ initD i 
+     r2 <- liftIO $ newIORef $ initD i 
+     let sum = Sum { sumInit     = i, 
+                     sumExternal = r1,
+                     sumInternal = r2 }
+         z = Dynamics $ \p -> 
+           do (Dynamics m) <- readIORef (sumInternal sum)
+              m p
+     y <- umemo0 z
+     liftIO $ writeIORef (sumExternal sum) y
+     return sum
+
+-- | Return the total sum defined by the difference equation.
+sumValue :: Sum a -> Dynamics a
+sumValue sum = 
+  Dynamics $ \p ->
+  do (Dynamics m) <- readIORef (sumExternal sum)
+     m p
+
+-- | Set the difference equation for the sum.
+sumDiff :: (MArray IOUArray a IO, Num a) => Sum a -> Dynamics a -> Dynamics ()
+sumDiff sum (Dynamics diff) =
+  do let z = Dynamics $ \p ->
+           case pointIteration p of
+             0 -> do
+               let Dynamics i = sumInit sum
+               i p
+             n -> do 
+               Dynamics y <- readIORef (sumExternal sum)
+               let sc = pointSpecs p
+                   ty = basicTime sc (n - 1) 0
+                   py = p { pointTime = ty, 
+                            pointIteration = n - 1, 
+                            pointPhase = 0 }
+               a <- y py
+               b <- diff py
+               let !v = a + b
+               return v
+     liftIO $ writeIORef (sumInternal sum) z
+
+--
+-- Table Functions
+--
+
+-- | Lookup @x@ in a table of pairs @(x, y)@ using linear interpolation.
+lookupD :: Dynamics Double -> Array Int (Double, Double) -> Dynamics Double
+lookupD (Dynamics m) tbl =
+  Dynamics (\p -> do a <- m p; return $ find first last a) where
+    (first, last) = bounds tbl
+    find left right x =
+      if left > right then
+        error "Incorrect index: table"
+      else
+        let index = (left + 1 + right) `div` 2
+            x1    = fst $ tbl ! index
+        in if x1 <= x then 
+             let y | index < right = find index right x
+                   | right == last  = snd $ tbl ! right
+                   | otherwise     = 
+                     let x2 = fst $ tbl ! (index + 1)
+                         y1 = snd $ tbl ! index
+                         y2 = snd $ tbl ! (index + 1)
+                     in y1 + (y2 - y1) * (x - x1) / (x2 - x1) 
+             in y
+           else
+             let y | left < index = find left (index - 1) x
+                   | left == first = snd $ tbl ! left
+                   | otherwise    = error "Incorrect index: table"
+             in y
+
+-- | Lookup @x@ in a table of pairs @(x, y)@ using stepwise function.
+lookupStepwiseD :: Dynamics Double -> Array Int (Double, Double)
+                  -> Dynamics Double
+lookupStepwiseD (Dynamics m) tbl =
+  Dynamics (\p -> do a <- m p; return $ find first last a) where
+    (first, last) = bounds tbl
+    find left right x =
+      if left > right then
+        error "Incorrect index: table"
+      else
+        let index = (left + 1 + right) `div` 2
+            x1    = fst $ tbl ! index
+        in if x1 <= x then 
+             let y | index < right = find index right x
+                   | right == last  = snd $ tbl ! right
+                   | otherwise     = snd $ tbl ! right
+             in y
+           else
+             let y | left < index = find left (index - 1) x
+                   | left == first = snd $ tbl ! left
+                   | otherwise    = error "Incorrect index: table"
+             in y
+
+-- --
+-- -- Discrete Functions
+-- --
+    
+-- delayTrans :: Dynamics a -> Dynamics Double -> Dynamics a 
+--               -> (Dynamics a -> Dynamics a) -> Dynamics a
+-- delayTrans (Dynamics x) (Dynamics d) (Dynamics i) tr = tr $ Dynamics r 
+--   where
+--     r p = do 
+--       let t  = parTime p
+--           sc = parSpecs p
+--           n  = parIteration p
+--       a <- d p
+--       let t' = (t - a) - spcStartTime sc
+--           n' = fromInteger $ toInteger $ floor $ t' / spcDT sc
+--           y | n' < 0    = i $ p { pointTime = spcStartTime sc,
+--                                   pointIteration = 0, 
+--                                   pointPhase = 0 }
+--             | n' < n    = x $ p { pointTime = t',
+--                                   pointIteration = n',
+--                                   pointPhase = -1 }
+--             | n' > n    = error "Cannot return the future data: delay"
+--             | otherwise = error "Cannot return the current data: delay"
+--       y    
+
+-- delay :: (Memo a) => Dynamics a -> Dynamics Double -> Dynamics a
+-- delay x d = delayTrans x d x $ memo0 discrete
+
+-- delay' :: (UMemo a) => Dynamics a -> Dynamics Double -> Dynamics a
+-- delay' x d = delayTrans x d x $ memo0' discrete
+
+-- delayI :: (Memo a) => Dynamics a -> Dynamics Double -> Dynamics a -> Dynamics a
+-- delayI x d i = delayTrans x d i $ memo0 discrete         
+
+-- delayI' :: (UMemo a) => Dynamics a -> Dynamics Double -> Dynamics a -> Dynamics a
+-- delayI' x d i = delayTrans x d i $ memo0' discrete         
diff --git a/Simulation/Aivika/Dynamics/UVar.hs b/Simulation/Aivika/Dynamics/UVar.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Dynamics/UVar.hs
@@ -0,0 +1,126 @@
+
+{-# LANGUAGE FlexibleContexts #-}
+
+-- |
+-- Module     : Simulation.Aivika.Dynamics.UVar
+-- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.0.3
+--
+-- This module defines a variable that is bound to the event queue and 
+-- that keeps the history of changes storing the values in an unboxed array.
+--
+module Simulation.Aivika.Dynamics.UVar
+       (UVar,
+        newUVar,
+        uvarQueue,
+        readUVar,
+        writeUVar,
+        modifyUVar,
+        freezeUVar) where
+
+import Data.Array
+import Data.Array.IO
+import Data.IORef
+
+import Simulation.Aivika.Dynamics.Internal.Dynamics
+import Simulation.Aivika.Dynamics.EventQueue
+
+import qualified Simulation.Aivika.UVector as UV
+
+-- | A version of the 'Var' type which uses an unboxed array to store the values 
+-- in time points. You should prefer this type whenever possible.
+data UVar a = 
+  UVar { uvarQueue :: EventQueue, -- ^ Return the bound event queue.
+         uvarRun   :: Dynamics (),
+         uvarXS    :: UV.UVector Double, 
+         uvarYS    :: UV.UVector a}
+     
+-- | Create a new variable bound to the specified event queue.
+newUVar :: (MArray IOUArray a IO) => EventQueue -> a -> Dynamics (UVar a)
+newUVar q a =
+  Dynamics $ \p ->
+  do xs <- UV.newVector
+     ys <- UV.newVector
+     UV.appendVector xs $ spcStartTime $ pointSpecs p
+     UV.appendVector ys a
+     return UVar { uvarQueue = q,
+                   uvarRun   = queueRun q,
+                   uvarXS = xs,
+                   uvarYS = ys }
+
+-- | Read the value of a variable, forcing the bound event queue to raise 
+-- the events in case of need.
+readUVar :: (MArray IOUArray a IO) => UVar a -> Dynamics a
+readUVar v =
+  Dynamics $ \p ->
+  do let Dynamics m = uvarRun v
+     m p
+     let xs = uvarXS v
+         ys = uvarYS v
+         t  = pointTime p
+     count <- UV.vectorCount xs
+     let i = count - 1
+     x <- UV.readVector xs i
+     if x <= t 
+       then UV.readVector ys i
+       else do i <- UV.vectorBinarySearch xs t
+               if i >= 0
+                 then UV.readVector ys i
+                 else UV.readVector ys $ - (i + 1) - 1
+
+-- | Write a new value into the variable.
+writeUVar :: (MArray IOUArray a IO) => UVar a -> a -> Dynamics ()
+writeUVar v a =
+  Dynamics $ \p ->
+  do let xs = uvarXS v
+         ys = uvarYS v
+         t  = pointTime p
+     count <- UV.vectorCount xs
+     let i = count - 1
+     x <- UV.readVector xs i
+     if t < x 
+       then error "Cannot update the past data: writeUVar."
+       else if t == x
+            then UV.writeVector ys i $! a
+            else do UV.appendVector xs t
+                    UV.appendVector ys $! a
+
+-- | Mutate the contents of the variable, forcing the bound event queue to
+-- raise all pending events in case of need.
+modifyUVar :: (MArray IOUArray a IO) => UVar a -> (a -> a) -> Dynamics ()
+modifyUVar v f =
+  Dynamics $ \p ->
+  do let Dynamics m = uvarRun v
+     m p
+     let xs = uvarXS v
+         ys = uvarYS v
+         t  = pointTime p
+     count <- UV.vectorCount xs
+     let i = count - 1
+     x <- UV.readVector xs i
+     if t < x
+       then error "Cannot update the past data: modifyUVar."
+       else if t == x
+            then do a <- UV.readVector ys i
+                    UV.writeVector ys i $! f a
+            else do i <- UV.vectorBinarySearch xs t
+                    if i >= 0
+                      then do a <- UV.readVector ys i
+                              UV.appendVector xs t
+                              UV.appendVector ys $! f a
+                      else do a <- UV.readVector ys $ - (i + 1) - 1
+                              UV.appendVector xs t
+                              UV.appendVector ys $! f a
+
+-- | Freeze the variable and return in arrays the time points and corresponded 
+-- values when the variable had changed.
+freezeUVar :: (MArray IOUArray a IO) => 
+              UVar a -> Dynamics (Array Int Double, Array Int a)
+freezeUVar v =
+  Dynamics $ \p ->
+  do xs <- UV.freezeVector (uvarXS v)
+     ys <- UV.freezeVector (uvarYS v)
+     return (xs, ys)
diff --git a/Simulation/Aivika/Dynamics/Var.hs b/Simulation/Aivika/Dynamics/Var.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Dynamics/Var.hs
@@ -0,0 +1,126 @@
+
+-- |
+-- Module     : Simulation.Aivika.Dynamics.Var
+-- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.0.3
+--
+-- This module defines a variable that is bound to the event queue and 
+-- that keeps the history of changes storing the values in an array.
+--
+module Simulation.Aivika.Dynamics.Var
+       (Var,
+        newVar,
+        varQueue,
+        readVar,
+        writeVar,
+        modifyVar,
+        freezeVar) where
+
+import Data.Array
+import Data.Array.IO
+import Data.IORef
+
+import Simulation.Aivika.Dynamics.Internal.Dynamics
+import Simulation.Aivika.Dynamics.EventQueue
+
+import qualified Simulation.Aivika.Vector as V
+import qualified Simulation.Aivika.UVector as UV
+
+-- | Like the 'Ref' reference but keeps the history of changes in 
+-- different time points. The 'Var' variable is safe in the hybrid 
+-- simulation and when you use different event queues, but this variable is 
+-- slower than references.
+data Var a = 
+  Var { varQueue :: EventQueue,  -- ^ Return the bound event queue.
+        varRun   :: Dynamics (),
+        varXS    :: UV.UVector Double, 
+        varYS    :: V.Vector a}
+     
+-- | Create a new variable bound to the specified event queue.
+newVar :: EventQueue -> a -> Dynamics (Var a)
+newVar q a =
+  Dynamics $ \p ->
+  do xs <- UV.newVector
+     ys <- V.newVector
+     UV.appendVector xs $ spcStartTime $ pointSpecs p
+     V.appendVector ys a
+     return Var { varQueue = q,
+                  varRun   = queueRun q,
+                  varXS = xs,
+                  varYS = ys }
+
+-- | Read the value of a variable, forcing the bound event queue to raise 
+-- the events in case of need.
+readVar :: Var a -> Dynamics a
+readVar v =
+  Dynamics $ \p ->
+  do let Dynamics m = varRun v
+     m p
+     let xs = varXS v
+         ys = varYS v
+         t  = pointTime p
+     count <- UV.vectorCount xs
+     let i = count - 1
+     x <- UV.readVector xs i
+     if x <= t 
+       then V.readVector ys i
+       else do i <- UV.vectorBinarySearch xs t
+               if i >= 0
+                 then V.readVector ys i
+                 else V.readVector ys $ - (i + 1) - 1
+
+-- | Write a new value into the variable.
+writeVar :: Var a -> a -> Dynamics ()
+writeVar v a =
+  Dynamics $ \p ->
+  do let xs = varXS v
+         ys = varYS v
+         t  = pointTime p
+     count <- UV.vectorCount xs
+     let i = count - 1
+     x <- UV.readVector xs i
+     if t < x 
+       then error "Cannot update the past data: writeVar."
+       else if t == x
+            then V.writeVector ys i $! a
+            else do UV.appendVector xs t
+                    V.appendVector ys $! a
+
+-- | Mutate the contents of the variable, forcing the bound event queue to
+-- raise all pending events in case of need.
+modifyVar :: Var a -> (a -> a) -> Dynamics ()
+modifyVar v f =
+  Dynamics $ \p ->
+  do let Dynamics m = varRun v
+     m p
+     let xs = varXS v
+         ys = varYS v
+         t  = pointTime p
+     count <- UV.vectorCount xs
+     let i = count - 1
+     x <- UV.readVector xs i
+     if t < x
+       then error "Cannot update the past data: modifyVar."
+       else if t == x
+            then do a <- V.readVector ys i
+                    V.writeVector ys i $! f a
+            else do i <- UV.vectorBinarySearch xs t
+                    if i >= 0
+                      then do a <- V.readVector ys i
+                              UV.appendVector xs t
+                              V.appendVector ys $! f a
+                      else do a <- V.readVector ys $ - (i + 1) - 1
+                              UV.appendVector xs t
+                              V.appendVector ys $! f a
+
+-- | Freeze the variable and return in arrays the time points and corresponded 
+-- values when the variable had changed.
+freezeVar :: Var a -> Dynamics (Array Int Double, Array Int a)
+freezeVar v =
+  Dynamics $ \p ->
+  do xs <- UV.freezeVector (varXS v)
+     ys <- V.freezeVector (varYS v)
+     return (xs, ys)
diff --git a/Simulation/Aivika/PriorityQueue.hs b/Simulation/Aivika/PriorityQueue.hs
--- a/Simulation/Aivika/PriorityQueue.hs
+++ b/Simulation/Aivika/PriorityQueue.hs
@@ -1,40 +1,18 @@
 
--- Copyright (c) 2009, 2010, 2011 David Sorokin <david.sorokin@gmail.com>
--- 
--- All rights reserved.
--- 
--- Redistribution and use in source and binary forms, with or without
--- modification, are permitted provided that the following conditions
--- are met:
--- 
--- 1. Redistributions of source code must retain the above copyright
---    notice, this list of conditions and the following disclaimer.
--- 
--- 2. Redistributions in binary form must reproduce the above copyright
---    notice, this list of conditions and the following disclaimer in the
---    documentation and/or other materials provided with the distribution.
--- 
--- 3. Neither the name of the author nor the names of his contributors
---    may be used to endorse or promote products derived from this software
---    without specific prior written permission.
--- 
--- THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND
--- ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
--- IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
--- ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE
--- FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
--- DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
--- OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
--- HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
--- LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
--- OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
--- SUCH DAMAGE.
-
--- This is an imperative version of the heap-based priority queue in monad IO.
-
+-- |
+-- Module     : Simulation.Aivika.PriorityQueue
+-- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.0.3
+--
+-- An imperative heap-based priority queue.
+--
 module Simulation.Aivika.PriorityQueue 
        (PriorityQueue, 
         queueNull, 
+        queueCount,
         newQueue, 
         enqueue, 
         dequeue, 
@@ -46,12 +24,12 @@
 import Data.IORef
 import Control.Monad
 
+-- | The 'PriorityQueue' type represents an imperative heap-based 
+-- priority queue.
 data PriorityQueue a = 
   PriorityQueue { pqKeys  :: IORef (IOUArray Int Double),
                   pqVals  :: IORef (IOArray Int a),
-                  pqSize  :: IORef Int,
-                  pqNoVal :: a    -- to release references                 
-                }
+                  pqSize  :: IORef Int }
 
 increase :: PriorityQueue a -> Int -> IO ()
 increase pq capacity = 
@@ -120,13 +98,19 @@
                               writeArray vals i vn''
                               siftDown keys vals size n'' k v
 
+-- | Test whether the priority queue is empty.
 queueNull :: PriorityQueue a -> IO Bool
 queueNull pq =
   do size <- readIORef (pqSize pq)
      return $ size == 0
 
-newQueue :: a -> IO (PriorityQueue a)
-newQueue defaultValue =
+-- | Return the number of elements in the priority queue.
+queueCount :: PriorityQueue a -> IO Int
+queueCount pq = readIORef (pqSize pq)
+
+-- | Create a new priority queue.
+newQueue :: IO (PriorityQueue a)
+newQueue =
   do keys <- newArray_ (0, 10)
      vals <- newArray_ (0, 10)
      keyRef  <- newIORef keys
@@ -134,20 +118,21 @@
      sizeRef <- newIORef 0
      return PriorityQueue { pqKeys = keyRef, 
                             pqVals = valRef, 
-                            pqSize = sizeRef,
-                            pqNoVal = defaultValue }
+                            pqSize = sizeRef }
 
+-- | Enqueue a new element with the specified priority.
 enqueue :: PriorityQueue a -> Double -> a -> IO ()
 enqueue pq k v =
   do i <- readIORef (pqSize pq)
      keys <- readIORef (pqKeys pq)
      (il, iu) <- getBounds keys
-     when (i >= iu - il + 1) $ increase pq (i + 1)
+     when (i >= iu - il) $ increase pq (i + 2)  -- plus one element on the end
      writeIORef (pqSize pq) (i + 1)
      keys <- readIORef (pqKeys pq)  -- it can be another! (side-effect)
      vals <- readIORef (pqVals pq)
      siftUp keys vals i k v
 
+-- | Dequeue the element with the minimal priority.
 dequeue :: PriorityQueue a -> IO ()
 dequeue pq =
   do size <- readIORef (pqSize pq)
@@ -156,12 +141,15 @@
      writeIORef (pqSize pq) i
      keys <- readIORef (pqKeys pq)
      vals <- readIORef (pqVals pq)
-     k <- readArray keys i
-     v <- readArray vals i
-     writeArray keys i 0.0
-     writeArray vals i (pqNoVal pq)    -- to release the reference!
+     k  <- readArray keys i
+     v  <- readArray vals i
+     k0 <- readArray keys size
+     v0 <- readArray vals size
+     writeArray keys i k0
+     writeArray vals i v0
      siftDown keys vals i 0 k v
 
+-- | Return the element with the minimal priority.
 queueFront :: PriorityQueue a -> IO (Double, a)
 queueFront pq =
   do size <- readIORef (pqSize pq)
diff --git a/Simulation/Aivika/Queue.hs b/Simulation/Aivika/Queue.hs
--- a/Simulation/Aivika/Queue.hs
+++ b/Simulation/Aivika/Queue.hs
@@ -1,40 +1,18 @@
 
--- Copyright (c) 2009, 2010, 2011 David Sorokin <david.sorokin@gmail.com>
--- 
--- All rights reserved.
--- 
--- Redistribution and use in source and binary forms, with or without
--- modification, are permitted provided that the following conditions
--- are met:
--- 
--- 1. Redistributions of source code must retain the above copyright
---    notice, this list of conditions and the following disclaimer.
--- 
--- 2. Redistributions in binary form must reproduce the above copyright
---    notice, this list of conditions and the following disclaimer in the
---    documentation and/or other materials provided with the distribution.
--- 
--- 3. Neither the name of the author nor the names of his contributors
---    may be used to endorse or promote products derived from this software
---    without specific prior written permission.
--- 
--- THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND
--- ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
--- IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
--- ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE
--- FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
--- DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
--- OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
--- HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
--- LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
--- OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
--- SUCH DAMAGE.
-
--- This is an imperative version of the queue in monad IO.
-
+-- |
+-- Module     : Simulation.Aivika.Queue
+-- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.0.3
+--
+-- An imperative double-linked queue.
+--
 module Simulation.Aivika.Queue 
        (Queue, 
         queueNull, 
+        queueCount,
         newQueue, 
         enqueue, 
         dequeue, 
@@ -43,15 +21,19 @@
 import Data.IORef
 import Control.Monad
 
+-- | A cell of the double-linked queue.
 data QueueItem a = 
   QueueItem { qiVal  :: a,
               qiPrev :: IORef (Maybe (QueueItem a)),
               qiNext :: IORef (Maybe (QueueItem a)) }
   
+-- | The 'Queue' type represents an imperative double-linked queue.
 data Queue a =  
   Queue { qHead :: IORef (Maybe (QueueItem a)),
-          qTail :: IORef (Maybe (QueueItem a)) }
+          qTail :: IORef (Maybe (QueueItem a)), 
+          qSize :: IORef Int }
 
+-- | Test whether the queue is empty.
 queueNull :: Queue a -> IO Bool
 queueNull q =
   do head <- readIORef (qHead q) 
@@ -59,15 +41,24 @@
        Nothing -> return True
        Just _  -> return False
     
+-- | Return the number of elements in the queue.
+queueCount :: Queue a -> IO Int
+queueCount q = readIORef (qSize q)
+
+-- | Create a new queue.
 newQueue :: IO (Queue a)
 newQueue =
   do head <- newIORef Nothing 
      tail <- newIORef Nothing
-     return Queue { qHead = head, qTail = tail }
+     size <- newIORef 0
+     return Queue { qHead = head, qTail = tail, qSize = size }
 
+-- | Enqueue a new element.
 enqueue :: Queue a -> a -> IO ()
 enqueue q v =
-  do head <- readIORef (qHead q)
+  do size <- readIORef (qSize q)
+     writeIORef (qSize q) (size + 1)
+     head <- readIORef (qHead q)
      case head of
        Nothing ->
          do prev <- newIORef Nothing
@@ -86,6 +77,7 @@
             writeIORef (qiPrev h) item
             writeIORef (qHead q) item
 
+-- | Dequeue the first element.
 dequeue :: Queue a -> IO ()
 dequeue q =
   do tail <- readIORef (qTail q) 
@@ -93,7 +85,9 @@
        Nothing ->
          error "Empty queue: dequeue"
        Just t ->
-         do tail' <- readIORef (qiPrev t)
+         do size  <- readIORef (qSize q)
+            writeIORef (qSize q) (size - 1)
+            tail' <- readIORef (qiPrev t)
             case tail' of
               Nothing ->
                 do writeIORef (qHead q) Nothing
@@ -102,6 +96,7 @@
                 do writeIORef (qiNext t') Nothing
                    writeIORef (qTail q) tail'
 
+-- | Return the first element.
 queueFront :: Queue a -> IO a
 queueFront q =
   do tail <- readIORef (qTail q)
diff --git a/Simulation/Aivika/Statistics.hs b/Simulation/Aivika/Statistics.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Statistics.hs
@@ -0,0 +1,118 @@
+
+{-# LANGUAGE FlexibleContexts #-}
+
+-- |
+-- Module     : Simulation.Aivika.Statistics
+-- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.0.3
+--
+-- Represents statistics and results.
+--
+module Simulation.Aivika.Statistics
+       (Statistics, 
+        newStatistics,
+        addStatistics,
+        statisticsData,
+        AnalysisResults(..),
+        analyzeData,
+        showResults) where 
+
+import Data.Foldable
+import Data.Array
+import Data.Array.IO
+import Control.Monad
+import Control.Concurrent.MVar
+
+import Simulation.Aivika.UVector
+
+-- | Represents statistics. 
+-- 
+-- All functions with the statistics in this module are thread-safe. Therefore 
+-- you can use them in experiments when parallel simulations execute simultaneously.
+data Statistics a = Statistics { statData :: UVector a, 
+                                 statLock :: MVar () }
+
+-- | Create new statistics.
+newStatistics :: (MArray IOUArray a IO) => IO (Statistics a)
+newStatistics = 
+  do v <- newVector
+     l <- newMVar ()
+     return Statistics { statData = v, 
+                         statLock = l }
+
+-- | Add data to the statistics. It is thread-safe.
+addStatistics :: (MArray IOUArray a IO) => Statistics a -> a -> IO ()
+addStatistics s x = 
+  withMVar (statLock s) $ \() ->
+  appendVector (statData s) x
+
+-- | Return the statistics data. It is thread-safe.
+statisticsData :: (MArray IOUArray a IO) => Statistics a -> IO (Array Int a)
+statisticsData s =
+  withMVar (statLock s) $ \() -> freezeVector (statData s)
+       
+-- | Represents the results of the statistic analysis.
+data AnalysisResults a = 
+  AnalysisResults { resultsData     :: Array Int a,
+                    -- ^ Statistic data.
+                    resultsMean     :: Double,
+                    -- ^ The average value.
+                    resultsVariance :: Double,
+                    -- ^ The variance.
+                    resultsMin      :: a,
+                    -- ^ The minimum value.
+                    resultsMax      :: a 
+                    -- ^ The maximum value.
+                  } deriving (Eq, Ord, Show)
+
+-- | Analyze data.
+analyzeData :: Real a => Array Int a -> AnalysisResults a
+analyzeData xs =
+  let (i1, i2) = bounds xs
+      meanx = foldl' (\y i -> y * (1 - k i) + f i * k i) 0 [i1 .. i2]
+      sqrx  = foldl' (\y i -> y * (1 - k i) + g i * k i) 0 [i1 .. i2]
+      minx  = foldl' (\y i -> if i == 0 then x i else min y (x i)) 0 [i1 .. i2]
+      maxx  = foldl' (\y i -> if i == 0 then x i else max y (x i)) 0 [i1 .. i2]
+      x i = xs ! i
+      f i = fromRational (toRational (x i))
+      g i = let y = f i in y * y
+      k i = 1 / fromInteger (toInteger (i - i1 + 1))
+  in AnalysisResults { resultsData = xs,
+                       resultsMean = meanx,
+                       resultsVariance = sqrx - meanx * meanx,
+                       resultsMin = minx,
+                       resultsMax = maxx }
+       
+-- | Show the results of analysis with the specified indent.       
+showResults :: (Show a) => AnalysisResults a -> Int -> ShowS
+showResults rs indent =
+  let (i1, i2) = bounds (resultsData rs)
+      tab = replicate indent ' '
+  in if i1 <= i2
+     then
+       showString tab .
+       showString "mean      = " . shows (resultsMean rs) . 
+       showString "\n" . 
+       showString tab .
+       showString "deviation = " . shows (sqrt (resultsVariance rs)) . 
+       showString "\n" .
+       showString tab .
+       showString "minimum   = " . shows (resultsMin rs) . 
+       showString "\n" .
+       showString tab .
+       showString "maximum   = " . shows (resultsMax rs)
+     else
+       showString tab .
+       showString "mean      = ---" .
+       showString "\n" . 
+       showString tab .
+       showString "deviation = ---" .
+       showString "\n" . 
+       showString tab .
+       showString "minimum   = ---" .
+       showString "\n" . 
+       showString tab .
+       showString "maximum   = ---"
diff --git a/Simulation/Aivika/UVector.hs b/Simulation/Aivika/UVector.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/UVector.hs
@@ -0,0 +1,130 @@
+
+{-# LANGUAGE FlexibleContexts #-}
+
+-- |
+-- Module     : Simulation.Aivika.UVector
+-- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.0.3
+--
+-- An imperative unboxed vector.
+--
+module Simulation.Aivika.UVector
+       (UVector, 
+        newVector, 
+        copyVector, 
+        vectorCount, 
+        appendVector, 
+        readVector, 
+        writeVector, 
+        vectorBinarySearch,
+        freezeVector) where 
+
+import Data.Array
+import Data.Array.MArray
+import Data.Array.IO
+import Data.IORef
+import Control.Monad
+
+-- | Represents an unboxed resizable vector.
+data UVector a = UVector { vectorArrayRef :: IORef (IOUArray Int a),
+                           vectorCountRef :: IORef Int, 
+                           vectorCapacityRef :: IORef Int }
+
+-- | Create a new vector.
+newVector :: MArray IOUArray a IO => IO (UVector a)
+newVector = 
+  do array <- newArray_ (0, 4 - 1)
+     arrayRef <- newIORef array
+     countRef <- newIORef 0
+     capacityRef <- newIORef 4
+     return UVector { vectorArrayRef = arrayRef,
+                      vectorCountRef = countRef,
+                      vectorCapacityRef = capacityRef }
+
+-- | Copy the vector.
+copyVector :: (MArray IOUArray a IO) => UVector a -> IO (UVector a)
+copyVector vector =
+  do array <- readIORef (vectorArrayRef vector)
+     count <- readIORef (vectorCountRef vector)
+     array' <- newArray_ (0, count - 1)
+     arrayRef' <- newIORef array'
+     countRef' <- newIORef count
+     capacityRef' <- newIORef count
+     forM_ [0 .. count - 1] $ \i ->
+       do x <- readArray array i
+          writeArray array' i x
+     return UVector { vectorArrayRef = arrayRef',
+                      vectorCountRef = countRef',
+                      vectorCapacityRef = capacityRef' }
+
+-- | Ensure that the vector has the specified capacity.
+vectorEnsureCapacity :: MArray IOUArray a IO => UVector a -> Int -> IO ()
+vectorEnsureCapacity vector capacity =
+  do capacity' <- readIORef (vectorCapacityRef vector)
+     when (capacity' < capacity) $
+       do array' <- readIORef (vectorArrayRef vector)
+          count' <- readIORef (vectorCountRef vector)
+          let capacity'' = max (2 * capacity') capacity
+          array'' <- newArray_ (0, capacity'' - 1)
+          forM_ [0 .. count' - 1] $ \i ->
+            do x <- readArray array' i
+               writeArray array'' i x
+          writeIORef (vectorArrayRef vector) array''
+          writeIORef (vectorCapacityRef vector) capacity''
+          
+-- | Return the element count.
+vectorCount :: MArray IOUArray a IO => UVector a -> IO Int
+vectorCount vector = readIORef (vectorCountRef vector)
+          
+-- | Add the specified element to the end of the vector.
+appendVector :: MArray IOUArray a IO => UVector a -> a -> IO ()          
+appendVector vector item =
+  do count <- readIORef (vectorCountRef vector)
+     vectorEnsureCapacity vector (count + 1)
+     array <- readIORef (vectorArrayRef vector)
+     writeArray array count item
+     writeIORef (vectorCountRef vector) (count + 1)
+     
+-- | Read a value from the vector, where indices are started from 0.
+readVector :: MArray IOUArray a IO => UVector a -> Int -> IO a
+readVector vector index =
+  do array <- readIORef (vectorArrayRef vector)
+     readArray array index
+          
+-- | Set an array item at the specified index which is started from 0.
+writeVector :: MArray IOUArray a IO => UVector a -> Int -> a -> IO ()
+writeVector vector index item =
+  do array <- readIORef (vectorArrayRef vector)
+     writeArray array index item
+          
+vectorBinarySearch' :: (MArray IOUArray a IO, Ord a) => 
+                      IOUArray Int a -> a -> Int -> Int -> IO Int
+vectorBinarySearch' array item left right =
+  if left > right 
+  then return $ - (right + 1) - 1
+  else
+    do let index = (left + right) `div` 2
+       curr <- readArray array index
+       if item < curr 
+         then vectorBinarySearch' array item left (index - 1)
+         else if item == curr
+              then return index
+              else vectorBinarySearch' array item (index + 1) right
+                   
+-- | Return the index of the specified element using binary search; otherwise, 
+-- a negated insertion index minus one: 0 -> -0 - 1, ..., i -> -i - 1, ....
+vectorBinarySearch :: (MArray IOUArray a IO, Ord a) => UVector a -> a -> IO Int
+vectorBinarySearch vector item =
+  do array <- readIORef (vectorArrayRef vector)
+     count <- readIORef (vectorCountRef vector)
+     vectorBinarySearch' array item 0 (count - 1)
+
+freezeVector :: (MArray IOUArray a IO) => UVector a -> IO (Array Int a)
+freezeVector vector = 
+  do vector' <- copyVector vector
+     array   <- readIORef (vectorArrayRef vector')
+     freeze array
+     
diff --git a/Simulation/Aivika/Vector.hs b/Simulation/Aivika/Vector.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Vector.hs
@@ -0,0 +1,127 @@
+
+-- |
+-- Module     : Simulation.Aivika.Vector
+-- Copyright  : Copyright (c) 2009-2011, David Sorokin <david.sorokin@gmail.com>
+-- License    : BSD3
+-- Maintainer : David Sorokin <david.sorokin@gmail.com>
+-- Stability  : experimental
+-- Tested with: GHC 7.0.3
+--
+-- An imperative vector.
+--
+module Simulation.Aivika.Vector
+       (Vector, 
+        newVector, 
+        copyVector,
+        vectorCount, 
+        appendVector, 
+        readVector, 
+        writeVector,
+        vectorBinarySearch,
+        freezeVector) where 
+
+import Data.Array
+import Data.Array.MArray
+import Data.Array.IO
+import Data.IORef
+import Control.Monad
+
+-- | Represents a resizable vector.
+data Vector a = Vector { vectorArrayRef :: IORef (IOArray Int a),
+                         vectorCountRef :: IORef Int, 
+                         vectorCapacityRef :: IORef Int }
+
+-- | Create a new vector.
+newVector :: IO (Vector a)
+newVector = 
+  do array <- newArray_ (0, 4 - 1)
+     arrayRef <- newIORef array
+     countRef <- newIORef 0
+     capacityRef <- newIORef 4
+     return Vector { vectorArrayRef = arrayRef,
+                     vectorCountRef = countRef,
+                     vectorCapacityRef = capacityRef }
+
+-- | Copy the vector.
+copyVector :: Vector a -> IO (Vector a)
+copyVector vector =
+  do array <- readIORef (vectorArrayRef vector)
+     count <- readIORef (vectorCountRef vector)
+     array' <- newArray_ (0, count - 1)
+     arrayRef' <- newIORef array'
+     countRef' <- newIORef count
+     capacityRef' <- newIORef count
+     forM_ [0 .. count - 1] $ \i ->
+       do x <- readArray array i
+          writeArray array' i x
+     return Vector { vectorArrayRef = arrayRef',
+                     vectorCountRef = countRef',
+                     vectorCapacityRef = capacityRef' }
+
+-- | Ensure that the vector has the specified capacity.
+vectorEnsureCapacity :: Vector a -> Int -> IO ()
+vectorEnsureCapacity vector capacity =
+  do capacity' <- readIORef (vectorCapacityRef vector)
+     when (capacity' < capacity) $
+       do array' <- readIORef (vectorArrayRef vector)
+          count' <- readIORef (vectorCountRef vector)
+          let capacity'' = max (2 * capacity') capacity
+          array'' <- newArray_ (0, capacity'' - 1)
+          forM_ [0 .. count' - 1] $ \i ->
+            do x <- readArray array' i
+               writeArray array'' i x
+          writeIORef (vectorArrayRef vector) array''
+          writeIORef (vectorCapacityRef vector) capacity''
+          
+-- | Return the element count.
+vectorCount :: Vector a -> IO Int
+vectorCount vector = readIORef (vectorCountRef vector)
+          
+-- | Add the specified element to the end of the vector.
+appendVector :: Vector a -> a -> IO ()          
+appendVector vector item =
+  do count <- readIORef (vectorCountRef vector)
+     vectorEnsureCapacity vector (count + 1)
+     array <- readIORef (vectorArrayRef vector)
+     writeArray array count item
+     writeIORef (vectorCountRef vector) (count + 1)
+     
+-- | Read a value from the vector, where indices are started from 0.
+readVector :: Vector a -> Int -> IO a
+readVector vector index =
+  do array <- readIORef (vectorArrayRef vector)
+     readArray array index
+          
+-- | Set an array item at the specified index which is started from 0.
+writeVector :: Vector a -> Int -> a -> IO ()
+writeVector vector index item =
+  do array <- readIORef (vectorArrayRef vector)
+     writeArray array index item
+
+vectorBinarySearch' :: Ord a => IOArray Int a -> a -> Int -> Int -> IO Int
+vectorBinarySearch' array item left right =
+  if left > right 
+  then return $ - (right + 1) - 1
+  else
+    do let index = (left + right) `div` 2
+       curr <- readArray array index
+       if item < curr 
+         then vectorBinarySearch' array item left (index - 1)
+         else if item == curr
+              then return index
+              else vectorBinarySearch' array item (index + 1) right
+                   
+-- | Return the index of the specified element using binary search; otherwise, 
+-- a negated insertion index minus one: 0 -> -0 - 1, ..., i -> -i - 1, ....
+vectorBinarySearch :: Ord a => Vector a -> a -> IO Int
+vectorBinarySearch vector item =
+  do array <- readIORef (vectorArrayRef vector)
+     count <- readIORef (vectorCountRef vector)
+     vectorBinarySearch' array item 0 (count - 1)
+
+freezeVector :: Vector a -> IO (Array Int a)
+freezeVector vector = 
+  do vector' <- copyVector vector
+     array   <- readIORef (vectorArrayRef vector')
+     freeze array
+     
diff --git a/aivika.cabal b/aivika.cabal
--- a/aivika.cabal
+++ b/aivika.cabal
@@ -1,28 +1,34 @@
 name:            aivika
-version:         0.1
+version:         0.2
 synopsis:        A multi-paradigm simulation library
 description:
-    Aivika is a multi-paradigm simulation library. It allows us to integrate 
-    a system of ordinary differential equations. Also it can be applied to
-    the Discrete Event Simulation. It supports the event-oriented, 
-    process-oriented and activity-oriented paradigms. Aivika also supports 
-    the Agent-based Modeling. Finally, it can be applied to System Dynamics.
+    Aivika is a small simulation library that covers many paradigms. 
+    It allows integrating a system of ordinary differential equations. 
+    Also it can be applied to the Discrete Event Simulation. It supports 
+    the event-oriented, process-oriented and activity-oriented paradigms. 
+    Aivika also supports the Agent-based Modeling. Finally, it can be applied 
+    to System Dynamics.
     .
     The library widely uses monads. The dynamic system is represented as 
-    a computation in the Dynamics monad. There is also the DynamicsProc
-    monad to represent the discontinuous processes which can be suspended
-    at any time and then resumed later. Everything else is expressed through 
-    these two monads, including the events, agent handlers and even integrals.
+    a computation in the Dynamics monad. There is also the Process
+    monad to represent the discontinuous processes which can suspend
+    at any time and then resume later. Everything else is expressed through 
+    these two monads, including the event handlers, agent handlers and even 
+    integrals.
     .
+    The PDF documentation is available at 
+    <https://github.com/dsorokin/aivika/blob/master/doc/aivika.pdf>
+    .
 category:        Simulation
 license:         BSD3
 license-file:    LICENSE
 copyright:       (c) 2009-2011. David Sorokin <david.sorokin@gmail.com>
 author:          David Sorokin
 maintainer:      David Sorokin <david.sorokin@gmail.com>
+homepage:        http://github.com/dsorokin/aivika
 cabal-version:   >= 1.2.0
 build-type:      Simple
-tested-with:     GHC == 6.12.1
+tested-with:     GHC == 7.0.3
 
 extra-source-files:  examples/BassDiffusion.hs
                      examples/ChemicalReaction.hs
@@ -32,16 +38,46 @@
                      examples/MachRep1TimeDriven.hs
                      examples/MachRep2.hs
                      examples/MachRep3.hs
+                     examples/Furnace.hs
 
 data-files:          doc/aivika.pdf
 
 library
+
     exposed-modules: Simulation.Aivika.Dynamics
-    other-modules:   Simulation.Aivika.Queue
+                     Simulation.Aivika.Dynamics.Agent
+                     Simulation.Aivika.Dynamics.Base
+                     Simulation.Aivika.Dynamics.Cont
+                     Simulation.Aivika.Dynamics.EventQueue
+                     Simulation.Aivika.Dynamics.Lift
+                     Simulation.Aivika.Dynamics.Process
+                     Simulation.Aivika.Dynamics.Random
+                     Simulation.Aivika.Dynamics.Ref
+                     Simulation.Aivika.Dynamics.Resource
+                     Simulation.Aivika.Dynamics.SystemDynamics
+                     Simulation.Aivika.Dynamics.UVar
+                     Simulation.Aivika.Dynamics.Var
+                     Simulation.Aivika.Dynamics.Parameter
                      Simulation.Aivika.PriorityQueue
+                     Simulation.Aivika.Queue
+                     Simulation.Aivika.Statistics
+
+    other-modules:   Simulation.Aivika.Dynamics.Internal.Dynamics
+                     Simulation.Aivika.Dynamics.Internal.Cont
+                     Simulation.Aivika.Dynamics.Internal.Process
+                     Simulation.Aivika.Dynamics.Internal.Time
+                     Simulation.Aivika.Dynamics.Internal.Memo
+                     Simulation.Aivika.Dynamics.Internal.Interpolate
+                     Simulation.Aivika.Dynamics.Internal.Fold
+                     Simulation.Aivika.Vector
+                     Simulation.Aivika.UVector
                      
-    build-depends:   base >= 3 && < 5,
+    build-depends:   base >= 3 && < 6,
+                     haskell98,
                      mtl >= 1.1.0.2,
-                     array >= 0.3.0.0
+                     array >= 0.3.0.0,
+                     containers >= 0.4.0.0
+
+    extensions:      FlexibleContexts
                      
     ghc-options:     -O2
diff --git a/doc/aivika.pdf b/doc/aivika.pdf
Binary files a/doc/aivika.pdf and b/doc/aivika.pdf differ
diff --git a/examples/BassDiffusion.hs b/examples/BassDiffusion.hs
--- a/examples/BassDiffusion.hs
+++ b/examples/BassDiffusion.hs
@@ -5,6 +5,9 @@
 import Control.Monad.Trans
 
 import Simulation.Aivika.Dynamics
+import Simulation.Aivika.Dynamics.EventQueue
+import Simulation.Aivika.Dynamics.Agent
+import Simulation.Aivika.Dynamics.Ref
 
 n = 500    -- the number of agents
 
@@ -31,7 +34,7 @@
                        personPotentialAdopter :: AgentState,
                        personAdopter :: AgentState }
               
-createPerson :: DynamicsQueue -> Dynamics Person              
+createPerson :: EventQueue -> Dynamics Person              
 createPerson q =    
   do agent <- newAgent q
      potentialAdopter <- newState agent
@@ -40,29 +43,27 @@
                      personPotentialAdopter = potentialAdopter,
                      personAdopter = adopter }
        
-createPersons :: DynamicsQueue -> Dynamics (Array Int Person)
+createPersons :: EventQueue -> Dynamics (Array Int Person)
 createPersons q =
   do list <- forM [1 .. n] $ \i ->
        do p <- createPerson q
           return (i, p)
      return $ array (1, n) list
      
-definePerson :: Person -> Array Int Person 
-               -> DynamicsRef Int -> DynamicsRef Int
-               -> Dynamics ()
+definePerson :: Person -> Array Int Person -> Ref Int -> Ref Int -> Dynamics ()
 definePerson p ps potentialAdopters adopters =
   do stateActivation (personPotentialAdopter p) $
-       do modifyRef' potentialAdopters $ \a -> a + 1
+       do modifyRef potentialAdopters $ \a -> a + 1
           -- add a timeout
           t <- liftIO $ exprnd advertisingEffectiveness 
           let st  = personPotentialAdopter p
               st' = personAdopter p
           addTimeout st t $ activateState st'
      stateActivation (personAdopter p) $ 
-       do modifyRef' adopters  $ \a -> a + 1
+       do modifyRef adopters  $ \a -> a + 1
           -- add a timer that works while the state is active
           let t = liftIO $ exprnd contactRate    -- many times!
-          addTimerD (personAdopter p) t $
+          addTimer (personAdopter p) t $
             do i <- liftIO $ getStdRandom $ randomR (1, n)
                let p' = ps ! i
                st <- agentState (personAgent p')
@@ -70,14 +71,11 @@
                  do b <- liftIO $ boolrnd adoptionFraction
                     when b $ activateState (personAdopter p')
      stateDeactivation (personPotentialAdopter p) $
-       modifyRef' potentialAdopters $ \a -> a - 1
+       modifyRef potentialAdopters $ \a -> a - 1
      stateDeactivation (personAdopter p) $
-       modifyRef' adopters $ \a -> a - 1
+       modifyRef adopters $ \a -> a - 1
         
-definePersons :: Array Int Person 
-                -> DynamicsRef Int 
-                -> DynamicsRef Int 
-                -> Dynamics ()
+definePersons :: Array Int Person -> Ref Int -> Ref Int -> Dynamics ()
 definePersons ps potentialAdopters adopters =
   forM_ (elems ps) $ \p -> 
   definePerson p ps potentialAdopters adopters
@@ -102,5 +100,4 @@
                  return [i1, i2]
 
 main =
-  do xs <- runDynamics model specs
-     print xs
+  printDynamics model specs
diff --git a/examples/ChemicalReaction.hs b/examples/ChemicalReaction.hs
--- a/examples/ChemicalReaction.hs
+++ b/examples/ChemicalReaction.hs
@@ -1,5 +1,6 @@
 
 import Simulation.Aivika.Dynamics
+import Simulation.Aivika.Dynamics.SystemDynamics
 
 specs = Specs { spcStartTime = 0, 
                 spcStopTime = 13, 
diff --git a/examples/FishBank.hs b/examples/FishBank.hs
--- a/examples/FishBank.hs
+++ b/examples/FishBank.hs
@@ -2,6 +2,7 @@
 import Data.Array
 
 import Simulation.Aivika.Dynamics
+import Simulation.Aivika.Dynamics.SystemDynamics
 
 specs = Specs { spcStartTime = 0, 
                 spcStopTime = 13, 
diff --git a/examples/Furnace.hs b/examples/Furnace.hs
new file mode 100644
--- /dev/null
+++ b/examples/Furnace.hs
@@ -0,0 +1,417 @@
+
+import Data.Maybe
+import Random
+import Control.Monad
+import Control.Monad.Trans
+
+import Simulation.Aivika.Dynamics
+import Simulation.Aivika.Dynamics.Base
+import Simulation.Aivika.Dynamics.Lift
+import Simulation.Aivika.Dynamics.EventQueue
+import Simulation.Aivika.Dynamics.Ref
+import Simulation.Aivika.Dynamics.UVar
+import Simulation.Aivika.Dynamics.Process
+import Simulation.Aivika.Dynamics.Random
+import Simulation.Aivika.Statistics
+
+import qualified Simulation.Aivika.Queue as Q
+
+-- | The simulation specs.
+specs = Specs { spcStartTime = 0.0,
+                -- spcStopTime = 1000.0,
+                spcStopTime = 300.0,
+                spcDT = 0.1,
+                spcMethod = RungeKutta4 }
+        
+-- | Return an exponentially distributed random value with mean 
+-- 1 / @lambda@, where @lambda@ is a parameter of the function.
+exprnd :: Double -> IO Double
+exprnd lambda =
+  do x <- getStdRandom random
+     return (- log x / lambda)
+     
+-- | Return a random initial temperature of the item.     
+temprnd :: IO Double
+temprnd =
+  do x <- getStdRandom random
+     return (400.0 + (600.0 - 400.0) * x)
+
+-- | Represents the furnace.
+data Furnace = 
+  Furnace { furnaceQueue :: EventQueue,
+            -- ^ The event queue.
+            furnaceNormalGen :: IO Double,
+            -- ^ The normal random number generator.
+            furnacePits :: [Pit],
+            -- ^ The pits for ingots.
+            furnacePitCount :: UVar Int,
+            -- ^ The count of active pits with ingots.
+            furnacePitCountStats :: Statistics Int,
+            -- ^ The statistics about the active pits.
+            furnaceAwaitingIngots :: Q.Queue Ingot,
+            -- ^ The awaiting ingots in the queue.
+            furnaceQueueCount :: UVar Int,
+            -- ^ The queue count.
+            furnaceQueueCountStats :: Statistics Int,
+            -- ^ The statistics about the queue count.
+            furnaceWaitCount :: Ref Int,
+            -- ^ The count of awaiting ingots.
+            furnaceWaitTime :: Ref Double,
+            -- ^ The wait time for all loaded ingots.
+            furnaceHeatingTime :: Ref Double,
+            -- ^ The heating time for all unloaded ingots.
+            furnaceTemp :: Ref Double,
+            -- ^ The furnace temperature.
+            furnaceTotalCount :: Ref Int,
+            -- ^ The total count of ingots.
+            furnaceLoadCount :: Ref Int,
+            -- ^ The count of loaded ingots.
+            furnaceUnloadCount :: Ref Int,
+            -- ^ The count of unloaded ingots.
+            furnaceUnloadTemps :: Ref [Double]
+            -- ^ The temperatures of all unloaded ingots.
+            }
+
+-- | A pit in the furnace to place the ingots.
+data Pit = 
+  Pit { pitQueue :: EventQueue,
+        -- ^ The bound dynamics queue.
+        pitIngot :: Ref (Maybe Ingot),
+        -- ^ The ingot in the pit.
+        pitTemp :: Ref Double
+        -- ^ The ingot temperature in the pit.
+        }
+
+data Ingot = 
+  Ingot { ingotFurnace :: Furnace,
+          -- ^ Return the furnace.
+          ingotReceiveTime :: Double,
+          -- ^ The time at which the ingot was received.
+          ingotReceiveTemp :: Double,
+          -- ^ The temperature with which the ingot was received.
+          ingotLoadTime :: Double,
+          -- ^ The time of loading in the furnace.
+          ingotLoadTemp :: Double,
+          -- ^ The temperature when the ingot was loaded in the furnace.
+          ingotCoeff :: Double
+          -- ^ The heating coefficient.
+          }
+
+-- | Create a furnace.
+newFurnace :: EventQueue -> Dynamics Furnace
+newFurnace queue =
+  do normalGen <- liftIO normalGen
+     pits <- sequence [newPit queue | i <- [1..10]]
+     pitCount <- newUVar queue 0
+     pitCountStats <- liftIO newStatistics
+     awaitingIngots <- liftIO Q.newQueue
+     queueCount <- newUVar queue 0
+     queueCountStats <- liftIO newStatistics
+     waitCount <- newRef queue 0
+     waitTime <- newRef queue 0.0
+     heatingTime <- newRef queue 0.0
+     h <- newRef queue 1650.0
+     totalCount <- newRef queue 0
+     loadCount <- newRef queue 0
+     unloadCount <- newRef queue 0
+     unloadTemps <- newRef queue []
+     return Furnace { furnaceQueue = queue,
+                      furnaceNormalGen = normalGen,
+                      furnacePits = pits,
+                      furnacePitCount = pitCount,
+                      furnacePitCountStats = pitCountStats,
+                      furnaceAwaitingIngots = awaitingIngots,
+                      furnaceQueueCount = queueCount,
+                      furnaceQueueCountStats = queueCountStats,
+                      furnaceWaitCount = waitCount,
+                      furnaceWaitTime = waitTime,
+                      furnaceHeatingTime = heatingTime,
+                      furnaceTemp = h,
+                      furnaceTotalCount = totalCount,
+                      furnaceLoadCount = loadCount, 
+                      furnaceUnloadCount = unloadCount, 
+                      furnaceUnloadTemps = unloadTemps }
+
+-- | Create a new pit.
+newPit :: EventQueue -> Dynamics Pit
+newPit queue =
+  do ingot <- newRef queue Nothing
+     h' <- newRef queue 0.0
+     return Pit { pitQueue = queue,
+                  pitIngot = ingot,
+                  pitTemp  = h' }
+
+-- | Create a new ingot.
+newIngot :: Furnace -> Dynamics Ingot
+newIngot furnace =
+  do t  <- time
+     xi <- liftIO $ furnaceNormalGen furnace
+     h' <- liftIO temprnd
+     let c = 0.1 + (0.05 + xi * 0.01)
+     return Ingot { ingotFurnace = furnace,
+                    ingotReceiveTime = t,
+                    ingotReceiveTemp = h',
+                    ingotLoadTime = t,
+                    ingotLoadTemp = h',
+                    ingotCoeff = c }
+
+-- | Heat the ingot up in the pit if there is such an ingot.
+heatPitUp :: Pit -> Dynamics ()
+heatPitUp pit =
+  do ingot <- readRef (pitIngot pit)
+     case ingot of
+       Nothing -> 
+         return ()
+       Just ingot -> do
+         
+         -- update the temperature of the ingot.
+         let furnace = ingotFurnace ingot
+         dt' <- dt
+         h'  <- readRef (pitTemp pit)
+         h   <- readRef (furnaceTemp furnace)
+         writeRef (pitTemp pit) $ 
+           h' + dt' * (h - h') * ingotCoeff ingot
+
+-- | Check whether there are ready ingots in the pits.
+ingotsReady :: Furnace -> Dynamics Bool
+ingotsReady furnace =
+  fmap (not . null) $ 
+  filterM (fmap (>= 2200.0) . readRef . pitTemp) $ 
+  furnacePits furnace
+
+-- | Try to unload the ready ingot from the specified pit.
+tryUnloadPit :: Furnace -> Pit -> Dynamics ()
+tryUnloadPit furnace pit =
+  do h' <- readRef (pitTemp pit)
+     when (h' >= 2000.0) $
+       do Just ingot <- readRef (pitIngot pit)  
+          unloadIngot ingot pit
+
+-- | Try to load an awaiting ingot in the specified empty pit.
+tryLoadPit :: Furnace -> Pit -> Dynamics ()       
+tryLoadPit furnace pit =
+  do let ingots = furnaceAwaitingIngots furnace
+     flag <- liftIO $ Q.queueNull ingots
+     unless flag $
+       do ingot <- liftIO $ Q.queueFront ingots
+          liftIO $ Q.dequeue ingots
+          t' <- time
+          modifyUVar (furnaceQueueCount furnace) (+ (-1))
+          c <- readUVar (furnaceQueueCount furnace)
+          liftIO $ addStatistics (furnaceQueueCountStats furnace) c
+          loadIngot (ingot { ingotLoadTime = t',
+                             ingotLoadTemp = 400.0 }) pit
+              
+-- | Unload the ingot from the specified pit.       
+unloadIngot :: Ingot -> Pit -> Dynamics ()
+unloadIngot ingot pit = 
+  do h' <- readRef (pitTemp pit)
+     writeRef (pitIngot pit) Nothing
+     writeRef (pitTemp pit) 0.0
+     
+     -- count the active pits
+     let furnace = ingotFurnace ingot
+     count <- readUVar (furnacePitCount furnace)
+     writeUVar (furnacePitCount furnace) (count - 1)
+     liftIO $ addStatistics (furnacePitCountStats furnace) (count - 1)
+     
+     -- how long did we heat the ingot up?
+     t' <- time
+     modifyRef (furnaceHeatingTime furnace)
+       (+ (t' - ingotLoadTime ingot))
+     
+     -- what is the temperature of the unloaded ingot?
+     modifyRef (furnaceUnloadTemps furnace) (h' :)
+     
+     -- count the unloaded ingots
+     modifyRef (furnaceUnloadCount furnace) (+ 1)
+     
+-- | Load the ingot in the specified pit
+loadIngot :: Ingot -> Pit -> Dynamics ()
+loadIngot ingot pit =
+  do writeRef (pitIngot pit) $ Just ingot
+     writeRef (pitTemp pit) $ ingotLoadTemp ingot
+     
+     -- count the active pits
+     let furnace = ingotFurnace ingot
+     count <- readUVar (furnacePitCount furnace)
+     writeUVar (furnacePitCount furnace) (count + 1)
+     liftIO $ addStatistics (furnacePitCountStats furnace) (count + 1)
+     
+     -- decrease the furnace temperature
+     h <- readRef (furnaceTemp furnace)
+     let h' = ingotLoadTemp ingot
+         dh = - (h - h') / fromInteger (toInteger (count + 1))
+     writeRef (furnaceTemp furnace) $ h + dh
+
+     -- how long did we keep the ingot in the queue?
+     t' <- time
+     when (ingotReceiveTime ingot < t') $
+       do modifyRef (furnaceWaitCount furnace) (+ 1) 
+          modifyRef (furnaceWaitTime furnace)
+            (+ (t' - ingotReceiveTime ingot))
+
+     -- count the loaded ingots
+     modifyRef (furnaceLoadCount furnace) (+ 1)
+  
+-- | Iterate the furnace processing.
+iterateFurnace :: Furnace -> Dynamics (Dynamics ())
+iterateFurnace furnace = 
+  let pits = furnacePits furnace
+  in iterateD $
+     do ready <- ingotsReady furnace
+        when ready $ 
+          do mapM_ (tryUnloadPit furnace) pits
+             pits' <- emptyPits furnace
+             mapM_ (tryLoadPit furnace) pits'
+        mapM_ heatPitUp pits
+        
+        -- update the temperature of the furnace
+        dt' <- dt
+        h   <- readRef (furnaceTemp furnace)
+        writeRef (furnaceTemp furnace) $
+          h + dt' * (2600.0 - h) * 0.2
+
+-- | Return all empty pits.
+emptyPits :: Furnace -> Dynamics [Pit]
+emptyPits furnace =
+  filterM (fmap isNothing . readRef . pitIngot) $
+  furnacePits furnace
+
+-- | Accept a new ingot.
+acceptIngot :: Furnace -> Dynamics ()
+acceptIngot furnace =
+  do ingot <- newIngot furnace
+     
+     -- counting
+     modifyRef (furnaceTotalCount furnace) (+ 1)
+     
+     -- check what to do with the new ingot
+     count <- readUVar (furnacePitCount furnace)
+     if count >= 10
+       then do let ingots = furnaceAwaitingIngots furnace
+               liftIO $ Q.enqueue ingots ingot
+               modifyUVar (furnaceQueueCount furnace) (+ 1)
+               c <- readUVar (furnaceQueueCount furnace)
+               liftIO $ addStatistics (furnaceQueueCountStats furnace) c
+       else do pit:_ <- emptyPits furnace
+               loadIngot ingot pit
+       
+-- | Process the furnace.
+processFurnace :: Furnace -> Process ()
+processFurnace furnace =
+  do delay <- liftIO $ exprnd (1.0 / 2.5)
+     holdProcess delay
+     -- we have got a new ingot
+     liftD $ acceptIngot furnace
+     -- repeat it again
+     processFurnace furnace
+
+-- | Initialize the furnace.
+initializeFurnace :: Furnace -> Dynamics ()
+initializeFurnace furnace =
+  do x1 <- newIngot furnace
+     x2 <- newIngot furnace
+     x3 <- newIngot furnace
+     x4 <- newIngot furnace
+     x5 <- newIngot furnace
+     x6 <- newIngot furnace
+     let p1 : p2 : p3 : p4 : p5 : p6 : ps = 
+           furnacePits furnace
+     loadIngot (x1 { ingotLoadTemp = 550.0 }) p1
+     loadIngot (x2 { ingotLoadTemp = 600.0 }) p2
+     loadIngot (x3 { ingotLoadTemp = 650.0 }) p3
+     loadIngot (x4 { ingotLoadTemp = 700.0 }) p4
+     loadIngot (x5 { ingotLoadTemp = 750.0 }) p5
+     loadIngot (x6 { ingotLoadTemp = 800.0 }) p6
+     writeRef (furnaceTotalCount furnace) 6
+     writeRef (furnaceTemp furnace) 1650.0
+     
+-- | Return a count, average and deviation.
+stats :: [Double] -> (Int, Double, Double)
+stats xs = (length xs, ex, sx)
+  where
+    n  = fromInteger $ toInteger $ length xs
+    ex = sum xs / n
+    dx = (sum . map rho) xs / (n - 1.0)
+    sx = sqrt dx
+    rho x = (x - ex) ^ 2
+
+-- | The simulation model.
+model :: Dynamics (Dynamics ())
+model =
+  do queue <- newQueue
+     furnace <- newFurnace queue
+     pid <- newProcessID queue
+     
+     initializeFurnace furnace
+     
+     -- get the furnace iterator
+     iterator <- iterateFurnace furnace
+     
+     -- accept input ingots
+     t0 <- starttime
+     runProcess (processFurnace furnace) pid t0
+     
+     let system :: Dynamics ()
+         system = 
+           do iterator   --  iterate in each time point
+         
+              -- the ingots
+              c0 <- readRef (furnaceTotalCount furnace)
+              c1 <- readRef (furnaceLoadCount furnace)
+              c2 <- readRef (furnaceUnloadCount furnace)
+              c3 <- readRef (furnaceWaitCount furnace)
+              
+              liftIO $ do
+                putStrLn "The count of ingots:"
+                putStrLn $ "  total  = " ++ show c0
+                putStrLn $ "  loaded = " ++ show c1
+                putStrLn $ "  ready  = " ++ show c2
+                putStrLn $ "  awaited in the queue = " ++ show c3
+                putStrLn ""
+         
+              -- the temperature of the ready ingots
+              (n1, e1, d1) <- 
+                fmap stats $ readRef (furnaceUnloadTemps furnace)
+                
+              liftIO $ do 
+                putStrLn "The temperature of the ready ingots:"
+                putStrLn $ "  average   = " ++ show e1
+                putStrLn $ "  deviation = " ++ show d1
+                putStrLn ""
+                
+              -- the ingots in pits
+              r2 <- fmap analyzeData $ liftIO $ statisticsData (furnacePitCountStats furnace)
+              
+              liftIO $ do
+                putStrLn "The ingots in pits: "
+                putStrLn $ showResults r2 2 []
+                putStrLn ""
+              
+              -- the queue size
+              r3 <- fmap analyzeData $ liftIO $ statisticsData (furnaceQueueCountStats furnace)
+     
+              liftIO $ do
+                putStrLn "The queue size: "
+                putStrLn $ showResults r3 2 []
+                putStrLn ""
+              
+              -- the mean wait time in the queue
+              t4 <- readRef (furnaceWaitTime furnace) /
+                   fmap (fromInteger . toInteger)
+                   (readRef (furnaceWaitCount furnace))
+              
+              -- the mean heating time
+              t5 <- readRef (furnaceHeatingTime furnace) /
+                   fmap (fromInteger . toInteger)
+                   (readRef (furnaceUnloadCount furnace))
+                    
+              liftIO $ do
+                putStrLn $ "The mean wait time: " ++ show t4
+                putStrLn $ "The mean heating time: " ++ show t5
+         
+     return system
+
+-- | The main program.
+main = runDynamics1 model specs
diff --git a/examples/MachRep1.hs b/examples/MachRep1.hs
--- a/examples/MachRep1.hs
+++ b/examples/MachRep1.hs
@@ -19,6 +19,11 @@
 import Control.Monad.Trans
 
 import Simulation.Aivika.Dynamics
+import Simulation.Aivika.Dynamics.Base
+import Simulation.Aivika.Dynamics.Lift
+import Simulation.Aivika.Dynamics.EventQueue
+import Simulation.Aivika.Dynamics.Ref
+import Simulation.Aivika.Dynamics.Process
 
 upRate = 1.0 / 1.0       -- reciprocal of mean up time
 repairRate = 1.0 / 0.5   -- reciprocal of mean repair time
@@ -38,23 +43,25 @@
   do queue <- newQueue
      totalUpTime <- newRef queue 0.0
      
-     pid1 <- newPID queue
-     pid2 <- newPID queue
+     pid1 <- newProcessID queue
+     pid2 <- newProcessID queue
      
-     let machine :: DynamicsProc ()
+     let machine :: Process ()
          machine =
            do startUpTime <- liftD time
               upTime <- liftIO $ exprnd upRate
-              holdProc upTime
+              holdProcess upTime
               finishUpTime <- liftD time
-              liftD $ modifyRef' totalUpTime
+              liftD $ modifyRef totalUpTime
                 (+ (finishUpTime - startUpTime))
               repairTime <- liftIO $ exprnd repairRate
-              holdProc repairTime
+              holdProcess repairTime
               machine
          
-     runProc machine pid1 starttime
-     runProc machine pid2 starttime
+     t0 <- starttime
+     
+     runProcess machine pid1 t0
+     runProcess machine pid2 t0
      
      let system :: Dynamics Double
          system =
diff --git a/examples/MachRep1EventDriven.hs b/examples/MachRep1EventDriven.hs
--- a/examples/MachRep1EventDriven.hs
+++ b/examples/MachRep1EventDriven.hs
@@ -19,6 +19,9 @@
 import Control.Monad.Trans
 
 import Simulation.Aivika.Dynamics
+import Simulation.Aivika.Dynamics.Base
+import Simulation.Aivika.Dynamics.EventQueue
+import Simulation.Aivika.Dynamics.Ref
 
 upRate = 1.0 / 1.0       -- reciprocal of mean up time
 repairRate = 1.0 / 0.5   -- reciprocal of mean repair time
@@ -42,12 +45,12 @@
          machineBroken startUpTime =
            
            do finishUpTime <- time
-              modifyRef' totalUpTime (+ (finishUpTime - startUpTime))
+              modifyRef totalUpTime (+ (finishUpTime - startUpTime))
               repairTime <- liftIO $ exprnd repairRate
               
               -- enqueue a new event
-              let t = return $ finishUpTime + repairTime
-              enqueueD queue t machineRepaired
+              let t = finishUpTime + repairTime
+              enqueue queue t machineRepaired
               
          machineRepaired :: Dynamics ()
          machineRepaired =
@@ -56,11 +59,13 @@
               upTime <- liftIO $ exprnd upRate
               
               -- enqueue a new event
-              let t = return $ startUpTime + upTime
-              enqueueD queue t $ machineBroken startUpTime
+              let t = startUpTime + upTime
+              enqueue queue t $ machineBroken startUpTime
      
-     enqueueD queue starttime machineRepaired    -- start the first machine
-     enqueueD queue starttime machineRepaired    -- start the second machine
+     t0 <- starttime
+     
+     enqueue queue t0 machineRepaired    -- start the first machine
+     enqueue queue t0 machineRepaired    -- start the second machine
      
      let system :: Dynamics Double
          system =
diff --git a/examples/MachRep1TimeDriven.hs b/examples/MachRep1TimeDriven.hs
--- a/examples/MachRep1TimeDriven.hs
+++ b/examples/MachRep1TimeDriven.hs
@@ -19,6 +19,9 @@
 import Control.Monad.Trans
 
 import Simulation.Aivika.Dynamics
+import Simulation.Aivika.Dynamics.Base
+import Simulation.Aivika.Dynamics.EventQueue
+import Simulation.Aivika.Dynamics.Ref
 
 upRate = 1.0 / 1.0       -- reciprocal of mean up time
 repairRate = 1.0 / 0.5   -- reciprocal of mean repair time
@@ -56,34 +59,34 @@
                    repairNum' <- readRef repairNum
                    
                    let untilBroken = 
-                         modifyRef' upNum $ \a -> a - 1
+                         modifyRef upNum $ \a -> a - 1
                                                   
                        untilRepaired =
-                         modifyRef' repairNum $ \a -> a - 1
+                         modifyRef repairNum $ \a -> a - 1
                                                       
                        broken =
-                         do writeRef' upNum (-1)
+                         do writeRef upNum (-1)
                             -- the machine is broken
                             startUpTime' <- readRef startUpTime
                             finishUpTime' <- time
                             dt' <- dt
-                            modifyRef' totalUpTime $ 
+                            modifyRef totalUpTime $ 
                               \a -> a +
                               (finishUpTime' - startUpTime')
                             repairTime' <- 
                               liftIO $ exprnd repairRate
-                            writeRef' repairNum $
+                            writeRef repairNum $
                               round (repairTime' / dt')
                               
                        repaired =
-                         do writeRef' repairNum (-1)
+                         do writeRef repairNum (-1)
                             -- the machine is repaired
                             t'  <- time
                             dt' <- dt
-                            writeRef' startUpTime t'
+                            writeRef startUpTime t'
                             upTime' <- 
                               liftIO $ exprnd upRate
-                            writeRef' upNum $
+                            writeRef upNum $
                               round (upTime' / dt')
                               
                        result | upNum' > 0     = untilBroken
@@ -98,8 +101,8 @@
      m2 <- machine
      
      -- create strictly sequential computations
-     c1 <- memo0 discrete m1
-     c2 <- memo0 discrete m2
+     c1 <- iterateD m1
+     c2 <- iterateD m2
        
      let system :: Dynamics Double
          system =
diff --git a/examples/MachRep2.hs b/examples/MachRep2.hs
--- a/examples/MachRep2.hs
+++ b/examples/MachRep2.hs
@@ -22,6 +22,12 @@
 import Control.Monad.Trans
 
 import Simulation.Aivika.Dynamics
+import Simulation.Aivika.Dynamics.Base
+import Simulation.Aivika.Dynamics.Lift
+import Simulation.Aivika.Dynamics.EventQueue
+import Simulation.Aivika.Dynamics.Ref
+import Simulation.Aivika.Dynamics.Resource
+import Simulation.Aivika.Dynamics.Process
 
 upRate = 1.0 / 1.0       -- reciprocal of mean up time
 repairRate = 1.0 / 0.5   -- reciprocal of mean repair time
@@ -52,33 +58,35 @@
      
      repairPerson <- newResource queue 1
      
-     pid1 <- newPID queue
-     pid2 <- newPID queue
+     pid1 <- newProcessID queue
+     pid2 <- newProcessID queue
      
-     let machine :: DynamicsProc ()
+     let machine :: Process ()
          machine =
            do startUpTime <- liftD time
               upTime <- liftIO $ exprnd upRate
-              holdProc upTime
+              holdProcess upTime
               finishUpTime <- liftD time
-              liftD $ modifyRef' totalUpTime 
+              liftD $ modifyRef totalUpTime 
                 (+ (finishUpTime - startUpTime))
               
               -- check the resource availability
-              liftD $ modifyRef' nRep (+ 1)
+              liftD $ modifyRef nRep (+ 1)
               n <- resourceCount repairPerson
               when (n == 1) $
-                liftD $ modifyRef' nImmedRep (+ 1)
+                liftD $ modifyRef nImmedRep (+ 1)
                 
               requestResource repairPerson
               repairTime <- liftIO $ exprnd repairRate
-              holdProc repairTime
+              holdProcess repairTime
               releaseResource repairPerson
               
               machine
          
-     runProc machine pid1 starttime
-     runProc machine pid2 starttime
+     t0 <- starttime
+     
+     runProcess machine pid1 t0
+     runProcess machine pid2 t0
      
      let system :: Dynamics (Double, Double)
          system =
diff --git a/examples/MachRep3.hs b/examples/MachRep3.hs
--- a/examples/MachRep3.hs
+++ b/examples/MachRep3.hs
@@ -18,6 +18,12 @@
 import Control.Monad.Trans
 
 import Simulation.Aivika.Dynamics
+import Simulation.Aivika.Dynamics.Base
+import Simulation.Aivika.Dynamics.Lift
+import Simulation.Aivika.Dynamics.EventQueue
+import Simulation.Aivika.Dynamics.Ref
+import Simulation.Aivika.Dynamics.Resource
+import Simulation.Aivika.Dynamics.Process
 
 upRate = 1.0 / 1.0       -- reciprocal of mean up time
 repairRate = 1.0 / 0.5   -- reciprocal of mean repair time
@@ -44,35 +50,37 @@
      
      repairPerson <- newResource queue 1
      
-     pid1 <- newPID queue
-     pid2 <- newPID queue
+     pid1 <- newProcessID queue
+     pid2 <- newProcessID queue
      
-     let machine :: DynamicsPID -> DynamicsProc ()
+     let machine :: ProcessID -> Process ()
          machine pid =
            do startUpTime <- liftD time
               upTime <- liftIO $ exprnd upRate
-              holdProc upTime
+              holdProcess upTime
               finishUpTime <- liftD time
-              liftD $ modifyRef' totalUpTime 
+              liftD $ modifyRef totalUpTime 
                 (+ (finishUpTime - startUpTime))
                 
-              liftD $ modifyRef' nUp $ \a -> a - 1
+              liftD $ modifyRef nUp $ \a -> a - 1
               nUp' <- liftD $ readRef nUp
               if nUp' == 1
-                then passivateProc
+                then passivateProcess
                 else do n <- resourceCount repairPerson
-                        when (n == 1) $ reactivateProc pid
+                        when (n == 1) $ reactivateProcess pid
               
               requestResource repairPerson
               repairTime <- liftIO $ exprnd repairRate
-              holdProc repairTime
-              liftD $ modifyRef' nUp $ \a -> a + 1
+              holdProcess repairTime
+              liftD $ modifyRef nUp $ \a -> a + 1
               releaseResource repairPerson
               
               machine pid
 
-     runProc (machine pid2) pid1 starttime
-     runProc (machine pid1) pid2 starttime
+     t0 <- starttime
+     
+     runProcess (machine pid2) pid1 t0
+     runProcess (machine pid1) pid2 t0
      
      let system :: Dynamics Double
          system =