{-# LANGUAGE BangPatterns, RecursiveDo #-}

-- |
-- Module     : Simulation.Aivika.SystemDynamics
-- Copyright  : Copyright (c) 2009-2017, David Sorokin <david.sorokin@gmail.com>
-- License    : BSD3
-- Maintainer : David Sorokin <david.sorokin@gmail.com>
-- Stability  : experimental
-- Tested with: GHC 8.0.1
--
-- This module defines integrals and other functions of System Dynamics.
--

module Simulation.Aivika.SystemDynamics
       (-- * Equality and Ordering
        (.==.),
        (./=.),
        (.<.),
        (.>=.),
        (.>.),
        (.<=.),
        maxDynamics,
        minDynamics,
        ifDynamics,
        -- * Ordinary Differential Equations
        integ,
        integEither,
        smoothI,
        smooth,
        smooth3I,
        smooth3,
        smoothNI,
        smoothN,
        delay1I,
        delay1,
        delay3I,
        delay3,
        delayNI,
        delayN,
        forecast,
        trend,
        -- * Difference Equations
        diffsum,
        diffsumEither,
        -- * Table Functions
        lookupDynamics,
        lookupStepwiseDynamics,
        -- * Discrete Functions
        delay,
        delayI,
        delayByDT,
        delayIByDT,
        step,
        pulse,
        pulseP,
        ramp,
        -- * Financial Functions
        npv,
        npve) where

import Data.Array
import Data.Array.IO.Safe
import Data.IORef
import Control.Monad
import Control.Monad.Trans

import Simulation.Aivika.Internal.Specs
import Simulation.Aivika.Internal.Parameter
import Simulation.Aivika.Internal.Simulation
import Simulation.Aivika.Internal.Dynamics
import Simulation.Aivika.Dynamics.Extra
import Simulation.Aivika.Unboxed
import Simulation.Aivika.Table

import qualified Simulation.Aivika.Dynamics.Memo as M
import qualified Simulation.Aivika.Dynamics.Memo.Unboxed as MU

--
-- Equality and Ordering
--

-- | Compare for equality.
(.==.) :: (Eq a) => Dynamics a -> Dynamics a -> Dynamics Bool
.==. :: forall a. Eq a => Dynamics a -> Dynamics a -> Dynamics Bool
(.==.) = forall (m :: * -> *) a1 a2 r.
Monad m =>
(a1 -> a2 -> r) -> m a1 -> m a2 -> m r
liftM2 forall a. Eq a => a -> a -> Bool
(==)

-- | Compare for inequality.
(./=.) :: (Eq a) => Dynamics a -> Dynamics a -> Dynamics Bool
./=. :: forall a. Eq a => Dynamics a -> Dynamics a -> Dynamics Bool
(./=.) = forall (m :: * -> *) a1 a2 r.
Monad m =>
(a1 -> a2 -> r) -> m a1 -> m a2 -> m r
liftM2 forall a. Eq a => a -> a -> Bool
(/=)

-- | Compare for ordering.
(.<.) :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics Bool
.<. :: forall a. Ord a => Dynamics a -> Dynamics a -> Dynamics Bool
(.<.) = forall (m :: * -> *) a1 a2 r.
Monad m =>
(a1 -> a2 -> r) -> m a1 -> m a2 -> m r
liftM2 forall a. Ord a => a -> a -> Bool
(<)

-- | Compare for ordering.
(.>=.) :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics Bool
.>=. :: forall a. Ord a => Dynamics a -> Dynamics a -> Dynamics Bool
(.>=.) = forall (m :: * -> *) a1 a2 r.
Monad m =>
(a1 -> a2 -> r) -> m a1 -> m a2 -> m r
liftM2 forall a. Ord a => a -> a -> Bool
(>=)

-- | Compare for ordering.
(.>.) :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics Bool
.>. :: forall a. Ord a => Dynamics a -> Dynamics a -> Dynamics Bool
(.>.) = forall (m :: * -> *) a1 a2 r.
Monad m =>
(a1 -> a2 -> r) -> m a1 -> m a2 -> m r
liftM2 forall a. Ord a => a -> a -> Bool
(>)

-- | Compare for ordering.
(.<=.) :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics Bool
.<=. :: forall a. Ord a => Dynamics a -> Dynamics a -> Dynamics Bool
(.<=.) = forall (m :: * -> *) a1 a2 r.
Monad m =>
(a1 -> a2 -> r) -> m a1 -> m a2 -> m r
liftM2 forall a. Ord a => a -> a -> Bool
(<=)

-- | Return the maximum.
maxDynamics :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics a
maxDynamics :: forall a. Ord a => Dynamics a -> Dynamics a -> Dynamics a
maxDynamics = forall (m :: * -> *) a1 a2 r.
Monad m =>
(a1 -> a2 -> r) -> m a1 -> m a2 -> m r
liftM2 forall a. Ord a => a -> a -> a
max

-- | Return the minimum.
minDynamics :: (Ord a) => Dynamics a -> Dynamics a -> Dynamics a
minDynamics :: forall a. Ord a => Dynamics a -> Dynamics a -> Dynamics a
minDynamics = forall (m :: * -> *) a1 a2 r.
Monad m =>
(a1 -> a2 -> r) -> m a1 -> m a2 -> m r
liftM2 forall a. Ord a => a -> a -> a
min

-- | Implement the if-then-else operator.
ifDynamics :: Dynamics Bool -> Dynamics a -> Dynamics a -> Dynamics a
ifDynamics :: forall a. Dynamics Bool -> Dynamics a -> Dynamics a -> Dynamics a
ifDynamics Dynamics Bool
cond Dynamics a
x Dynamics a
y =
  do Bool
a <- Dynamics Bool
cond
     if Bool
a then Dynamics a
x else Dynamics a
y

--
-- Ordinary Differential Equations
--

integEuler :: Dynamics Double
             -> Dynamics Double 
             -> Dynamics Double 
             -> Point -> IO Double
integEuler :: Dynamics Double
-> Dynamics Double -> Dynamics Double -> Point -> IO Double
integEuler (Dynamics Point -> IO Double
f) (Dynamics Point -> IO Double
i) (Dynamics Point -> IO Double
y) Point
p = 
  case Point -> Int
pointIteration Point
p of
    Int
0 -> 
      Point -> IO Double
i Point
p
    Int
n -> do 
      let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
          ty :: Double
ty = Specs -> Int -> Int -> Double
basicTime Specs
sc (Int
n forall a. Num a => a -> a -> a
- Int
1) Int
0
          py :: Point
py = Point
p { pointTime :: Double
pointTime = Double
ty, pointIteration :: Int
pointIteration = Int
n forall a. Num a => a -> a -> a
- Int
1, pointPhase :: Int
pointPhase = Int
0 }
      Double
a <- Point -> IO Double
y Point
py
      Double
b <- Point -> IO Double
f Point
py
      let !v :: Double
v = Double
a forall a. Num a => a -> a -> a
+ Specs -> Double
spcDT (Point -> Specs
pointSpecs Point
p) forall a. Num a => a -> a -> a
* Double
b
      forall (m :: * -> *) a. Monad m => a -> m a
return Double
v

integRK2 :: Dynamics Double
           -> Dynamics Double
           -> Dynamics Double
           -> Point -> IO Double
integRK2 :: Dynamics Double
-> Dynamics Double -> Dynamics Double -> Point -> IO Double
integRK2 (Dynamics Point -> IO Double
f) (Dynamics Point -> IO Double
i) (Dynamics Point -> IO Double
y) Point
p =
  case Point -> Int
pointPhase Point
p of
    Int
0 -> case Point -> Int
pointIteration Point
p of
      Int
0 ->
        Point -> IO Double
i Point
p
      Int
n -> do
        let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
            ty :: Double
ty = Specs -> Int -> Int -> Double
basicTime Specs
sc (Int
n forall a. Num a => a -> a -> a
- Int
1) Int
0
            t1 :: Double
t1 = Double
ty
            t2 :: Double
t2 = Specs -> Int -> Int -> Double
basicTime Specs
sc (Int
n forall a. Num a => a -> a -> a
- Int
1) Int
1
            py :: Point
py = Point
p { pointTime :: Double
pointTime = Double
ty, pointIteration :: Int
pointIteration = Int
n forall a. Num a => a -> a -> a
- Int
1, pointPhase :: Int
pointPhase = Int
0 }
            p1 :: Point
p1 = Point
py
            p2 :: Point
p2 = Point
p { pointTime :: Double
pointTime = Double
t2, pointIteration :: Int
pointIteration = Int
n forall a. Num a => a -> a -> a
- Int
1, pointPhase :: Int
pointPhase = Int
1 }
        Double
vy <- Point -> IO Double
y Point
py
        Double
k1 <- Point -> IO Double
f Point
p1
        Double
k2 <- Point -> IO Double
f Point
p2
        let !v :: Double
v = Double
vy forall a. Num a => a -> a -> a
+ Specs -> Double
spcDT Specs
sc forall a. Fractional a => a -> a -> a
/ Double
2.0 forall a. Num a => a -> a -> a
* (Double
k1 forall a. Num a => a -> a -> a
+ Double
k2)
        forall (m :: * -> *) a. Monad m => a -> m a
return Double
v
    Int
1 -> do
      let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
          n :: Int
n  = Point -> Int
pointIteration Point
p
          ty :: Double
ty = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
0
          t1 :: Double
t1 = Double
ty
          py :: Point
py = Point
p { pointTime :: Double
pointTime = Double
ty, pointIteration :: Int
pointIteration = Int
n, pointPhase :: Int
pointPhase = Int
0 }
          p1 :: Point
p1 = Point
py
      Double
vy <- Point -> IO Double
y Point
py
      Double
k1 <- Point -> IO Double
f Point
p1
      let !v :: Double
v = Double
vy forall a. Num a => a -> a -> a
+ Specs -> Double
spcDT Specs
sc forall a. Num a => a -> a -> a
* Double
k1
      forall (m :: * -> *) a. Monad m => a -> m a
return Double
v
    Int
_ -> 
      forall a. HasCallStack => [Char] -> a
error [Char]
"Incorrect phase: integRK2"

integRK4 :: Dynamics Double
           -> Dynamics Double
           -> Dynamics Double
           -> Point -> IO Double
integRK4 :: Dynamics Double
-> Dynamics Double -> Dynamics Double -> Point -> IO Double
integRK4 (Dynamics Point -> IO Double
f) (Dynamics Point -> IO Double
i) (Dynamics Point -> IO Double
y) Point
p =
  case Point -> Int
pointPhase Point
p of
    Int
0 -> case Point -> Int
pointIteration Point
p of
      Int
0 -> 
        Point -> IO Double
i Point
p
      Int
n -> do
        let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
            ty :: Double
ty = Specs -> Int -> Int -> Double
basicTime Specs
sc (Int
n forall a. Num a => a -> a -> a
- Int
1) Int
0
            t1 :: Double
t1 = Double
ty
            t2 :: Double
t2 = Specs -> Int -> Int -> Double
basicTime Specs
sc (Int
n forall a. Num a => a -> a -> a
- Int
1) Int
1
            t3 :: Double
t3 = Specs -> Int -> Int -> Double
basicTime Specs
sc (Int
n forall a. Num a => a -> a -> a
- Int
1) Int
2
            t4 :: Double
t4 = Specs -> Int -> Int -> Double
basicTime Specs
sc (Int
n forall a. Num a => a -> a -> a
- Int
1) Int
3
            py :: Point
py = Point
p { pointTime :: Double
pointTime = Double
ty, pointIteration :: Int
pointIteration = Int
n forall a. Num a => a -> a -> a
- Int
1, pointPhase :: Int
pointPhase = Int
0 }
            p1 :: Point
p1 = Point
py
            p2 :: Point
p2 = Point
p { pointTime :: Double
pointTime = Double
t2, pointIteration :: Int
pointIteration = Int
n forall a. Num a => a -> a -> a
- Int
1, pointPhase :: Int
pointPhase = Int
1 }
            p3 :: Point
p3 = Point
p { pointTime :: Double
pointTime = Double
t3, pointIteration :: Int
pointIteration = Int
n forall a. Num a => a -> a -> a
- Int
1, pointPhase :: Int
pointPhase = Int
2 }
            p4 :: Point
p4 = Point
p { pointTime :: Double
pointTime = Double
t4, pointIteration :: Int
pointIteration = Int
n forall a. Num a => a -> a -> a
- Int
1, pointPhase :: Int
pointPhase = Int
3 }
        Double
vy <- Point -> IO Double
y Point
py
        Double
k1 <- Point -> IO Double
f Point
p1
        Double
k2 <- Point -> IO Double
f Point
p2
        Double
k3 <- Point -> IO Double
f Point
p3
        Double
k4 <- Point -> IO Double
f Point
p4
        let !v :: Double
v = Double
vy forall a. Num a => a -> a -> a
+ Specs -> Double
spcDT Specs
sc forall a. Fractional a => a -> a -> a
/ Double
6.0 forall a. Num a => a -> a -> a
* (Double
k1 forall a. Num a => a -> a -> a
+ Double
2.0 forall a. Num a => a -> a -> a
* Double
k2 forall a. Num a => a -> a -> a
+ Double
2.0 forall a. Num a => a -> a -> a
* Double
k3 forall a. Num a => a -> a -> a
+ Double
k4)
        forall (m :: * -> *) a. Monad m => a -> m a
return Double
v
    Int
1 -> do
      let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
          n :: Int
n  = Point -> Int
pointIteration Point
p
          ty :: Double
ty = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
0
          t1 :: Double
t1 = Double
ty
          py :: Point
py = Point
p { pointTime :: Double
pointTime = Double
ty, pointIteration :: Int
pointIteration = Int
n, pointPhase :: Int
pointPhase = Int
0 }
          p1 :: Point
p1 = Point
py
      Double
vy <- Point -> IO Double
y Point
py
      Double
k1 <- Point -> IO Double
f Point
p1
      let !v :: Double
v = Double
vy forall a. Num a => a -> a -> a
+ Specs -> Double
spcDT Specs
sc forall a. Fractional a => a -> a -> a
/ Double
2.0 forall a. Num a => a -> a -> a
* Double
k1
      forall (m :: * -> *) a. Monad m => a -> m a
return Double
v
    Int
2 -> do
      let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
          n :: Int
n  = Point -> Int
pointIteration Point
p
          ty :: Double
ty = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
0
          t2 :: Double
t2 = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
1
          py :: Point
py = Point
p { pointTime :: Double
pointTime = Double
ty, pointIteration :: Int
pointIteration = Int
n, pointPhase :: Int
pointPhase = Int
0 }
          p2 :: Point
p2 = Point
p { pointTime :: Double
pointTime = Double
t2, pointIteration :: Int
pointIteration = Int
n, pointPhase :: Int
pointPhase = Int
1 }
      Double
vy <- Point -> IO Double
y Point
py
      Double
k2 <- Point -> IO Double
f Point
p2
      let !v :: Double
v = Double
vy forall a. Num a => a -> a -> a
+ Specs -> Double
spcDT Specs
sc forall a. Fractional a => a -> a -> a
/ Double
2.0 forall a. Num a => a -> a -> a
* Double
k2
      forall (m :: * -> *) a. Monad m => a -> m a
return Double
v
    Int
3 -> do
      let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
          n :: Int
n  = Point -> Int
pointIteration Point
p
          ty :: Double
ty = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
0
          t3 :: Double
t3 = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
2
          py :: Point
py = Point
p { pointTime :: Double
pointTime = Double
ty, pointIteration :: Int
pointIteration = Int
n, pointPhase :: Int
pointPhase = Int
0 }
          p3 :: Point
p3 = Point
p { pointTime :: Double
pointTime = Double
t3, pointIteration :: Int
pointIteration = Int
n, pointPhase :: Int
pointPhase = Int
2 }
      Double
vy <- Point -> IO Double
y Point
py
      Double
k3 <- Point -> IO Double
f Point
p3
      let !v :: Double
v = Double
vy forall a. Num a => a -> a -> a
+ Specs -> Double
spcDT Specs
sc forall a. Num a => a -> a -> a
* Double
k3
      forall (m :: * -> *) a. Monad m => a -> m a
return Double
v
    Int
_ -> 
      forall a. HasCallStack => [Char] -> a
error [Char]
"Incorrect phase: integRK4"

integRK4b :: Dynamics Double
             -> Dynamics Double
             -> Dynamics Double
             -> Point -> IO Double
integRK4b :: Dynamics Double
-> Dynamics Double -> Dynamics Double -> Point -> IO Double
integRK4b (Dynamics Point -> IO Double
f) (Dynamics Point -> IO Double
i) (Dynamics Point -> IO Double
y) Point
p =
  case Point -> Int
pointPhase Point
p of
    Int
0 -> case Point -> Int
pointIteration Point
p of
      Int
0 -> 
        Point -> IO Double
i Point
p
      Int
n -> do
        let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
            ty :: Double
ty = Specs -> Int -> Int -> Double
basicTime Specs
sc (Int
n forall a. Num a => a -> a -> a
- Int
1) Int
0
            t1 :: Double
t1 = Double
ty
            t2 :: Double
t2 = Specs -> Int -> Int -> Double
basicTime Specs
sc (Int
n forall a. Num a => a -> a -> a
- Int
1) Int
1
            t3 :: Double
t3 = Specs -> Int -> Int -> Double
basicTime Specs
sc (Int
n forall a. Num a => a -> a -> a
- Int
1) Int
2
            t4 :: Double
t4 = Specs -> Int -> Int -> Double
basicTime Specs
sc (Int
n forall a. Num a => a -> a -> a
- Int
1) Int
3
            py :: Point
py = Point
p { pointTime :: Double
pointTime = Double
ty, pointIteration :: Int
pointIteration = Int
n forall a. Num a => a -> a -> a
- Int
1, pointPhase :: Int
pointPhase = Int
0 }
            p1 :: Point
p1 = Point
py
            p2 :: Point
p2 = Point
p { pointTime :: Double
pointTime = Double
t2, pointIteration :: Int
pointIteration = Int
n forall a. Num a => a -> a -> a
- Int
1, pointPhase :: Int
pointPhase = Int
1 }
            p3 :: Point
p3 = Point
p { pointTime :: Double
pointTime = Double
t3, pointIteration :: Int
pointIteration = Int
n forall a. Num a => a -> a -> a
- Int
1, pointPhase :: Int
pointPhase = Int
2 }
            p4 :: Point
p4 = Point
p { pointTime :: Double
pointTime = Double
t4, pointIteration :: Int
pointIteration = Int
n forall a. Num a => a -> a -> a
- Int
1, pointPhase :: Int
pointPhase = Int
3 }
        Double
vy <- Point -> IO Double
y Point
py
        Double
k1 <- Point -> IO Double
f Point
p1
        Double
k2 <- Point -> IO Double
f Point
p2
        Double
k3 <- Point -> IO Double
f Point
p3
        Double
k4 <- Point -> IO Double
f Point
p4
        let !v :: Double
v = Double
vy forall a. Num a => a -> a -> a
+ Specs -> Double
spcDT Specs
sc forall a. Fractional a => a -> a -> a
/ Double
8.0 forall a. Num a => a -> a -> a
* (Double
k1 forall a. Num a => a -> a -> a
+ Double
3.0 forall a. Num a => a -> a -> a
* (Double
k2 forall a. Num a => a -> a -> a
+ Double
k3) forall a. Num a => a -> a -> a
+ Double
k4)
        forall (m :: * -> *) a. Monad m => a -> m a
return Double
v
    Int
1 -> do
      let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
          n :: Int
n  = Point -> Int
pointIteration Point
p
          ty :: Double
ty = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
0
          t1 :: Double
t1 = Double
ty
          py :: Point
py = Point
p { pointTime :: Double
pointTime = Double
ty, pointIteration :: Int
pointIteration = Int
n, pointPhase :: Int
pointPhase = Int
0 }
          p1 :: Point
p1 = Point
py
      Double
vy <- Point -> IO Double
y Point
py
      Double
k1 <- Point -> IO Double
f Point
p1
      let !v :: Double
v = Double
vy forall a. Num a => a -> a -> a
+ Specs -> Double
spcDT Specs
sc forall a. Fractional a => a -> a -> a
/ Double
3.0 forall a. Num a => a -> a -> a
* Double
k1
      forall (m :: * -> *) a. Monad m => a -> m a
return Double
v
    Int
2 -> do
      let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
          n :: Int
n  = Point -> Int
pointIteration Point
p
          ty :: Double
ty = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
0
          t1 :: Double
t1 = Double
ty
          t2 :: Double
t2 = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
1
          py :: Point
py = Point
p { pointTime :: Double
pointTime = Double
ty, pointIteration :: Int
pointIteration = Int
n, pointPhase :: Int
pointPhase = Int
0 }
          p1 :: Point
p1 = Point
py
          p2 :: Point
p2 = Point
p { pointTime :: Double
pointTime = Double
t2, pointIteration :: Int
pointIteration = Int
n, pointPhase :: Int
pointPhase = Int
1 }
      Double
vy <- Point -> IO Double
y Point
py
      Double
k1 <- Point -> IO Double
f Point
p1
      Double
k2 <- Point -> IO Double
f Point
p2
      let !v :: Double
v = Double
vy forall a. Num a => a -> a -> a
+ Specs -> Double
spcDT Specs
sc forall a. Num a => a -> a -> a
* (- Double
k1 forall a. Fractional a => a -> a -> a
/ Double
3.0 forall a. Num a => a -> a -> a
+ Double
k2)
      forall (m :: * -> *) a. Monad m => a -> m a
return Double
v
    Int
3 -> do
      let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
          n :: Int
n  = Point -> Int
pointIteration Point
p
          ty :: Double
ty = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
0
          t1 :: Double
t1 = Double
ty
          t2 :: Double
t2 = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
1
          t3 :: Double
t3 = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
2
          py :: Point
py = Point
p { pointTime :: Double
pointTime = Double
ty, pointIteration :: Int
pointIteration = Int
n, pointPhase :: Int
pointPhase = Int
0 }
          p1 :: Point
p1 = Point
py
          p2 :: Point
p2 = Point
p { pointTime :: Double
pointTime = Double
t2, pointIteration :: Int
pointIteration = Int
n, pointPhase :: Int
pointPhase = Int
1 }
          p3 :: Point
p3 = Point
p { pointTime :: Double
pointTime = Double
t3, pointIteration :: Int
pointIteration = Int
n, pointPhase :: Int
pointPhase = Int
2 }
      Double
vy <- Point -> IO Double
y Point
py
      Double
k1 <- Point -> IO Double
f Point
p1
      Double
k2 <- Point -> IO Double
f Point
p2
      Double
k3 <- Point -> IO Double
f Point
p3
      let !v :: Double
v = Double
vy forall a. Num a => a -> a -> a
+ Specs -> Double
spcDT Specs
sc forall a. Num a => a -> a -> a
* (Double
k1 forall a. Num a => a -> a -> a
- Double
k2 forall a. Num a => a -> a -> a
+ Double
k3)
      forall (m :: * -> *) a. Monad m => a -> m a
return Double
v
    Int
_ -> 
      forall a. HasCallStack => [Char] -> a
error [Char]
"Incorrect phase: integRK4b"

-- | Return an integral with the specified derivative and initial value.
--
-- To create a loopback, you should use the recursive do-notation.
-- It allows defining the differential equations unordered as
-- in mathematics:
--
-- @
-- model :: Simulation [Double]
-- model = 
--   mdo a <- integ (- ka * a) 100
--       b <- integ (ka * a - kb * b) 0
--       c <- integ (kb * b) 0
--       let ka = 1
--           kb = 1
--       runDynamicsInStopTime $ sequence [a, b, c]
-- @
integ :: Dynamics Double                  -- ^ the derivative
         -> Dynamics Double               -- ^ the initial value
         -> Simulation (Dynamics Double)  -- ^ the integral
integ :: Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
integ Dynamics Double
diff Dynamics Double
i =
  mdo Dynamics Double
y <- forall e. Unboxed e => Dynamics e -> Simulation (Dynamics e)
MU.memoDynamics Dynamics Double
z
      Dynamics Double
z <- forall a. (Run -> IO a) -> Simulation a
Simulation forall a b. (a -> b) -> a -> b
$ \Run
r ->
        case Specs -> Method
spcMethod (Run -> Specs
runSpecs Run
r) of
          Method
Euler -> forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ forall a. (Point -> IO a) -> Dynamics a
Dynamics forall a b. (a -> b) -> a -> b
$ Dynamics Double
-> Dynamics Double -> Dynamics Double -> Point -> IO Double
integEuler Dynamics Double
diff Dynamics Double
i Dynamics Double
y
          Method
RungeKutta2 -> forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ forall a. (Point -> IO a) -> Dynamics a
Dynamics forall a b. (a -> b) -> a -> b
$ Dynamics Double
-> Dynamics Double -> Dynamics Double -> Point -> IO Double
integRK2 Dynamics Double
diff Dynamics Double
i Dynamics Double
y
          Method
RungeKutta4 -> forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ forall a. (Point -> IO a) -> Dynamics a
Dynamics forall a b. (a -> b) -> a -> b
$ Dynamics Double
-> Dynamics Double -> Dynamics Double -> Point -> IO Double
integRK4 Dynamics Double
diff Dynamics Double
i Dynamics Double
y
          Method
RungeKutta4b -> forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ forall a. (Point -> IO a) -> Dynamics a
Dynamics forall a b. (a -> b) -> a -> b
$ Dynamics Double
-> Dynamics Double -> Dynamics Double -> Point -> IO Double
integRK4b Dynamics Double
diff Dynamics Double
i Dynamics Double
y
      forall (m :: * -> *) a. Monad m => a -> m a
return Dynamics Double
y

integEulerEither :: Dynamics (Either Double Double)
                    -> Dynamics Double 
                    -> Dynamics Double 
                    -> Point -> IO Double
integEulerEither :: Dynamics (Either Double Double)
-> Dynamics Double -> Dynamics Double -> Point -> IO Double
integEulerEither (Dynamics Point -> IO (Either Double Double)
f) (Dynamics Point -> IO Double
i) (Dynamics Point -> IO Double
y) Point
p = 
  case Point -> Int
pointIteration Point
p of
    Int
0 -> 
      Point -> IO Double
i Point
p
    Int
n -> do 
      let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
          ty :: Double
ty = Specs -> Int -> Int -> Double
basicTime Specs
sc (Int
n forall a. Num a => a -> a -> a
- Int
1) Int
0
          py :: Point
py = Point
p { pointTime :: Double
pointTime = Double
ty, pointIteration :: Int
pointIteration = Int
n forall a. Num a => a -> a -> a
- Int
1, pointPhase :: Int
pointPhase = Int
0 }
      Either Double Double
b <- Point -> IO (Either Double Double)
f Point
py
      case Either Double Double
b of
        Left Double
v ->
          forall (m :: * -> *) a. Monad m => a -> m a
return Double
v
        Right Double
b -> do
          Double
a <- Point -> IO Double
y Point
py
          let !v :: Double
v = Double
a forall a. Num a => a -> a -> a
+ Specs -> Double
spcDT (Point -> Specs
pointSpecs Point
p) forall a. Num a => a -> a -> a
* Double
b
          forall (m :: * -> *) a. Monad m => a -> m a
return Double
v

-- | Like 'integ' but allows either setting a new 'Left' integral value,
-- or integrating using the 'Right' derivative directly within computation.
--
-- This function always uses Euler's method.
integEither :: Dynamics (Either Double Double)
               -- ^ either set a new 'Left' integral value, or use a 'Right' derivative
               -> Dynamics Double
               -- ^ the initial value
               -> Simulation (Dynamics Double)
integEither :: Dynamics (Either Double Double)
-> Dynamics Double -> Simulation (Dynamics Double)
integEither Dynamics (Either Double Double)
diff Dynamics Double
i =
  mdo Dynamics Double
y <- forall e. Unboxed e => Dynamics e -> Simulation (Dynamics e)
MU.memoDynamics Dynamics Double
z
      Dynamics Double
z <- forall a. (Run -> IO a) -> Simulation a
Simulation forall a b. (a -> b) -> a -> b
$ \Run
r ->
        forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ forall a. (Point -> IO a) -> Dynamics a
Dynamics forall a b. (a -> b) -> a -> b
$ Dynamics (Either Double Double)
-> Dynamics Double -> Dynamics Double -> Point -> IO Double
integEulerEither Dynamics (Either Double Double)
diff Dynamics Double
i Dynamics Double
y
      forall (m :: * -> *) a. Monad m => a -> m a
return Dynamics Double
y

-- | Return the first order exponential smooth.
--
-- To create a loopback, you should use the recursive do-notation
-- with help of which the function itself is defined:
--
-- @
-- smoothI x t i =
--   mdo y <- integ ((x - y) \/ t) i
--       return y
-- @     
smoothI :: Dynamics Double                  -- ^ the value to smooth over time
           -> Dynamics Double               -- ^ time
           -> Dynamics Double               -- ^ the initial value
           -> Simulation (Dynamics Double)  -- ^ the first order exponential smooth
smoothI :: Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
smoothI Dynamics Double
x Dynamics Double
t Dynamics Double
i =
  mdo Dynamics Double
y <- Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
integ ((Dynamics Double
x forall a. Num a => a -> a -> a
- Dynamics Double
y) forall a. Fractional a => a -> a -> a
/ Dynamics Double
t) Dynamics Double
i
      forall (m :: * -> *) a. Monad m => a -> m a
return Dynamics Double
y

-- | Return the first order exponential smooth.
--
-- This is a simplified version of the 'smoothI' function
-- without specifing the initial value.
smooth :: Dynamics Double                  -- ^ the value to smooth over time
          -> Dynamics Double               -- ^ time
          -> Simulation (Dynamics Double)  -- ^ the first order exponential smooth
smooth :: Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
smooth Dynamics Double
x Dynamics Double
t = Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
smoothI Dynamics Double
x Dynamics Double
t Dynamics Double
x

-- | Return the third order exponential smooth.
--
-- To create a loopback, you should use the recursive do-notation
-- with help of which the function itself is defined:
--
-- @
-- smooth3I x t i =
--   mdo y  <- integ ((s2 - y) \/ t') i
--       s2 <- integ ((s1 - s2) \/ t') i
--       s1 <- integ ((x - s1) \/ t') i
--       let t' = t \/ 3.0
--       return y
-- @     
smooth3I :: Dynamics Double                  -- ^ the value to smooth over time
            -> Dynamics Double               -- ^ time
            -> Dynamics Double               -- ^ the initial value
            -> Simulation (Dynamics Double)  -- ^ the third order exponential smooth
smooth3I :: Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
smooth3I Dynamics Double
x Dynamics Double
t Dynamics Double
i =
  mdo Dynamics Double
y  <- Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
integ ((Dynamics Double
s2 forall a. Num a => a -> a -> a
- Dynamics Double
y) forall a. Fractional a => a -> a -> a
/ Dynamics Double
t') Dynamics Double
i
      Dynamics Double
s2 <- Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
integ ((Dynamics Double
s1 forall a. Num a => a -> a -> a
- Dynamics Double
s2) forall a. Fractional a => a -> a -> a
/ Dynamics Double
t') Dynamics Double
i
      Dynamics Double
s1 <- Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
integ ((Dynamics Double
x forall a. Num a => a -> a -> a
- Dynamics Double
s1) forall a. Fractional a => a -> a -> a
/ Dynamics Double
t') Dynamics Double
i
      let t' :: Dynamics Double
t' = Dynamics Double
t forall a. Fractional a => a -> a -> a
/ Dynamics Double
3.0
      forall (m :: * -> *) a. Monad m => a -> m a
return Dynamics Double
y

-- | Return the third order exponential smooth.
-- 
-- This is a simplified version of the 'smooth3I' function
-- without specifying the initial value.
smooth3 :: Dynamics Double                  -- ^ the value to smooth over time
           -> Dynamics Double               -- ^ time
           -> Simulation (Dynamics Double)  -- ^ the third order exponential smooth
smooth3 :: Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
smooth3 Dynamics Double
x Dynamics Double
t = Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
smooth3I Dynamics Double
x Dynamics Double
t Dynamics Double
x

-- | Return the n'th order exponential smooth.
--
-- The result is not discrete in that sense that it may change within the integration time
-- interval depending on the integration method used. Probably, you should apply
-- the 'discreteDynamics' function to the result if you want to achieve an effect when
-- the value is not changed within the time interval, which is used sometimes.
smoothNI :: Dynamics Double                  -- ^ the value to smooth over time
            -> Dynamics Double               -- ^ time
            -> Int                           -- ^ the order
            -> Dynamics Double               -- ^ the initial value
            -> Simulation (Dynamics Double)  -- ^ the n'th order exponential smooth
smoothNI :: Dynamics Double
-> Dynamics Double
-> Int
-> Dynamics Double
-> Simulation (Dynamics Double)
smoothNI Dynamics Double
x Dynamics Double
t Int
n Dynamics Double
i =
  mdo [Dynamics Double]
s <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [Int
1 .. Int
n] forall a b. (a -> b) -> a -> b
$ \Int
k ->
        if Int
k forall a. Eq a => a -> a -> Bool
== Int
1
        then Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
integ ((Dynamics Double
x forall a. Num a => a -> a -> a
- Array Int (Dynamics Double)
a forall i e. Ix i => Array i e -> i -> e
! Int
1) forall a. Fractional a => a -> a -> a
/ Dynamics Double
t') Dynamics Double
i
        else Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
integ ((Array Int (Dynamics Double)
a forall i e. Ix i => Array i e -> i -> e
! (Int
k forall a. Num a => a -> a -> a
- Int
1) forall a. Num a => a -> a -> a
- Array Int (Dynamics Double)
a forall i e. Ix i => Array i e -> i -> e
! Int
k) forall a. Fractional a => a -> a -> a
/ Dynamics Double
t') Dynamics Double
i
      let a :: Array Int (Dynamics Double)
a  = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (Int
1, Int
n) [Dynamics Double]
s 
          t' :: Dynamics Double
t' = Dynamics Double
t forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
      forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ Array Int (Dynamics Double)
a forall i e. Ix i => Array i e -> i -> e
! Int
n

-- | Return the n'th order exponential smooth.
--
-- This is a simplified version of the 'smoothNI' function
-- without specifying the initial value.
smoothN :: Dynamics Double                  -- ^ the value to smooth over time
           -> Dynamics Double               -- ^ time
           -> Int                           -- ^ the order
           -> Simulation (Dynamics Double)  -- ^ the n'th order exponential smooth
smoothN :: Dynamics Double
-> Dynamics Double -> Int -> Simulation (Dynamics Double)
smoothN Dynamics Double
x Dynamics Double
t Int
n = Dynamics Double
-> Dynamics Double
-> Int
-> Dynamics Double
-> Simulation (Dynamics Double)
smoothNI Dynamics Double
x Dynamics Double
t Int
n Dynamics Double
x

-- | Return the first order exponential delay.
--
-- To create a loopback, you should use the recursive do-notation
-- with help of which the function itself is defined:
--
-- @
-- delay1I x t i =
--   mdo y <- integ (x - y \/ t) (i * t)
--       return $ y \/ t
-- @     
delay1I :: Dynamics Double                  -- ^ the value to conserve
           -> Dynamics Double               -- ^ time
           -> Dynamics Double               -- ^ the initial value
           -> Simulation (Dynamics Double)  -- ^ the first order exponential delay
delay1I :: Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
delay1I Dynamics Double
x Dynamics Double
t Dynamics Double
i =
  mdo Dynamics Double
y <- Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
integ (Dynamics Double
x forall a. Num a => a -> a -> a
- Dynamics Double
y forall a. Fractional a => a -> a -> a
/ Dynamics Double
t) (Dynamics Double
i forall a. Num a => a -> a -> a
* Dynamics Double
t)
      forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ Dynamics Double
y forall a. Fractional a => a -> a -> a
/ Dynamics Double
t

-- | Return the first order exponential delay.
--
-- This is a simplified version of the 'delay1I' function
-- without specifying the initial value.
delay1 :: Dynamics Double                  -- ^ the value to conserve
          -> Dynamics Double               -- ^ time
          -> Simulation (Dynamics Double)  -- ^ the first order exponential delay
delay1 :: Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
delay1 Dynamics Double
x Dynamics Double
t = Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
delay1I Dynamics Double
x Dynamics Double
t Dynamics Double
x

-- | Return the third order exponential delay.
delay3I :: Dynamics Double                  -- ^ the value to conserve
           -> Dynamics Double               -- ^ time
           -> Dynamics Double               -- ^ the initial value
           -> Simulation (Dynamics Double)  -- ^ the third order exponential delay
delay3I :: Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
delay3I Dynamics Double
x Dynamics Double
t Dynamics Double
i =
  mdo Dynamics Double
y  <- Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
integ (Dynamics Double
s2 forall a. Fractional a => a -> a -> a
/ Dynamics Double
t' forall a. Num a => a -> a -> a
- Dynamics Double
y forall a. Fractional a => a -> a -> a
/ Dynamics Double
t') (Dynamics Double
i forall a. Num a => a -> a -> a
* Dynamics Double
t')
      Dynamics Double
s2 <- Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
integ (Dynamics Double
s1 forall a. Fractional a => a -> a -> a
/ Dynamics Double
t' forall a. Num a => a -> a -> a
- Dynamics Double
s2 forall a. Fractional a => a -> a -> a
/ Dynamics Double
t') (Dynamics Double
i forall a. Num a => a -> a -> a
* Dynamics Double
t')
      Dynamics Double
s1 <- Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
integ (Dynamics Double
x forall a. Num a => a -> a -> a
- Dynamics Double
s1 forall a. Fractional a => a -> a -> a
/ Dynamics Double
t') (Dynamics Double
i forall a. Num a => a -> a -> a
* Dynamics Double
t')
      let t' :: Dynamics Double
t' = Dynamics Double
t forall a. Fractional a => a -> a -> a
/ Dynamics Double
3.0
      forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ Dynamics Double
y forall a. Fractional a => a -> a -> a
/ Dynamics Double
t'         

-- | Return the third order exponential delay.
--
-- This is a simplified version of the 'delay3I' function
-- without specifying the initial value.
delay3 :: Dynamics Double                  -- ^ the value to conserve
          -> Dynamics Double               -- ^ time
          -> Simulation (Dynamics Double)  -- ^ the third order exponential delay
delay3 :: Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
delay3 Dynamics Double
x Dynamics Double
t = Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
delay3I Dynamics Double
x Dynamics Double
t Dynamics Double
x

-- | Return the n'th order exponential delay.
delayNI :: Dynamics Double                  -- ^ the value to conserve
           -> Dynamics Double               -- ^ time
           -> Int                           -- ^ the order
           -> Dynamics Double               -- ^ the initial value
           -> Simulation (Dynamics Double)  -- ^ the n'th order exponential delay
delayNI :: Dynamics Double
-> Dynamics Double
-> Int
-> Dynamics Double
-> Simulation (Dynamics Double)
delayNI Dynamics Double
x Dynamics Double
t Int
n Dynamics Double
i =
  mdo [Dynamics Double]
s <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [Int
1 .. Int
n] forall a b. (a -> b) -> a -> b
$ \Int
k ->
        if Int
k forall a. Eq a => a -> a -> Bool
== Int
1
        then Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
integ (Dynamics Double
x forall a. Num a => a -> a -> a
- (Array Int (Dynamics Double)
a forall i e. Ix i => Array i e -> i -> e
! Int
1) forall a. Fractional a => a -> a -> a
/ Dynamics Double
t') (Dynamics Double
i forall a. Num a => a -> a -> a
* Dynamics Double
t')
        else Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
integ ((Array Int (Dynamics Double)
a forall i e. Ix i => Array i e -> i -> e
! (Int
k forall a. Num a => a -> a -> a
- Int
1)) forall a. Fractional a => a -> a -> a
/ Dynamics Double
t' forall a. Num a => a -> a -> a
- (Array Int (Dynamics Double)
a forall i e. Ix i => Array i e -> i -> e
! Int
k) forall a. Fractional a => a -> a -> a
/ Dynamics Double
t') (Dynamics Double
i forall a. Num a => a -> a -> a
* Dynamics Double
t')
      let a :: Array Int (Dynamics Double)
a  = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (Int
1, Int
n) [Dynamics Double]
s
          t' :: Dynamics Double
t' = Dynamics Double
t forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
      forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ (Array Int (Dynamics Double)
a forall i e. Ix i => Array i e -> i -> e
! Int
n) forall a. Fractional a => a -> a -> a
/ Dynamics Double
t'

-- | Return the n'th order exponential delay.
--
-- This is a simplified version of the 'delayNI' function
-- without specifying the initial value.
delayN :: Dynamics Double                  -- ^ the value to conserve
          -> Dynamics Double               -- ^ time
          -> Int                           -- ^ the order
          -> Simulation (Dynamics Double)  -- ^ the n'th order exponential delay
delayN :: Dynamics Double
-> Dynamics Double -> Int -> Simulation (Dynamics Double)
delayN Dynamics Double
x Dynamics Double
t Int
n = Dynamics Double
-> Dynamics Double
-> Int
-> Dynamics Double
-> Simulation (Dynamics Double)
delayNI Dynamics Double
x Dynamics Double
t Int
n Dynamics Double
x

-- | Return the forecast.
--
-- The function has the following definition:
--
-- @
-- forecast x at hz =
--   do y <- smooth x at
--      return $ x * (1.0 + (x \/ y - 1.0) \/ at * hz)
-- @
forecast :: Dynamics Double                  -- ^ the value to forecast
            -> Dynamics Double               -- ^ the average time
            -> Dynamics Double               -- ^ the time horizon
            -> Simulation (Dynamics Double)  -- ^ the forecast
forecast :: Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
forecast Dynamics Double
x Dynamics Double
at Dynamics Double
hz =
  do Dynamics Double
y <- Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
smooth Dynamics Double
x Dynamics Double
at
     forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ Dynamics Double
x forall a. Num a => a -> a -> a
* (Dynamics Double
1.0 forall a. Num a => a -> a -> a
+ (Dynamics Double
x forall a. Fractional a => a -> a -> a
/ Dynamics Double
y forall a. Num a => a -> a -> a
- Dynamics Double
1.0) forall a. Fractional a => a -> a -> a
/ Dynamics Double
at forall a. Num a => a -> a -> a
* Dynamics Double
hz)

-- | Return the trend.
--
-- The function has the following definition:
--
-- @
-- trend x at i =
--   do y <- smoothI x at (x \/ (1.0 + i * at))
--      return $ (x \/ y - 1.0) \/ at
-- @
trend :: Dynamics Double                  -- ^ the value for which the trend is calculated
         -> Dynamics Double               -- ^ the average time
         -> Dynamics Double               -- ^ the initial value
         -> Simulation (Dynamics Double)  -- ^ the fractional change rate
trend :: Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
trend Dynamics Double
x Dynamics Double
at Dynamics Double
i =
  do Dynamics Double
y <- Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
smoothI Dynamics Double
x Dynamics Double
at (Dynamics Double
x forall a. Fractional a => a -> a -> a
/ (Dynamics Double
1.0 forall a. Num a => a -> a -> a
+ Dynamics Double
i forall a. Num a => a -> a -> a
* Dynamics Double
at))
     forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ (Dynamics Double
x forall a. Fractional a => a -> a -> a
/ Dynamics Double
y forall a. Num a => a -> a -> a
- Dynamics Double
1.0) forall a. Fractional a => a -> a -> a
/ Dynamics Double
at

--
-- Difference Equations
--

-- | Retun the sum for the difference equation.
-- It is like an integral returned by the 'integ' function, only now
-- the difference is used instead of derivative.
--
-- As usual, to create a loopback, you should use the recursive do-notation.
diffsum :: (Num a, Unboxed a)
           => Dynamics a               -- ^ the difference
           -> Dynamics a               -- ^ the initial value
           -> Simulation (Dynamics a)  -- ^ the sum
diffsum :: forall a.
(Num a, Unboxed a) =>
Dynamics a -> Dynamics a -> Simulation (Dynamics a)
diffsum (Dynamics Point -> IO a
diff) (Dynamics Point -> IO a
i) =
  mdo Dynamics a
y <-
        forall e. Unboxed e => Dynamics e -> Simulation (Dynamics e)
MU.memo0Dynamics forall a b. (a -> b) -> a -> b
$
        forall a. (Point -> IO a) -> Dynamics a
Dynamics forall a b. (a -> b) -> a -> b
$ \Point
p ->
        case Point -> Int
pointIteration Point
p of
          Int
0 -> Point -> IO a
i Point
p
          Int
n -> do 
            let Dynamics Point -> IO a
m = Dynamics a
y
                sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
                ty :: Double
ty = Specs -> Int -> Int -> Double
basicTime Specs
sc (Int
n forall a. Num a => a -> a -> a
- Int
1) Int
0
                py :: Point
py = Point
p { pointTime :: Double
pointTime = Double
ty, 
                         pointIteration :: Int
pointIteration = Int
n forall a. Num a => a -> a -> a
- Int
1, 
                         pointPhase :: Int
pointPhase = Int
0 }
            a
a <- Point -> IO a
m Point
py
            a
b <- Point -> IO a
diff Point
py
            let !v :: a
v = a
a forall a. Num a => a -> a -> a
+ a
b
            forall (m :: * -> *) a. Monad m => a -> m a
return a
v
      forall (m :: * -> *) a. Monad m => a -> m a
return Dynamics a
y

-- | Like 'diffsum' but allows either setting a new 'Left' sum value, or adding the 'Right' difference.
diffsumEither :: (Num a, Unboxed a)
                 => Dynamics (Either a a)
                 -- ^ either set the 'Left' value for the sum, or add the 'Right' difference to the sum
                 -> Dynamics a
                 -- ^ the initial value
                 -> Simulation (Dynamics a)
                 -- ^ the sum
diffsumEither :: forall a.
(Num a, Unboxed a) =>
Dynamics (Either a a) -> Dynamics a -> Simulation (Dynamics a)
diffsumEither (Dynamics Point -> IO (Either a a)
diff) (Dynamics Point -> IO a
i) =
  mdo Dynamics a
y <-
        forall e. Unboxed e => Dynamics e -> Simulation (Dynamics e)
MU.memo0Dynamics forall a b. (a -> b) -> a -> b
$
        forall a. (Point -> IO a) -> Dynamics a
Dynamics forall a b. (a -> b) -> a -> b
$ \Point
p ->
        case Point -> Int
pointIteration Point
p of
          Int
0 -> Point -> IO a
i Point
p
          Int
n -> do 
            let Dynamics Point -> IO a
m = Dynamics a
y
                sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
                ty :: Double
ty = Specs -> Int -> Int -> Double
basicTime Specs
sc (Int
n forall a. Num a => a -> a -> a
- Int
1) Int
0
                py :: Point
py = Point
p { pointTime :: Double
pointTime = Double
ty, 
                         pointIteration :: Int
pointIteration = Int
n forall a. Num a => a -> a -> a
- Int
1, 
                         pointPhase :: Int
pointPhase = Int
0 }
            Either a a
b <- Point -> IO (Either a a)
diff Point
py
            case Either a a
b of
              Left a
v ->
                forall (m :: * -> *) a. Monad m => a -> m a
return a
v
              Right a
b -> do
                a
a <- Point -> IO a
m Point
py
                let !v :: a
v = a
a forall a. Num a => a -> a -> a
+ a
b
                forall (m :: * -> *) a. Monad m => a -> m a
return a
v
      forall (m :: * -> *) a. Monad m => a -> m a
return Dynamics a
y

--
-- Table Functions
--

-- | Lookup @x@ in a table of pairs @(x, y)@ using linear interpolation.
lookupDynamics :: Dynamics Double -> Array Int (Double, Double) -> Dynamics Double
lookupDynamics :: Dynamics Double -> Array Int (Double, Double) -> Dynamics Double
lookupDynamics (Dynamics Point -> IO Double
m) Array Int (Double, Double)
tbl =
  forall a. (Point -> IO a) -> Dynamics a
Dynamics forall a b. (a -> b) -> a -> b
$ \Point
p ->
  do Double
a <- Point -> IO Double
m Point
p
     forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ Double -> Array Int (Double, Double) -> Double
tableLookup Double
a Array Int (Double, Double)
tbl

-- | Lookup @x@ in a table of pairs @(x, y)@ using stepwise function.
lookupStepwiseDynamics :: Dynamics Double -> Array Int (Double, Double) -> Dynamics Double
lookupStepwiseDynamics :: Dynamics Double -> Array Int (Double, Double) -> Dynamics Double
lookupStepwiseDynamics (Dynamics Point -> IO Double
m) Array Int (Double, Double)
tbl =
  forall a. (Point -> IO a) -> Dynamics a
Dynamics forall a b. (a -> b) -> a -> b
$ \Point
p ->
  do Double
a <- Point -> IO Double
m Point
p
     forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ Double -> Array Int (Double, Double) -> Double
tableLookupStepwise Double
a Array Int (Double, Double)
tbl

--
-- Discrete Functions
--

-- | Return the delayed value using the specified lag time.
-- This function is less accurate than 'delayByDT'.
delay :: Dynamics a          -- ^ the value to delay
         -> Dynamics Double  -- ^ the lag time
         -> Dynamics a       -- ^ the delayed value
delay :: forall a. Dynamics a -> Dynamics Double -> Dynamics a
delay (Dynamics Point -> IO a
x) (Dynamics Point -> IO Double
d) = forall a. Dynamics a -> Dynamics a
discreteDynamics forall a b. (a -> b) -> a -> b
$ forall a. (Point -> IO a) -> Dynamics a
Dynamics Point -> IO a
r 
  where
    r :: Point -> IO a
r Point
p = do 
      let t :: Double
t  = Point -> Double
pointTime Point
p
          sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
          n :: Int
n  = Point -> Int
pointIteration Point
p
      Double
a <- Point -> IO Double
d Point
p
      let t' :: Double
t' = Double
t forall a. Num a => a -> a -> a
- Double
a
          n' :: Int
n' = forall a b. (Integral a, Num b) => a -> b
fromIntegral forall a b. (a -> b) -> a -> b
$ forall a b. (RealFrac a, Integral b) => a -> b
floor forall a b. (a -> b) -> a -> b
$ (Double
t' forall a. Num a => a -> a -> a
- Specs -> Double
spcStartTime Specs
sc) forall a. Fractional a => a -> a -> a
/ Specs -> Double
spcDT Specs
sc
          y :: IO a
y | Int
n' forall a. Ord a => a -> a -> Bool
< Int
0    = Point -> IO a
x forall a b. (a -> b) -> a -> b
$ Point
p { pointTime :: Double
pointTime = Specs -> Double
spcStartTime Specs
sc,
                                  pointIteration :: Int
pointIteration = Int
0, 
                                  pointPhase :: Int
pointPhase = Int
0 }
            | Int
n' forall a. Ord a => a -> a -> Bool
< Int
n    = Point -> IO a
x forall a b. (a -> b) -> a -> b
$ Point
p { pointTime :: Double
pointTime = Double
t',
                                  pointIteration :: Int
pointIteration = Int
n',
                                  pointPhase :: Int
pointPhase = -Int
1 }
            | Int
n' forall a. Ord a => a -> a -> Bool
> Int
n    = forall a. HasCallStack => [Char] -> a
error forall a b. (a -> b) -> a -> b
$
                          [Char]
"Cannot return the future data: delay. " forall a. [a] -> [a] -> [a]
++
                          [Char]
"The lag time cannot be negative."
            | Bool
otherwise = forall a. HasCallStack => [Char] -> a
error forall a b. (a -> b) -> a -> b
$
                          [Char]
"Cannot return the current data: delay. " forall a. [a] -> [a] -> [a]
++
                          [Char]
"The lag time is too small."
      IO a
y

-- | Return the delayed value using the specified lag time and initial value.
-- Because of the latter, it allows creating a loop back.
-- This function is less accurate than 'delayIByDT'.
delayI :: Dynamics a          -- ^ the value to delay
          -> Dynamics Double  -- ^ the lag time
          -> Dynamics a       -- ^ the initial value
          -> Simulation (Dynamics a)    -- ^ the delayed value
delayI :: forall a.
Dynamics a
-> Dynamics Double -> Dynamics a -> Simulation (Dynamics a)
delayI (Dynamics Point -> IO a
x) (Dynamics Point -> IO Double
d) (Dynamics Point -> IO a
i) = forall e. Dynamics e -> Simulation (Dynamics e)
M.memo0Dynamics forall a b. (a -> b) -> a -> b
$ forall a. (Point -> IO a) -> Dynamics a
Dynamics Point -> IO a
r 
  where
    r :: Point -> IO a
r Point
p = do 
      let t :: Double
t  = Point -> Double
pointTime Point
p
          sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
          n :: Int
n  = Point -> Int
pointIteration Point
p
      Double
a <- Point -> IO Double
d Point
p
      let t' :: Double
t' = Double
t forall a. Num a => a -> a -> a
- Double
a
          n' :: Int
n' = forall a b. (Integral a, Num b) => a -> b
fromIntegral forall a b. (a -> b) -> a -> b
$ forall a b. (RealFrac a, Integral b) => a -> b
floor forall a b. (a -> b) -> a -> b
$ (Double
t' forall a. Num a => a -> a -> a
- Specs -> Double
spcStartTime Specs
sc) forall a. Fractional a => a -> a -> a
/ Specs -> Double
spcDT Specs
sc
          y :: IO a
y | Int
n' forall a. Ord a => a -> a -> Bool
< Int
0    = Point -> IO a
i forall a b. (a -> b) -> a -> b
$ Point
p { pointTime :: Double
pointTime = Specs -> Double
spcStartTime Specs
sc,
                                  pointIteration :: Int
pointIteration = Int
0, 
                                  pointPhase :: Int
pointPhase = Int
0 }
            | Int
n' forall a. Ord a => a -> a -> Bool
< Int
n    = Point -> IO a
x forall a b. (a -> b) -> a -> b
$ Point
p { pointTime :: Double
pointTime = Double
t',
                                  pointIteration :: Int
pointIteration = Int
n',
                                  pointPhase :: Int
pointPhase = -Int
1 }
            | Int
n' forall a. Ord a => a -> a -> Bool
> Int
n    = forall a. HasCallStack => [Char] -> a
error forall a b. (a -> b) -> a -> b
$
                          [Char]
"Cannot return the future data: delayI. " forall a. [a] -> [a] -> [a]
++
                          [Char]
"The lag time cannot be negative."
            | Bool
otherwise = forall a. HasCallStack => [Char] -> a
error forall a b. (a -> b) -> a -> b
$
                          [Char]
"Cannot return the current data: delayI. " forall a. [a] -> [a] -> [a]
++
                          [Char]
"The lag time is too small."
      IO a
y

-- | Return the delayed value by the specified positive number of
-- integration time steps used for calculating the lag time.
delayByDT :: Dynamics a
             -- ^ the value to delay
             -> Dynamics Int
             -- ^ the delay as a multiplication of the corresponding number
             -- and the integration time step
             -> Dynamics a
             -- ^ the delayed value
delayByDT :: forall a. Dynamics a -> Dynamics Int -> Dynamics a
delayByDT (Dynamics Point -> IO a
x) (Dynamics Point -> IO Int
d) = forall a. Dynamics a -> Dynamics a
discreteDynamics forall a b. (a -> b) -> a -> b
$ forall a. (Point -> IO a) -> Dynamics a
Dynamics Point -> IO a
r 
  where
    r :: Point -> IO a
r Point
p = do 
      let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
          n :: Int
n  = Point -> Int
pointIteration Point
p
      Int
a <- Point -> IO Int
d Point
p
      let p' :: Point
p' = Point -> Int -> Point
delayPoint Point
p Int
a
          n' :: Int
n' = Point -> Int
pointIteration Point
p'
          y :: IO a
y | Int
n' forall a. Ord a => a -> a -> Bool
< Int
0    = Point -> IO a
x forall a b. (a -> b) -> a -> b
$ Point
p { pointTime :: Double
pointTime = Specs -> Double
spcStartTime Specs
sc,
                                  pointIteration :: Int
pointIteration = Int
0, 
                                  pointPhase :: Int
pointPhase = Int
0 }
            | Int
n' forall a. Ord a => a -> a -> Bool
< Int
n    = Point -> IO a
x Point
p'
            | Int
n' forall a. Ord a => a -> a -> Bool
> Int
n    = forall a. HasCallStack => [Char] -> a
error forall a b. (a -> b) -> a -> b
$
                          [Char]
"Cannot return the future data: delayByDT. " forall a. [a] -> [a] -> [a]
++
                          [Char]
"The lag time cannot be negative."
            | Bool
otherwise = forall a. HasCallStack => [Char] -> a
error forall a b. (a -> b) -> a -> b
$
                          [Char]
"Cannot return the current data: delayByDT. " forall a. [a] -> [a] -> [a]
++
                          [Char]
"The lag time is too small."
      IO a
y
      
-- | Return the delayed value by the specified initial value and
-- a positive number of integration time steps used for calculating
-- the lag time. It allows creating a loop back.
delayIByDT :: Dynamics a
              -- ^ the value to delay
              -> Dynamics Int
              -- ^ the delay as a multiplication of the corresponding number
              -- and the integration time step
              -> Dynamics a
              -- ^ the initial value
              -> Simulation (Dynamics a)
              -- ^ the delayed value
delayIByDT :: forall a.
Dynamics a -> Dynamics Int -> Dynamics a -> Simulation (Dynamics a)
delayIByDT (Dynamics Point -> IO a
x) (Dynamics Point -> IO Int
d) (Dynamics Point -> IO a
i) = forall e. Dynamics e -> Simulation (Dynamics e)
M.memoDynamics forall a b. (a -> b) -> a -> b
$ forall a. (Point -> IO a) -> Dynamics a
Dynamics Point -> IO a
r 
  where
    r :: Point -> IO a
r Point
p = do 
      let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
          n :: Int
n  = Point -> Int
pointIteration Point
p
      Int
a <- Point -> IO Int
d Point
p
      let p' :: Point
p' = Point -> Int -> Point
delayPoint Point
p Int
a
          n' :: Int
n' = Point -> Int
pointIteration Point
p'
          y :: IO a
y | Int
n' forall a. Ord a => a -> a -> Bool
< Int
0    = Point -> IO a
i forall a b. (a -> b) -> a -> b
$ Point
p { pointTime :: Double
pointTime = Specs -> Double
spcStartTime Specs
sc,
                                  pointIteration :: Int
pointIteration = Int
0, 
                                  pointPhase :: Int
pointPhase = Int
0 }
            | Int
n' forall a. Ord a => a -> a -> Bool
< Int
n    = Point -> IO a
x Point
p'
            | Int
n' forall a. Ord a => a -> a -> Bool
> Int
n    = forall a. HasCallStack => [Char] -> a
error forall a b. (a -> b) -> a -> b
$
                          [Char]
"Cannot return the future data: delayIByDT. " forall a. [a] -> [a] -> [a]
++
                          [Char]
"The lag time cannot be negative."
            | Bool
otherwise = forall a. HasCallStack => [Char] -> a
error forall a b. (a -> b) -> a -> b
$
                          [Char]
"Cannot return the current data: delayIByDT. " forall a. [a] -> [a] -> [a]
++
                          [Char]
"The lag time is too small."
      IO a
y

--
-- Financial Functions
--

-- | Return the Net Present Value (NPV) of the stream computed using the specified
-- discount rate, the initial value and some factor (usually 1).
--
-- It is defined in the following way:
--
-- @
-- npv stream rate init factor =
--   mdo let dt' = liftParameter dt
--       df <- integ (- df * rate) 1
--       accum <- integ (stream * df) init
--       return $ (accum + dt' * stream * df) * factor
-- @
npv :: Dynamics Double                  -- ^ the stream
       -> Dynamics Double               -- ^ the discount rate
       -> Dynamics Double               -- ^ the initial value
       -> Dynamics Double               -- ^ factor
       -> Simulation (Dynamics Double)  -- ^ the Net Present Value (NPV)
npv :: Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
npv Dynamics Double
stream Dynamics Double
rate Dynamics Double
init Dynamics Double
factor =
  mdo let dt' :: Dynamics Double
dt' = forall (m :: * -> *) a. ParameterLift m => Parameter a -> m a
liftParameter Parameter Double
dt
      Dynamics Double
df <- Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
integ (- Dynamics Double
df forall a. Num a => a -> a -> a
* Dynamics Double
rate) Dynamics Double
1
      Dynamics Double
accum <- Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
integ (Dynamics Double
stream forall a. Num a => a -> a -> a
* Dynamics Double
df) Dynamics Double
init
      forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ (Dynamics Double
accum forall a. Num a => a -> a -> a
+ Dynamics Double
dt' forall a. Num a => a -> a -> a
* Dynamics Double
stream forall a. Num a => a -> a -> a
* Dynamics Double
df) forall a. Num a => a -> a -> a
* Dynamics Double
factor

-- | Return the Net Present Value End of period (NPVE) of the stream computed
-- using the specified discount rate, the initial value and some factor.
--
-- It is defined in the following way:
--
-- @
-- npve stream rate init factor =
--   mdo let dt' = liftParameter dt
--       df <- integ (- df * rate \/ (1 + rate * dt')) (1 \/ (1 + rate * dt'))
--       accum <- integ (stream * df) init
--       return $ (accum + dt' * stream * df) * factor
-- @
npve :: Dynamics Double                  -- ^ the stream
        -> Dynamics Double               -- ^ the discount rate
        -> Dynamics Double               -- ^ the initial value
        -> Dynamics Double               -- ^ factor
        -> Simulation (Dynamics Double)  -- ^ the Net Present Value End (NPVE)
npve :: Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Dynamics Double
-> Simulation (Dynamics Double)
npve Dynamics Double
stream Dynamics Double
rate Dynamics Double
init Dynamics Double
factor =
  mdo let dt' :: Dynamics Double
dt' = forall (m :: * -> *) a. ParameterLift m => Parameter a -> m a
liftParameter Parameter Double
dt
      Dynamics Double
df <- Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
integ (- Dynamics Double
df forall a. Num a => a -> a -> a
* Dynamics Double
rate forall a. Fractional a => a -> a -> a
/ (Dynamics Double
1 forall a. Num a => a -> a -> a
+ Dynamics Double
rate forall a. Num a => a -> a -> a
* Dynamics Double
dt')) (Dynamics Double
1 forall a. Fractional a => a -> a -> a
/ (Dynamics Double
1 forall a. Num a => a -> a -> a
+ Dynamics Double
rate forall a. Num a => a -> a -> a
* Dynamics Double
dt'))
      Dynamics Double
accum <- Dynamics Double -> Dynamics Double -> Simulation (Dynamics Double)
integ (Dynamics Double
stream forall a. Num a => a -> a -> a
* Dynamics Double
df) Dynamics Double
init
      forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ (Dynamics Double
accum forall a. Num a => a -> a -> a
+ Dynamics Double
dt' forall a. Num a => a -> a -> a
* Dynamics Double
stream forall a. Num a => a -> a -> a
* Dynamics Double
df) forall a. Num a => a -> a -> a
* Dynamics Double
factor

-- | Computation that returns 0 until the step time and then returns the specified height.
step :: Dynamics Double
        -- ^ the height
        -> Dynamics Double
        -- ^ the step time
        -> Dynamics Double
step :: Dynamics Double -> Dynamics Double -> Dynamics Double
step Dynamics Double
h Dynamics Double
st =
  forall a. Dynamics a -> Dynamics a
discreteDynamics forall a b. (a -> b) -> a -> b
$
  forall a. (Point -> IO a) -> Dynamics a
Dynamics forall a b. (a -> b) -> a -> b
$ \Point
p ->
  do let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
         t :: Double
t  = Point -> Double
pointTime Point
p
     Double
st' <- forall a. Point -> Dynamics a -> IO a
invokeDynamics Point
p Dynamics Double
st
     let t' :: Double
t' = Double
t forall a. Num a => a -> a -> a
+ Specs -> Double
spcDT Specs
sc forall a. Fractional a => a -> a -> a
/ Double
2
     if Double
st' forall a. Ord a => a -> a -> Bool
< Double
t'
       then forall a. Point -> Dynamics a -> IO a
invokeDynamics Point
p Dynamics Double
h
       else forall (m :: * -> *) a. Monad m => a -> m a
return Double
0

-- | Computation that returns 1, starting at the time start, and lasting for the interval
-- width; 0 is returned at all other times.
pulse :: Dynamics Double
         -- ^ the time start
         -> Dynamics Double
         -- ^ the interval width
         -> Dynamics Double
pulse :: Dynamics Double -> Dynamics Double -> Dynamics Double
pulse Dynamics Double
st Dynamics Double
w =
  forall a. Dynamics a -> Dynamics a
discreteDynamics forall a b. (a -> b) -> a -> b
$
  forall a. (Point -> IO a) -> Dynamics a
Dynamics forall a b. (a -> b) -> a -> b
$ \Point
p ->
  do let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
         t :: Double
t  = Point -> Double
pointTime Point
p
     Double
st' <- forall a. Point -> Dynamics a -> IO a
invokeDynamics Point
p Dynamics Double
st
     let t' :: Double
t' = Double
t forall a. Num a => a -> a -> a
+ Specs -> Double
spcDT Specs
sc forall a. Fractional a => a -> a -> a
/ Double
2
     if Double
st' forall a. Ord a => a -> a -> Bool
< Double
t'
       then do Double
w' <- forall a. Point -> Dynamics a -> IO a
invokeDynamics Point
p Dynamics Double
w
               forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ if Double
t' forall a. Ord a => a -> a -> Bool
< Double
st' forall a. Num a => a -> a -> a
+ Double
w' then Double
1 else Double
0
       else forall (m :: * -> *) a. Monad m => a -> m a
return Double
0

-- | Computation that returns 1, starting at the time start, and lasting for the interval
-- width and then repeats this pattern with the specified period; 0 is returned at all
-- other times.
pulseP :: Dynamics Double
          -- ^ the time start
          -> Dynamics Double
          -- ^ the interval width
          -> Dynamics Double
          -- ^ the time period
          -> Dynamics Double
pulseP :: Dynamics Double
-> Dynamics Double -> Dynamics Double -> Dynamics Double
pulseP Dynamics Double
st Dynamics Double
w Dynamics Double
period =
  forall a. Dynamics a -> Dynamics a
discreteDynamics forall a b. (a -> b) -> a -> b
$
  forall a. (Point -> IO a) -> Dynamics a
Dynamics forall a b. (a -> b) -> a -> b
$ \Point
p ->
  do let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
         t :: Double
t  = Point -> Double
pointTime Point
p
     Double
p'  <- forall a. Point -> Dynamics a -> IO a
invokeDynamics Point
p Dynamics Double
period
     Double
st' <- forall a. Point -> Dynamics a -> IO a
invokeDynamics Point
p Dynamics Double
st
     let y' :: Double
y' = if (Double
p' forall a. Ord a => a -> a -> Bool
> Double
0) Bool -> Bool -> Bool
&& (Double
t forall a. Ord a => a -> a -> Bool
> Double
st')
              then forall a b. (Integral a, Num b) => a -> b
fromIntegral (forall a b. (RealFrac a, Integral b) => a -> b
floor forall a b. (a -> b) -> a -> b
$ (Double
t forall a. Num a => a -> a -> a
- Double
st') forall a. Fractional a => a -> a -> a
/ Double
p') forall a. Num a => a -> a -> a
* Double
p'
              else Double
0
     let st' :: Double
st' = Double
st' forall a. Num a => a -> a -> a
+ Double
y'
     let t' :: Double
t' = Double
t forall a. Num a => a -> a -> a
+ Specs -> Double
spcDT Specs
sc forall a. Fractional a => a -> a -> a
/ Double
2
     if Double
st' forall a. Ord a => a -> a -> Bool
< Double
t'
       then do Double
w' <- forall a. Point -> Dynamics a -> IO a
invokeDynamics Point
p Dynamics Double
w
               forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ if Double
t' forall a. Ord a => a -> a -> Bool
< Double
st' forall a. Num a => a -> a -> a
+ Double
w' then Double
1 else Double
0
       else forall (m :: * -> *) a. Monad m => a -> m a
return Double
0

-- | Computation that returns 0 until the specified time start and then
-- slopes upward until the end time and then holds constant.
ramp :: Dynamics Double
        -- ^ the slope parameter
        -> Dynamics Double
        -- ^ the time start
        -> Dynamics Double
        -- ^ the end time
        -> Dynamics Double
ramp :: Dynamics Double
-> Dynamics Double -> Dynamics Double -> Dynamics Double
ramp Dynamics Double
slope Dynamics Double
st Dynamics Double
e =
  forall a. Dynamics a -> Dynamics a
discreteDynamics forall a b. (a -> b) -> a -> b
$
  forall a. (Point -> IO a) -> Dynamics a
Dynamics forall a b. (a -> b) -> a -> b
$ \Point
p ->
  do let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
         t :: Double
t  = Point -> Double
pointTime Point
p
     Double
st' <- forall a. Point -> Dynamics a -> IO a
invokeDynamics Point
p Dynamics Double
st
     if Double
st' forall a. Ord a => a -> a -> Bool
< Double
t
       then do Double
slope' <- forall a. Point -> Dynamics a -> IO a
invokeDynamics Point
p Dynamics Double
slope
               Double
e' <- forall a. Point -> Dynamics a -> IO a
invokeDynamics Point
p Dynamics Double
e
               if Double
t forall a. Ord a => a -> a -> Bool
< Double
e'
                 then forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ Double
slope' forall a. Num a => a -> a -> a
* (Double
t forall a. Num a => a -> a -> a
- Double
st')
                 else forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ Double
slope' forall a. Num a => a -> a -> a
* (Double
e' forall a. Num a => a -> a -> a
- Double
st')
       else forall (m :: * -> *) a. Monad m => a -> m a
return Double
0