module Integrator (
  newIntegrator, input, output, stop,
  Time, timeInterval
) where
import Control.Concurrent (forkIO, threadDelay)
import Control.Concurrent.MVar (MVar, newMVar, modifyMVar_, modifyMVar, readMVar)
import Control.Exception (evaluate)
import Data.Time (UTCTime)
import Data.Time.Clock (getCurrentTime, diffUTCTime)


main = do let f = 0.5 
          t0 <- getCurrentTime
          i <- newIntegrator
          input i (\t -> sin(2*pi * f * timeInterval t0 t)) 
          threadDelay 2000000                       
          input i (const 0)                                 
          threadDelay 500000                        
          result <- output i
          stop i
          print result




type Time = UTCTime
type Func a = Time -> a
timeInterval t0 t1 = realToFrac $ diffUTCTime t1 t0


newIntegrator :: Fractional a => IO (Integrator a) 
input  :: Integrator a -> Func a -> IO ()          
output :: Integrator a           -> IO a           
stop   :: Integrator a           -> IO ()          


data Integrator a = Integrator (MVar (IntState a)) 
  deriving Eq
data IntState a = IntState { func  :: Func a,      
                             run   :: Bool,        
                             value :: a,           
                             time  :: Time }       

newIntegrator = do
  now <- getCurrentTime
  state <- newMVar $ IntState { func  = const 0,
                                run   = True,
                                value = 0,
                                time  = now }
  thread <- forkIO (intThread state)  
  return (Integrator state)           

input  (Integrator stv) f = modifyMVar_ stv (\st -> return st { func = f })
output (Integrator stv)   = fmap value $ readMVar stv
stop   (Integrator stv)   = modifyMVar_ stv (\st -> return st { run = False })
  
  


intThread :: Fractional a => MVar (IntState a) -> IO ()
intThread stv = whileM $ modifyMVar stv updateAndCheckRun
  
  
  
  where updateAndCheckRun st = do
          now <- getCurrentTime
          let value' = integrate (func st) (value st) (time st) now
          evaluate value'                             
          return (st { value = value', time  = now }, 
                  run st)                             

integrate:: Fractional a => Func a -> a -> Time -> Time -> a
integrate f value t0 t1 = value + (f t0 + f t1)/2 * dt
  where dt = timeInterval t0 t1


whileM action = do b <- action; if b then whileM action else return ()