diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+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.
diff --git a/Setup.lhs b/Setup.lhs
new file mode 100644
--- /dev/null
+++ b/Setup.lhs
@@ -0,0 +1,3 @@
+#!/usr/bin/env runhaskell
+> import Distribution.Simple
+> main = defaultMain
diff --git a/Simulation/Aivika/Dynamics.hs b/Simulation/Aivika/Dynamics.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Dynamics.hs
@@ -0,0 +1,1628 @@
+
+-- 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
+  
+          
diff --git a/Simulation/Aivika/PriorityQueue.hs b/Simulation/Aivika/PriorityQueue.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/PriorityQueue.hs
@@ -0,0 +1,173 @@
+
+-- 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 
+       (PriorityQueue, 
+        queueNull, 
+        newQueue, 
+        enqueue, 
+        dequeue, 
+        queueFront) where 
+
+import Data.Array
+import Data.Array.MArray
+import Data.Array.IO
+import Data.IORef
+import Control.Monad
+
+data PriorityQueue a = 
+  PriorityQueue { pqKeys  :: IORef (IOUArray Int Double),
+                  pqVals  :: IORef (IOArray Int a),
+                  pqSize  :: IORef Int,
+                  pqNoVal :: a    -- to release references                 
+                }
+
+increase :: PriorityQueue a -> Int -> IO ()
+increase pq capacity = 
+  do let keyRef = pqKeys pq
+         valRef = pqVals pq
+     keys <- readIORef keyRef
+     vals <- readIORef valRef
+     (il, iu)  <- getBounds keys
+     let len = (iu - il) + 1
+         capacity' | len < 64  = max capacity ((len + 1) * 2)
+                   | otherwise = max capacity ((len `div` 2) * 3)
+         il' = il
+         iu' = il + capacity' - 1
+     keys' <- newArray_ (il', iu')
+     vals' <- newArray_ (il', iu')
+     mapM_ (\i -> do { k <- readArray keys i; writeArray keys' i k }) [il..iu]
+     mapM_ (\i -> do { v <- readArray vals i; writeArray vals' i v }) [il..iu]
+     writeIORef keyRef keys'
+     writeIORef valRef vals'
+
+siftUp :: IOUArray Int Double 
+         -> IOArray Int a
+         -> Int -> Double -> a 
+         -> IO ()
+siftUp keys vals i k v =
+  if i == 0 
+  then do writeArray keys i k
+          writeArray vals i v
+  else do let n = (i - 1) `div` 2
+          kn <- readArray keys n
+          if k >= kn 
+            then do writeArray keys i k
+                    writeArray vals i v
+            else do vn <- readArray vals n
+                    writeArray keys i kn
+                    writeArray vals i vn
+                    siftUp keys vals n k v
+
+siftDown :: IOUArray Int Double 
+           -> IOArray Int a -> Int
+           -> Int -> Double -> a 
+           -> IO ()
+siftDown keys vals size i k v =
+  if i >= (size `div` 2)
+  then do writeArray keys i k
+          writeArray vals i v
+  else do let n  = 2 * i + 1
+              n' = n + 1
+          kn  <- readArray keys n
+          if n' >= size 
+            then if k <= kn
+                 then do writeArray keys i k
+                         writeArray vals i v
+                 else do vn <- readArray vals n
+                         writeArray keys i kn
+                         writeArray vals i vn
+                         siftDown keys vals size n k v
+            else do kn' <- readArray keys n'
+                    let n''  = if kn > kn' then n' else n
+                        kn'' = min kn' kn
+                    if k <= kn''
+                      then do writeArray keys i k
+                              writeArray vals i v
+                      else do vn'' <- readArray vals n''
+                              writeArray keys i kn''
+                              writeArray vals i vn''
+                              siftDown keys vals size n'' k v
+
+queueNull :: PriorityQueue a -> IO Bool
+queueNull pq =
+  do size <- readIORef (pqSize pq)
+     return $ size == 0
+
+newQueue :: a -> IO (PriorityQueue a)
+newQueue defaultValue =
+  do keys <- newArray_ (0, 10)
+     vals <- newArray_ (0, 10)
+     keyRef  <- newIORef keys
+     valRef  <- newIORef vals
+     sizeRef <- newIORef 0
+     return PriorityQueue { pqKeys = keyRef, 
+                            pqVals = valRef, 
+                            pqSize = sizeRef,
+                            pqNoVal = defaultValue }
+
+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)
+     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 :: PriorityQueue a -> IO ()
+dequeue pq =
+  do size <- readIORef (pqSize pq)
+     when (size == 0) $ error "Empty priority queue: dequeue"
+     let i = size - 1
+     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!
+     siftDown keys vals i 0 k v
+
+queueFront :: PriorityQueue a -> IO (Double, a)
+queueFront pq =
+  do size <- readIORef (pqSize pq)
+     when (size == 0) $ error "Empty priority queue: front"
+     keys <- readIORef (pqKeys pq)
+     vals <- readIORef (pqVals pq)
+     k <- readArray keys 0
+     v <- readArray vals 0
+     return (k, v)
diff --git a/Simulation/Aivika/Queue.hs b/Simulation/Aivika/Queue.hs
new file mode 100644
--- /dev/null
+++ b/Simulation/Aivika/Queue.hs
@@ -0,0 +1,112 @@
+
+-- 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 
+       (Queue, 
+        queueNull, 
+        newQueue, 
+        enqueue, 
+        dequeue, 
+        queueFront) where 
+
+import Data.IORef
+import Control.Monad
+
+data QueueItem a = 
+  QueueItem { qiVal  :: a,
+              qiPrev :: IORef (Maybe (QueueItem a)),
+              qiNext :: IORef (Maybe (QueueItem a)) }
+  
+data Queue a =  
+  Queue { qHead :: IORef (Maybe (QueueItem a)),
+          qTail :: IORef (Maybe (QueueItem a)) }
+
+queueNull :: Queue a -> IO Bool
+queueNull q =
+  do head <- readIORef (qHead q) 
+     case head of
+       Nothing -> return True
+       Just _  -> return False
+    
+newQueue :: IO (Queue a)
+newQueue =
+  do head <- newIORef Nothing 
+     tail <- newIORef Nothing
+     return Queue { qHead = head, qTail = tail }
+
+enqueue :: Queue a -> a -> IO ()
+enqueue q v =
+  do head <- readIORef (qHead q)
+     case head of
+       Nothing ->
+         do prev <- newIORef Nothing
+            next <- newIORef Nothing
+            let item = Just QueueItem { qiVal = v, 
+                                        qiPrev = prev, 
+                                        qiNext = next }
+            writeIORef (qHead q) item
+            writeIORef (qTail q) item
+       Just h ->
+         do prev <- newIORef Nothing
+            next <- newIORef head
+            let item = Just QueueItem { qiVal = v,
+                                        qiPrev = prev,
+                                        qiNext = next }
+            writeIORef (qiPrev h) item
+            writeIORef (qHead q) item
+
+dequeue :: Queue a -> IO ()
+dequeue q =
+  do tail <- readIORef (qTail q) 
+     case tail of
+       Nothing ->
+         error "Empty queue: dequeue"
+       Just t ->
+         do tail' <- readIORef (qiPrev t)
+            case tail' of
+              Nothing ->
+                do writeIORef (qHead q) Nothing
+                   writeIORef (qTail q) Nothing
+              Just t' ->
+                do writeIORef (qiNext t') Nothing
+                   writeIORef (qTail q) tail'
+
+queueFront :: Queue a -> IO a
+queueFront q =
+  do tail <- readIORef (qTail q)
+     case tail of
+       Nothing ->
+         error "Empty queue: front"
+       Just t ->
+         return $ qiVal t
diff --git a/aivika.cabal b/aivika.cabal
new file mode 100644
--- /dev/null
+++ b/aivika.cabal
@@ -0,0 +1,47 @@
+name:            aivika
+version:         0.1
+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.
+    .
+    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.
+    .
+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>
+cabal-version:   >= 1.2.0
+build-type:      Simple
+tested-with:     GHC == 6.12.1
+
+extra-source-files:  examples/BassDiffusion.hs
+                     examples/ChemicalReaction.hs
+                     examples/FishBank.hs
+                     examples/MachRep1.hs
+                     examples/MachRep1EventDriven.hs
+                     examples/MachRep1TimeDriven.hs
+                     examples/MachRep2.hs
+                     examples/MachRep3.hs
+
+data-files:          doc/aivika.pdf
+
+library
+    exposed-modules: Simulation.Aivika.Dynamics
+    other-modules:   Simulation.Aivika.Queue
+                     Simulation.Aivika.PriorityQueue
+                     
+    build-depends:   base >= 3 && < 5,
+                     mtl >= 1.1.0.2,
+                     array >= 0.3.0.0
+                     
+    ghc-options:     -O2
diff --git a/doc/aivika.pdf b/doc/aivika.pdf
new file mode 100644
Binary files /dev/null and b/doc/aivika.pdf differ
diff --git a/examples/BassDiffusion.hs b/examples/BassDiffusion.hs
new file mode 100644
--- /dev/null
+++ b/examples/BassDiffusion.hs
@@ -0,0 +1,106 @@
+
+import Random
+import Data.Array
+import Control.Monad
+import Control.Monad.Trans
+
+import Simulation.Aivika.Dynamics
+
+n = 500    -- the number of agents
+
+advertisingEffectiveness = 0.011
+contactRate = 100.0
+adoptionFraction = 0.015
+
+specs = Specs { spcStartTime = 0.0, 
+                spcStopTime = 8.0,
+                spcDT = 0.1,
+                spcMethod = RungeKutta4 }
+
+exprnd :: Double -> IO Double
+exprnd lambda =
+  do x <- getStdRandom random
+     return (- log x / lambda)
+     
+boolrnd :: Double -> IO Bool
+boolrnd p =
+  do x <- getStdRandom random
+     return (x <= p)
+
+data Person = Person { personAgent :: Agent,
+                       personPotentialAdopter :: AgentState,
+                       personAdopter :: AgentState }
+              
+createPerson :: DynamicsQueue -> Dynamics Person              
+createPerson q =    
+  do agent <- newAgent q
+     potentialAdopter <- newState agent
+     adopter <- newState agent
+     return Person { personAgent = agent,
+                     personPotentialAdopter = potentialAdopter,
+                     personAdopter = adopter }
+       
+createPersons :: DynamicsQueue -> 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 p ps potentialAdopters adopters =
+  do stateActivation (personPotentialAdopter p) $
+       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
+          -- add a timer that works while the state is active
+          let t = liftIO $ exprnd contactRate    -- many times!
+          addTimerD (personAdopter p) t $
+            do i <- liftIO $ getStdRandom $ randomR (1, n)
+               let p' = ps ! i
+               st <- agentState (personAgent p')
+               when (st == Just (personPotentialAdopter p')) $
+                 do b <- liftIO $ boolrnd adoptionFraction
+                    when b $ activateState (personAdopter p')
+     stateDeactivation (personPotentialAdopter p) $
+       modifyRef' potentialAdopters $ \a -> a - 1
+     stateDeactivation (personAdopter p) $
+       modifyRef' adopters $ \a -> a - 1
+        
+definePersons :: Array Int Person 
+                -> DynamicsRef Int 
+                -> DynamicsRef Int 
+                -> Dynamics ()
+definePersons ps potentialAdopters adopters =
+  forM_ (elems ps) $ \p -> 
+  definePerson p ps potentialAdopters adopters
+                               
+activatePerson :: Person -> Dynamics ()
+activatePerson p = activateState (personPotentialAdopter p)
+
+activatePersons :: Array Int Person -> Dynamics ()
+activatePersons ps =
+  forM_ (elems ps) $ \p -> activatePerson p
+
+model :: Dynamics (Dynamics [Int])
+model =
+  do q <- newQueue
+     potentialAdopters <- newRef q 0
+     adopters <- newRef q 0
+     ps <- createPersons q
+     definePersons ps potentialAdopters adopters
+     activatePersons ps
+     return $ do i1 <- readRef potentialAdopters
+                 i2 <- readRef adopters
+                 return [i1, i2]
+
+main =
+  do xs <- runDynamics model specs
+     print xs
diff --git a/examples/ChemicalReaction.hs b/examples/ChemicalReaction.hs
new file mode 100644
--- /dev/null
+++ b/examples/ChemicalReaction.hs
@@ -0,0 +1,26 @@
+
+import Simulation.Aivika.Dynamics
+
+specs = Specs { spcStartTime = 0, 
+                spcStopTime = 13, 
+                spcDT = 0.01,
+                spcMethod = RungeKutta4 }
+
+model :: Dynamics (Dynamics [Double])
+model =
+  do integA <- newInteg 100
+     integB <- newInteg 0
+     integC <- newInteg 0
+     let a = integValue integA
+         b = integValue integB
+         c = integValue integC
+     let ka = 1
+         kb = 1
+     integDiff integA (- ka * a)
+     integDiff integB (ka * a - kb * b)
+     integDiff integC (kb * b)
+     return $ sequence [a, b, c]
+
+main = 
+  do a <- runDynamics1 model specs
+     print a
diff --git a/examples/FishBank.hs b/examples/FishBank.hs
new file mode 100644
--- /dev/null
+++ b/examples/FishBank.hs
@@ -0,0 +1,57 @@
+
+import Data.Array
+
+import Simulation.Aivika.Dynamics
+
+specs = Specs { spcStartTime = 0, 
+                spcStopTime = 13, 
+                spcDT = 0.01,
+                -- spcDT = 0.000005,
+                spcMethod = RungeKutta4 }
+
+model :: Dynamics (Dynamics Double)
+model =
+  do fishInteg <- newInteg 1000
+     shipsInteg <- newInteg 10
+     totalProfitInteg <- newInteg 0
+     -- integral values --
+     let fish = integValue fishInteg
+         ships = integValue shipsInteg
+         totalProfit = integValue totalProfitInteg
+     -- auxiliary values --
+     let annualProfit = profit
+         area = 100
+         carryingCapacity = 1000
+         catchPerShip = 
+           lookupD density $
+           listArray (1, 11) [(0.0, -0.048), (1.2, 10.875), (2.4, 17.194), 
+                              (3.6, 20.548), (4.8, 22.086), (6.0, 23.344), 
+                              (7.2, 23.903), (8.4, 24.462), (9.6, 24.882), 
+                              (10.8, 25.301), (12.0, 25.86)]
+         deathFraction = 
+           lookupD (fish / carryingCapacity) $
+           listArray (1, 11) [(0.0, 5.161), (0.1, 5.161), (0.2, 5.161), 
+                              (0.3, 5.161), (0.4, 5.161), (0.5, 5.161), 
+                              (0.6, 5.118), (0.7, 5.247), (0.8, 5.849), 
+                              (0.9, 6.151), (10.0, 6.194)]
+         density = fish / area
+         fishDeathRate = maxD 0 (fish * deathFraction)
+         fishHatchRate = maxD 0 (fish * hatchFraction)
+         fishPrice = 20
+         fractionInvested = 0.2
+         hatchFraction = 6
+         operatingCost = ships * 250
+         profit = revenue - operatingCost
+         revenue = totalCatchPerYear * fishPrice
+         shipBuildingRate = maxD 0 (profit * fractionInvested / shipCost)
+         shipCost = 300
+         totalCatchPerYear = maxD 0 (ships * catchPerShip)
+     -- derivatives --
+     integDiff fishInteg (fishHatchRate - fishDeathRate - totalCatchPerYear)
+     integDiff shipsInteg shipBuildingRate
+     integDiff totalProfitInteg annualProfit
+     -- results --
+     return annualProfit
+
+main = do a <- runDynamics1 model specs
+          print a
diff --git a/examples/MachRep1.hs b/examples/MachRep1.hs
new file mode 100644
--- /dev/null
+++ b/examples/MachRep1.hs
@@ -0,0 +1,69 @@
+
+-- It corresponds to model MachRep1 described in document 
+-- Introduction to Discrete-Event Simulation and the SimPy Language
+-- [http://heather.cs.ucdavis.edu/~matloff/156/PLN/DESimIntro.pdf]. 
+-- SimPy is available on [http://simpy.sourceforge.net/].
+--   
+-- The model description is as follows.
+--
+-- Two machines, which sometimes break down.
+-- Up time is exponentially distributed with mean 1.0, and repair time is
+-- exponentially distributed with mean 0.5. There are two repairpersons,
+-- so the two machines can be repaired simultaneously if they are down
+-- at the same time.
+--
+-- Output is long-run proportion of up time. Should get value of about
+-- 0.66.
+
+import Random
+import Control.Monad.Trans
+
+import Simulation.Aivika.Dynamics
+
+upRate = 1.0 / 1.0       -- reciprocal of mean up time
+repairRate = 1.0 / 0.5   -- reciprocal of mean repair time
+
+specs = Specs { spcStartTime = 0.0,
+                spcStopTime = 1000.0,
+                spcDT = 1.0,
+                spcMethod = RungeKutta4 }
+        
+exprnd :: Double -> IO Double
+exprnd lambda =
+  do x <- getStdRandom random
+     return (- log x / lambda)
+     
+model :: Dynamics (Dynamics Double)
+model =
+  do queue <- newQueue
+     totalUpTime <- newRef queue 0.0
+     
+     pid1 <- newPID queue
+     pid2 <- newPID queue
+     
+     let machine :: DynamicsProc ()
+         machine =
+           do startUpTime <- liftD time
+              upTime <- liftIO $ exprnd upRate
+              holdProc upTime
+              finishUpTime <- liftD time
+              liftD $ modifyRef' totalUpTime
+                (+ (finishUpTime - startUpTime))
+              repairTime <- liftIO $ exprnd repairRate
+              holdProc repairTime
+              machine
+         
+     runProc machine pid1 starttime
+     runProc machine pid2 starttime
+     
+     let system :: Dynamics Double
+         system =
+           do x <- readRef totalUpTime
+              y <- stoptime
+              return $ x / (2 * y)
+     
+     return system
+  
+main =         
+  do a <- runDynamics1 model specs
+     print a
diff --git a/examples/MachRep1EventDriven.hs b/examples/MachRep1EventDriven.hs
new file mode 100644
--- /dev/null
+++ b/examples/MachRep1EventDriven.hs
@@ -0,0 +1,75 @@
+
+-- It corresponds to model MachRep1 described in document 
+-- Introduction to Discrete-Event Simulation and the SimPy Language
+-- [http://heather.cs.ucdavis.edu/~matloff/156/PLN/DESimIntro.pdf]. 
+-- SimPy is available on [http://simpy.sourceforge.net/].
+--   
+-- The model description is as follows.
+--
+-- Two machines, which sometimes break down.
+-- Up time is exponentially distributed with mean 1.0, and repair time is
+-- exponentially distributed with mean 0.5. There are two repairpersons,
+-- so the two machines can be repaired simultaneously if they are down
+-- at the same time.
+--
+-- Output is long-run proportion of up time. Should get value of about
+-- 0.66.
+
+import Random
+import Control.Monad.Trans
+
+import Simulation.Aivika.Dynamics
+
+upRate = 1.0 / 1.0       -- reciprocal of mean up time
+repairRate = 1.0 / 0.5   -- reciprocal of mean repair time
+
+specs = Specs { spcStartTime = 0.0,
+                spcStopTime = 1000.0,
+                spcDT = 1.0,
+                spcMethod = RungeKutta4 }
+        
+exprnd :: Double -> IO Double
+exprnd lambda =
+  do x <- getStdRandom random
+     return (- log x / lambda)
+     
+model :: Dynamics (Dynamics Double)
+model =
+  do queue <- newQueue
+     totalUpTime <- newRef queue 0.0
+     
+     let machineBroken :: Double -> Dynamics ()
+         machineBroken startUpTime =
+           
+           do finishUpTime <- time
+              modifyRef' totalUpTime (+ (finishUpTime - startUpTime))
+              repairTime <- liftIO $ exprnd repairRate
+              
+              -- enqueue a new event
+              let t = return $ finishUpTime + repairTime
+              enqueueD queue t machineRepaired
+              
+         machineRepaired :: Dynamics ()
+         machineRepaired =
+           
+           do startUpTime <- time
+              upTime <- liftIO $ exprnd upRate
+              
+              -- enqueue a new event
+              let t = return $ startUpTime + upTime
+              enqueueD queue t $ machineBroken startUpTime
+     
+     enqueueD queue starttime machineRepaired    -- start the first machine
+     enqueueD queue starttime machineRepaired    -- start the second machine
+     
+     let system :: Dynamics Double
+         system =
+           do x <- readRef totalUpTime
+              y <- stoptime
+              return $ x / (2 * y)
+              
+     return system
+  
+main =         
+  do a <- runDynamics1 model specs
+     print a
diff --git a/examples/MachRep1TimeDriven.hs b/examples/MachRep1TimeDriven.hs
new file mode 100644
--- /dev/null
+++ b/examples/MachRep1TimeDriven.hs
@@ -0,0 +1,116 @@
+
+-- It corresponds to model MachRep1 described in document 
+-- Introduction to Discrete-Event Simulation and the SimPy Language
+-- [http://heather.cs.ucdavis.edu/~matloff/156/PLN/DESimIntro.pdf]. 
+-- SimPy is available on [http://simpy.sourceforge.net/].
+--   
+-- The model description is as follows.
+--
+-- Two machines, which sometimes break down.
+-- Up time is exponentially distributed with mean 1.0, and repair time is
+-- exponentially distributed with mean 0.5. There are two repairpersons,
+-- so the two machines can be repaired simultaneously if they are down
+-- at the same time.
+--
+-- Output is long-run proportion of up time. Should get value of about
+-- 0.66.
+
+import Random
+import Control.Monad.Trans
+
+import Simulation.Aivika.Dynamics
+
+upRate = 1.0 / 1.0       -- reciprocal of mean up time
+repairRate = 1.0 / 0.5   -- reciprocal of mean repair time
+
+specs = Specs { spcStartTime = 0.0,
+                spcStopTime = 1000.0,
+                spcDT = 0.05,
+                spcMethod = RungeKutta4 }
+        
+exprnd :: Double -> IO Double
+exprnd lambda =
+  do x <- getStdRandom random
+     return (- log x / lambda)
+     
+model :: Dynamics (Dynamics Double)
+model =
+  do queue <- newQueue
+     totalUpTime <- newRef queue 0.0
+     
+     let machine :: Dynamics (Dynamics ())
+         machine =
+           do startUpTime <- newRef queue 0.0 
+             
+              -- a number of iterations when 
+              -- the machine works
+              upNum <- newRef queue (-1)
+              
+              -- a number of iterations when 
+              -- the machine is broken
+              repairNum <- newRef queue (-1)
+              
+              -- create a simulation model
+              return $
+                do upNum' <- readRef upNum
+                   repairNum' <- readRef repairNum
+                   
+                   let untilBroken = 
+                         modifyRef' upNum $ \a -> a - 1
+                                                  
+                       untilRepaired =
+                         modifyRef' repairNum $ \a -> a - 1
+                                                      
+                       broken =
+                         do writeRef' upNum (-1)
+                            -- the machine is broken
+                            startUpTime' <- readRef startUpTime
+                            finishUpTime' <- time
+                            dt' <- dt
+                            modifyRef' totalUpTime $ 
+                              \a -> a +
+                              (finishUpTime' - startUpTime')
+                            repairTime' <- 
+                              liftIO $ exprnd repairRate
+                            writeRef' repairNum $
+                              round (repairTime' / dt')
+                              
+                       repaired =
+                         do writeRef' repairNum (-1)
+                            -- the machine is repaired
+                            t'  <- time
+                            dt' <- dt
+                            writeRef' startUpTime t'
+                            upTime' <- 
+                              liftIO $ exprnd upRate
+                            writeRef' upNum $
+                              round (upTime' / dt')
+                              
+                       result | upNum' > 0     = untilBroken
+                              | upNum' == 0     = broken
+                              | repairNum' > 0 = untilRepaired
+                              | repairNum' == 0 = repaired
+                              | otherwise      = repaired 
+                   result
+                            
+     -- create two machines with type Dynamics ()
+     m1 <- machine
+     m2 <- machine
+     
+     -- create strictly sequential computations
+     c1 <- memo0 discrete m1
+     c2 <- memo0 discrete m2
+       
+     let system :: Dynamics Double
+         system =
+           do c1    -- involve in the simulation
+              c2    -- involve in the simulation
+              x <- readRef totalUpTime
+              y <- stoptime
+              return $ x / (2 * y)
+     
+     return system
+  
+main =         
+  do a <- runDynamics1 model specs
+     print a
diff --git a/examples/MachRep2.hs b/examples/MachRep2.hs
new file mode 100644
--- /dev/null
+++ b/examples/MachRep2.hs
@@ -0,0 +1,96 @@
+
+-- It corresponds to model MachRep2 described in document 
+-- Introduction to Discrete-Event Simulation and the SimPy Language
+-- [http://heather.cs.ucdavis.edu/~matloff/156/PLN/DESimIntro.pdf]. 
+-- SimPy is available on [http://simpy.sourceforge.net/].
+--   
+-- The model description is as follows.
+--   
+-- Two machines, but sometimes break down. Up time is exponentially 
+-- distributed with mean 1.0, and repair time is exponentially distributed 
+-- with mean 0.5. In this example, there is only one repairperson, so 
+-- the two machines cannot be repaired simultaneously if they are down 
+-- at the same time.
+--
+-- In addition to finding the long-run proportion of up time as in
+-- model MachRep1, let’s also find the long-run proportion of the time 
+-- that a given machine does not have immediate access to the repairperson 
+-- when the machine breaks down. Output values should be about 0.6 and 0.67. 
+
+import Random
+import Control.Monad
+import Control.Monad.Trans
+
+import Simulation.Aivika.Dynamics
+
+upRate = 1.0 / 1.0       -- reciprocal of mean up time
+repairRate = 1.0 / 0.5   -- reciprocal of mean repair time
+
+specs = Specs { spcStartTime = 0.0,
+                spcStopTime = 1000.0,
+                spcDT = 1.0,
+                spcMethod = RungeKutta4 }
+        
+exprnd :: Double -> IO Double
+exprnd lambda =
+  do x <- getStdRandom random
+     return (- log x / lambda)
+     
+model :: Dynamics (Dynamics (Double, Double))
+model =
+  do queue <- newQueue
+     
+     -- number of times the machines have broken down
+     nRep <- newRef queue 0 
+     
+     -- number of breakdowns in which the machine 
+     -- started repair service right away
+     nImmedRep <- newRef queue 0
+     
+     -- total up time for all machines
+     totalUpTime <- newRef queue 0.0
+     
+     repairPerson <- newResource queue 1
+     
+     pid1 <- newPID queue
+     pid2 <- newPID queue
+     
+     let machine :: DynamicsProc ()
+         machine =
+           do startUpTime <- liftD time
+              upTime <- liftIO $ exprnd upRate
+              holdProc upTime
+              finishUpTime <- liftD time
+              liftD $ modifyRef' totalUpTime 
+                (+ (finishUpTime - startUpTime))
+              
+              -- check the resource availability
+              liftD $ modifyRef' nRep (+ 1)
+              n <- resourceCount repairPerson
+              when (n == 1) $
+                liftD $ modifyRef' nImmedRep (+ 1)
+                
+              requestResource repairPerson
+              repairTime <- liftIO $ exprnd repairRate
+              holdProc repairTime
+              releaseResource repairPerson
+              
+              machine
+         
+     runProc machine pid1 starttime
+     runProc machine pid2 starttime
+     
+     let system :: Dynamics (Double, Double)
+         system =
+           do x <- readRef totalUpTime
+              y <- stoptime
+              n <- readRef nRep
+              nImmed <- readRef nImmedRep
+              return (x / (2 * y), 
+                      fromIntegral nImmed / fromIntegral n)
+     
+     return system
+  
+main =         
+  do a <- runDynamics1 model specs
+     print a
diff --git a/examples/MachRep3.hs b/examples/MachRep3.hs
new file mode 100644
--- /dev/null
+++ b/examples/MachRep3.hs
@@ -0,0 +1,87 @@
+
+-- It corresponds to model MachRep3 described in document 
+-- Introduction to Discrete-Event Simulation and the SimPy Language
+-- [http://heather.cs.ucdavis.edu/~matloff/156/PLN/DESimIntro.pdf]. 
+-- SimPy is available on [http://simpy.sourceforge.net/].
+--   
+-- The model description is as follows.
+--
+-- Variation of models MachRep1, MachRep2. Two machines, but
+-- sometimes break down. Up time is exponentially distributed with mean
+-- 1.0, and repair time is exponentially distributed with mean 0.5. In
+-- this example, there is only one repairperson, and she is not summoned
+-- until both machines are down. We find the proportion of up time. It
+-- should come out to about 0.45.
+
+import Random
+import Control.Monad
+import Control.Monad.Trans
+
+import Simulation.Aivika.Dynamics
+
+upRate = 1.0 / 1.0       -- reciprocal of mean up time
+repairRate = 1.0 / 0.5   -- reciprocal of mean repair time
+
+specs = Specs { spcStartTime = 0.0,
+                spcStopTime = 1000.0,
+                spcDT = 1.0,
+                spcMethod = RungeKutta4 }
+        
+exprnd :: Double -> IO Double
+exprnd lambda =
+  do x <- getStdRandom random
+     return (- log x / lambda)
+     
+model :: Dynamics (Dynamics Double)
+model =
+  do queue <- newQueue
+     
+     -- number of machines currently up
+     nUp <- newRef queue 2
+     
+     -- total up time for all machines
+     totalUpTime <- newRef queue 0.0
+     
+     repairPerson <- newResource queue 1
+     
+     pid1 <- newPID queue
+     pid2 <- newPID queue
+     
+     let machine :: DynamicsPID -> DynamicsProc ()
+         machine pid =
+           do startUpTime <- liftD time
+              upTime <- liftIO $ exprnd upRate
+              holdProc upTime
+              finishUpTime <- liftD time
+              liftD $ modifyRef' totalUpTime 
+                (+ (finishUpTime - startUpTime))
+                
+              liftD $ modifyRef' nUp $ \a -> a - 1
+              nUp' <- liftD $ readRef nUp
+              if nUp' == 1
+                then passivateProc
+                else do n <- resourceCount repairPerson
+                        when (n == 1) $ reactivateProc pid
+              
+              requestResource repairPerson
+              repairTime <- liftIO $ exprnd repairRate
+              holdProc repairTime
+              liftD $ modifyRef' nUp $ \a -> a + 1
+              releaseResource repairPerson
+              
+              machine pid
+
+     runProc (machine pid2) pid1 starttime
+     runProc (machine pid1) pid2 starttime
+     
+     let system :: Dynamics Double
+         system =
+           do x <- readRef totalUpTime
+              y <- stoptime
+              return $ x / (2 * y)
+     
+     return system
+  
+main =         
+  do a <- runDynamics1 model specs
+     print a