```
{-# LANGUAGE BangPatterns, RecursiveDo, FlexibleContexts #-}

-- |
-- Module     : Simulation.Aivika.Trans.SystemDynamics
-- 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.Trans.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 Simulation.Aivika.Trans.Internal.Specs
import Simulation.Aivika.Trans.Internal.Parameter
import Simulation.Aivika.Trans.Internal.Simulation
import Simulation.Aivika.Trans.Internal.Dynamics
import Simulation.Aivika.Trans.Dynamics.Extra
import Simulation.Aivika.Trans.Table
import Simulation.Aivika.Trans.SD

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

--
-- Equality and Ordering
--

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

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

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

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

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

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

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

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

-- | Implement the if-then-else operator.
ifDynamics :: Monad m => Dynamics m Bool -> Dynamics m a -> Dynamics m a -> Dynamics m a
{-# INLINE ifDynamics #-}
ifDynamics :: Dynamics m Bool -> Dynamics m a -> Dynamics m a -> Dynamics m a
ifDynamics Dynamics m Bool
cond Dynamics m a
x Dynamics m a
y =
do Bool
a <- Dynamics m Bool
cond
if Bool
a then Dynamics m a
x else Dynamics m a
y

--
-- Ordinary Differential Equations
--

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

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

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

=> Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Point m
-> m Double
{-# INLINABLE integRK4b #-}
integRK4b :: Dynamics m Double
-> Dynamics m Double -> Dynamics m Double -> Point m -> m Double
integRK4b (Dynamics Point m -> m Double
f) (Dynamics Point m -> m Double
i) (Dynamics Point m -> m Double
y) Point m
p =
case Point m -> Int
forall (m :: * -> *). Point m -> Int
pointPhase Point m
p of
Int
0 -> case Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p of
Int
0 ->
Point m -> m Double
i Point m
p
Int
n -> do
let sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
ty :: Double
ty = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int
0
t1 :: Double
t1 = Double
ty
t2 :: Double
t2 = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int
1
t3 :: Double
t3 = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int
2
t4 :: Double
t4 = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int
3
py :: Point m
py = Point m
p { pointTime :: Double
pointTime = Double
ty, pointIteration :: Int
pointIteration = Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1, pointPhase :: Int
pointPhase = Int
0 }
p1 :: Point m
p1 = Point m
py
p2 :: Point m
p2 = Point m
p { pointTime :: Double
pointTime = Double
t2, pointIteration :: Int
pointIteration = Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1, pointPhase :: Int
pointPhase = Int
1 }
p3 :: Point m
p3 = Point m
p { pointTime :: Double
pointTime = Double
t3, pointIteration :: Int
pointIteration = Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1, pointPhase :: Int
pointPhase = Int
2 }
p4 :: Point m
p4 = Point m
p { pointTime :: Double
pointTime = Double
t4, pointIteration :: Int
pointIteration = Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1, pointPhase :: Int
pointPhase = Int
3 }
Double
vy <- Point m -> m Double
y Point m
py
Double
k1 <- Point m -> m Double
f Point m
p1
Double
k2 <- Point m -> m Double
f Point m
p2
Double
k3 <- Point m -> m Double
f Point m
p3
Double
k4 <- Point m -> m Double
f Point m
p4
let !v :: Double
v = Double
vy Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
8.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
k1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
3.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
k2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
k3) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
k4)
Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return Double
v
Int
1 -> do
let sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
n :: Int
n  = Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p
ty :: Double
ty = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc Int
n Int
0
t1 :: Double
t1 = Double
ty
py :: Point m
py = Point m
p { pointTime :: Double
pointTime = Double
ty, pointIteration :: Int
pointIteration = Int
n, pointPhase :: Int
pointPhase = Int
0 }
p1 :: Point m
p1 = Point m
py
Double
vy <- Point m -> m Double
y Point m
py
Double
k1 <- Point m -> m Double
f Point m
p1
let !v :: Double
v = Double
vy Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
k1
Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return Double
v
Int
2 -> do
let sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
n :: Int
n  = Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p
ty :: Double
ty = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc Int
n Int
0
t1 :: Double
t1 = Double
ty
t2 :: Double
t2 = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc Int
n Int
1
py :: Point m
py = Point m
p { pointTime :: Double
pointTime = Double
ty, pointIteration :: Int
pointIteration = Int
n, pointPhase :: Int
pointPhase = Int
0 }
p1 :: Point m
p1 = Point m
py
p2 :: Point m
p2 = Point m
p { pointTime :: Double
pointTime = Double
t2, pointIteration :: Int
pointIteration = Int
n, pointPhase :: Int
pointPhase = Int
1 }
Double
vy <- Point m -> m Double
y Point m
py
Double
k1 <- Point m -> m Double
f Point m
p1
Double
k2 <- Point m -> m Double
f Point m
p2
let !v :: Double
v = Double
vy Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc Double -> Double -> Double
forall a. Num a => a -> a -> a
* (- Double
k1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
k2)
Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return Double
v
Int
3 -> do
let sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
n :: Int
n  = Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p
ty :: Double
ty = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc Int
n Int
0
t1 :: Double
t1 = Double
ty
t2 :: Double
t2 = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc Int
n Int
1
t3 :: Double
t3 = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc Int
n Int
2
py :: Point m
py = Point m
p { pointTime :: Double
pointTime = Double
ty, pointIteration :: Int
pointIteration = Int
n, pointPhase :: Int
pointPhase = Int
0 }
p1 :: Point m
p1 = Point m
py
p2 :: Point m
p2 = Point m
p { pointTime :: Double
pointTime = Double
t2, pointIteration :: Int
pointIteration = Int
n, pointPhase :: Int
pointPhase = Int
1 }
p3 :: Point m
p3 = Point m
p { pointTime :: Double
pointTime = Double
t3, pointIteration :: Int
pointIteration = Int
n, pointPhase :: Int
pointPhase = Int
2 }
Double
vy <- Point m -> m Double
y Point m
py
Double
k1 <- Point m -> m Double
f Point m
p1
Double
k2 <- Point m -> m Double
f Point m
p2
Double
k3 <- Point m -> m Double
f Point m
p3
let !v :: Double
v = Double
vy Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
k1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
k2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
k3)
Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return Double
v
Int
_ ->
[Char] -> m Double
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 =
--   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]
-- @
=> Dynamics m Double                  -- ^ the derivative
-> Dynamics m Double                  -- ^ the initial value
-> Simulation m (Dynamics m Double)   -- ^ the integral
{-# INLINABLE integ #-}
integ :: Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
integ Dynamics m Double
diff Dynamics m Double
i =
mdo Dynamics m Double
y <- Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *) e.
Dynamics m e -> Simulation m (Dynamics m e)
MU.memoDynamics Dynamics m Double
z
Dynamics m Double
z <- (Run m -> m (Dynamics m Double))
-> Simulation m (Dynamics m Double)
forall (m :: * -> *) a. (Run m -> m a) -> Simulation m a
Simulation ((Run m -> m (Dynamics m Double))
-> Simulation m (Dynamics m Double))
-> (Run m -> m (Dynamics m Double))
-> Simulation m (Dynamics m Double)
forall a b. (a -> b) -> a -> b
\$ \Run m
r ->
case Specs m -> Method
forall (m :: * -> *). Specs m -> Method
spcMethod (Run m -> Specs m
forall (m :: * -> *). Run m -> Specs m
runSpecs Run m
r) of
Method
Euler -> Dynamics m Double -> m (Dynamics m Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics m Double -> m (Dynamics m Double))
-> Dynamics m Double -> m (Dynamics m Double)
forall a b. (a -> b) -> a -> b
\$ (Point m -> m Double) -> Dynamics m Double
forall (m :: * -> *) a. (Point m -> m a) -> Dynamics m a
Dynamics ((Point m -> m Double) -> Dynamics m Double)
-> (Point m -> m Double) -> Dynamics m Double
forall a b. (a -> b) -> a -> b
\$ Dynamics m Double
-> Dynamics m Double -> Dynamics m Double -> Point m -> m Double
forall (m :: * -> *).
Dynamics m Double
-> Dynamics m Double -> Dynamics m Double -> Point m -> m Double
integEuler Dynamics m Double
diff Dynamics m Double
i Dynamics m Double
y
Method
RungeKutta2 -> Dynamics m Double -> m (Dynamics m Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics m Double -> m (Dynamics m Double))
-> Dynamics m Double -> m (Dynamics m Double)
forall a b. (a -> b) -> a -> b
\$ (Point m -> m Double) -> Dynamics m Double
forall (m :: * -> *) a. (Point m -> m a) -> Dynamics m a
Dynamics ((Point m -> m Double) -> Dynamics m Double)
-> (Point m -> m Double) -> Dynamics m Double
forall a b. (a -> b) -> a -> b
\$ Dynamics m Double
-> Dynamics m Double -> Dynamics m Double -> Point m -> m Double
forall (m :: * -> *).
Dynamics m Double
-> Dynamics m Double -> Dynamics m Double -> Point m -> m Double
integRK2 Dynamics m Double
diff Dynamics m Double
i Dynamics m Double
y
Method
RungeKutta4 -> Dynamics m Double -> m (Dynamics m Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics m Double -> m (Dynamics m Double))
-> Dynamics m Double -> m (Dynamics m Double)
forall a b. (a -> b) -> a -> b
\$ (Point m -> m Double) -> Dynamics m Double
forall (m :: * -> *) a. (Point m -> m a) -> Dynamics m a
Dynamics ((Point m -> m Double) -> Dynamics m Double)
-> (Point m -> m Double) -> Dynamics m Double
forall a b. (a -> b) -> a -> b
\$ Dynamics m Double
-> Dynamics m Double -> Dynamics m Double -> Point m -> m Double
forall (m :: * -> *).
Dynamics m Double
-> Dynamics m Double -> Dynamics m Double -> Point m -> m Double
integRK4 Dynamics m Double
diff Dynamics m Double
i Dynamics m Double
y
Method
RungeKutta4b -> Dynamics m Double -> m (Dynamics m Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics m Double -> m (Dynamics m Double))
-> Dynamics m Double -> m (Dynamics m Double)
forall a b. (a -> b) -> a -> b
\$ (Point m -> m Double) -> Dynamics m Double
forall (m :: * -> *) a. (Point m -> m a) -> Dynamics m a
Dynamics ((Point m -> m Double) -> Dynamics m Double)
-> (Point m -> m Double) -> Dynamics m Double
forall a b. (a -> b) -> a -> b
\$ Dynamics m Double
-> Dynamics m Double -> Dynamics m Double -> Point m -> m Double
forall (m :: * -> *).
Dynamics m Double
-> Dynamics m Double -> Dynamics m Double -> Point m -> m Double
integRK4b Dynamics m Double
diff Dynamics m Double
i Dynamics m Double
y
Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *) a. Monad m => a -> m a
return Dynamics m Double
y

=> Dynamics m (Either Double Double)
-> Dynamics m Double
-> Dynamics m Double
-> Point m
-> m Double
{-# INLINABLE integEulerEither #-}
integEulerEither :: Dynamics m (Either Double Double)
-> Dynamics m Double -> Dynamics m Double -> Point m -> m Double
integEulerEither (Dynamics Point m -> m (Either Double Double)
f) (Dynamics Point m -> m Double
i) (Dynamics Point m -> m Double
y) Point m
p =
case Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p of
Int
0 ->
Point m -> m Double
i Point m
p
Int
n -> do
let sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
ty :: Double
ty = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int
0
py :: Point m
py = Point m
p { pointTime :: Double
pointTime = Double
ty, pointIteration :: Int
pointIteration = Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1, pointPhase :: Int
pointPhase = Int
0 }
Either Double Double
b <- Point m -> m (Either Double Double)
f Point m
py
case Either Double Double
b of
Left Double
v ->
Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return Double
v
Right Double
b -> do
Double
a <- Point m -> m Double
y Point m
py
let !v :: Double
v = Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT (Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b
Double -> m Double
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.
=> Dynamics m (Either Double Double)
-- ^ either set a new 'Left' integral value, or use a 'Right' derivative
-> Dynamics m Double
-- ^ the initial value
-> Simulation m (Dynamics m Double)
{-# INLINABLE integEither #-}
integEither :: Dynamics m (Either Double Double)
-> Dynamics m Double -> Simulation m (Dynamics m Double)
integEither Dynamics m (Either Double Double)
diff Dynamics m Double
i =
mdo Dynamics m Double
y <- Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *) e.
Dynamics m e -> Simulation m (Dynamics m e)
MU.memoDynamics Dynamics m Double
z
Dynamics m Double
z <- (Run m -> m (Dynamics m Double))
-> Simulation m (Dynamics m Double)
forall (m :: * -> *) a. (Run m -> m a) -> Simulation m a
Simulation ((Run m -> m (Dynamics m Double))
-> Simulation m (Dynamics m Double))
-> (Run m -> m (Dynamics m Double))
-> Simulation m (Dynamics m Double)
forall a b. (a -> b) -> a -> b
\$ \Run m
r ->
Dynamics m Double -> m (Dynamics m Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics m Double -> m (Dynamics m Double))
-> Dynamics m Double -> m (Dynamics m Double)
forall a b. (a -> b) -> a -> b
\$ (Point m -> m Double) -> Dynamics m Double
forall (m :: * -> *) a. (Point m -> m a) -> Dynamics m a
Dynamics ((Point m -> m Double) -> Dynamics m Double)
-> (Point m -> m Double) -> Dynamics m Double
forall a b. (a -> b) -> a -> b
\$ Dynamics m (Either Double Double)
-> Dynamics m Double -> Dynamics m Double -> Point m -> m Double
forall (m :: * -> *).
Dynamics m (Either Double Double)
-> Dynamics m Double -> Dynamics m Double -> Point m -> m Double
integEulerEither Dynamics m (Either Double Double)
diff Dynamics m Double
i Dynamics m Double
y
Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *) a. Monad m => a -> m a
return Dynamics m 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
-- @
=> Dynamics m Double                  -- ^ the value to smooth over time
-> Dynamics m Double                  -- ^ time
-> Dynamics m Double                  -- ^ the initial value
-> Simulation m (Dynamics m Double)   -- ^ the first order exponential smooth
{-# INLINABLE smoothI #-}
smoothI :: Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
smoothI Dynamics m Double
x Dynamics m Double
t Dynamics m Double
i =
mdo Dynamics m Double
y <- Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *).
Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
integ ((Dynamics m Double
x Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
- Dynamics m Double
y) Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
t) Dynamics m Double
i
Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *) a. Monad m => a -> m a
return Dynamics m Double
y

-- | Return the first order exponential smooth.
--
-- This is a simplified version of the 'smoothI' function
-- without specifing the initial value.
=> Dynamics m Double                  -- ^ the value to smooth over time
-> Dynamics m Double                  -- ^ time
-> Simulation m (Dynamics m Double)   -- ^ the first order exponential smooth
{-# INLINABLE smooth #-}
smooth :: Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
smooth Dynamics m Double
x Dynamics m Double
t = Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
forall (m :: * -> *).
Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
smoothI Dynamics m Double
x Dynamics m Double
t Dynamics m 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
-- @
=> Dynamics m Double                  -- ^ the value to smooth over time
-> Dynamics m Double                  -- ^ time
-> Dynamics m Double                  -- ^ the initial value
-> Simulation m (Dynamics m Double)   -- ^ the third order exponential smooth
{-# INLINABLE smooth3I #-}
smooth3I :: Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
smooth3I Dynamics m Double
x Dynamics m Double
t Dynamics m Double
i =
mdo Dynamics m Double
y  <- Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *).
Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
integ ((Dynamics m Double
s2 Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
- Dynamics m Double
y) Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
t') Dynamics m Double
i
Dynamics m Double
s2 <- Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *).
Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
integ ((Dynamics m Double
s1 Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
- Dynamics m Double
s2) Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
t') Dynamics m Double
i
Dynamics m Double
s1 <- Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *).
Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
integ ((Dynamics m Double
x Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
- Dynamics m Double
s1) Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
t') Dynamics m Double
i
let t' :: Dynamics m Double
t' = Dynamics m Double
t Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
3.0
Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *) a. Monad m => a -> m a
return Dynamics m Double
y

-- | Return the third order exponential smooth.
--
-- This is a simplified version of the 'smooth3I' function
-- without specifying the initial value.
=> Dynamics m Double                  -- ^ the value to smooth over time
-> Dynamics m Double                  -- ^ time
-> Simulation m (Dynamics m Double)   -- ^ the third order exponential smooth
{-# INLINABLE smooth3 #-}
smooth3 :: Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
smooth3 Dynamics m Double
x Dynamics m Double
t = Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
forall (m :: * -> *).
Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
smooth3I Dynamics m Double
x Dynamics m Double
t Dynamics m 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.
=> Dynamics m Double                  -- ^ the value to smooth over time
-> Dynamics m Double                  -- ^ time
-> Int                                -- ^ the order
-> Dynamics m Double                  -- ^ the initial value
-> Simulation m (Dynamics m Double)   -- ^ the n'th order exponential smooth
{-# INLINABLE smoothNI #-}
smoothNI :: Dynamics m Double
-> Dynamics m Double
-> Int
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
smoothNI Dynamics m Double
x Dynamics m Double
t Int
n Dynamics m Double
i =
mdo [Dynamics m Double]
s <- [Int]
-> (Int -> Simulation m (Dynamics m Double))
-> Simulation m [Dynamics m Double]
forall (t :: * -> *) (m :: * -> *) a b.
t a -> (a -> m b) -> m (t b)
forM [Int
1 .. Int
n] ((Int -> Simulation m (Dynamics m Double))
-> Simulation m [Dynamics m Double])
-> (Int -> Simulation m (Dynamics m Double))
-> Simulation m [Dynamics m Double]
forall a b. (a -> b) -> a -> b
\$ \Int
k ->
if Int
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1
then Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *).
Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
integ ((Dynamics m Double
x Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
- Array Int (Dynamics m Double)
a Array Int (Dynamics m Double) -> Int -> Dynamics m Double
forall i e. Ix i => Array i e -> i -> e
! Int
1) Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
t') Dynamics m Double
i
else Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *).
Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
integ ((Array Int (Dynamics m Double)
a Array Int (Dynamics m Double) -> Int -> Dynamics m Double
forall i e. Ix i => Array i e -> i -> e
! (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
- Array Int (Dynamics m Double)
a Array Int (Dynamics m Double) -> Int -> Dynamics m Double
forall i e. Ix i => Array i e -> i -> e
! Int
k) Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
t') Dynamics m Double
i
let a :: Array Int (Dynamics m Double)
a  = (Int, Int) -> [Dynamics m Double] -> Array Int (Dynamics m Double)
forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (Int
1, Int
n) [Dynamics m Double]
s
t' :: Dynamics m Double
t' = Dynamics m Double
t Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Int -> Dynamics m Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics m Double -> Simulation m (Dynamics m Double))
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall a b. (a -> b) -> a -> b
\$ Array Int (Dynamics m Double)
a Array Int (Dynamics m Double) -> Int -> Dynamics m Double
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.
=> Dynamics m Double                  -- ^ the value to smooth over time
-> Dynamics m Double                  -- ^ time
-> Int                                -- ^ the order
-> Simulation m (Dynamics m Double)   -- ^ the n'th order exponential smooth
{-# INLINABLE smoothN #-}
smoothN :: Dynamics m Double
-> Dynamics m Double -> Int -> Simulation m (Dynamics m Double)
smoothN Dynamics m Double
x Dynamics m Double
t Int
n = Dynamics m Double
-> Dynamics m Double
-> Int
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
forall (m :: * -> *).
Dynamics m Double
-> Dynamics m Double
-> Int
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
smoothNI Dynamics m Double
x Dynamics m Double
t Int
n Dynamics m 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
-- @
=> Dynamics m Double                  -- ^ the value to conserve
-> Dynamics m Double                  -- ^ time
-> Dynamics m Double                  -- ^ the initial value
-> Simulation m (Dynamics m Double)   -- ^ the first order exponential delay
{-# INLINABLE delay1I #-}
delay1I :: Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
delay1I Dynamics m Double
x Dynamics m Double
t Dynamics m Double
i =
mdo Dynamics m Double
y <- Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *).
Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
integ (Dynamics m Double
x Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
- Dynamics m Double
y Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
t) (Dynamics m Double
i Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m Double
t)
Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics m Double -> Simulation m (Dynamics m Double))
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall a b. (a -> b) -> a -> b
\$ Dynamics m Double
y Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
t

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

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

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

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

-- | Return the n'th order exponential delay.
--
-- This is a simplified version of the 'delayNI' function
-- without specifying the initial value.
=> Dynamics m Double                  -- ^ the value to conserve
-> Dynamics m Double                  -- ^ time
-> Int                                -- ^ the order
-> Simulation m (Dynamics m Double)   -- ^ the n'th order exponential delay
{-# INLINABLE delayN #-}
delayN :: Dynamics m Double
-> Dynamics m Double -> Int -> Simulation m (Dynamics m Double)
delayN Dynamics m Double
x Dynamics m Double
t Int
n = Dynamics m Double
-> Dynamics m Double
-> Int
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
forall (m :: * -> *).
Dynamics m Double
-> Dynamics m Double
-> Int
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
delayNI Dynamics m Double
x Dynamics m Double
t Int
n Dynamics m 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)
-- @
=> Dynamics m Double                  -- ^ the value to forecast
-> Dynamics m Double                  -- ^ the average time
-> Dynamics m Double                  -- ^ the time horizon
-> Simulation m (Dynamics m Double)   -- ^ the forecast
{-# INLINABLE forecast #-}
forecast :: Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
forecast Dynamics m Double
x Dynamics m Double
at Dynamics m Double
hz =
do Dynamics m Double
y <- Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *).
Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
smooth Dynamics m Double
x Dynamics m Double
at
Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics m Double -> Simulation m (Dynamics m Double))
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall a b. (a -> b) -> a -> b
\$ Dynamics m Double
x Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* (Dynamics m Double
1.0 Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
+ (Dynamics m Double
x Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
y Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
- Dynamics m Double
1.0) Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
at Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m 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
-- @
=> Dynamics m Double                  -- ^ the value for which the trend is calculated
-> Dynamics m Double                  -- ^ the average time
-> Dynamics m Double                  -- ^ the initial value
-> Simulation m (Dynamics m Double)   -- ^ the fractional change rate
{-# INLINABLE trend #-}
trend :: Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
trend Dynamics m Double
x Dynamics m Double
at Dynamics m Double
i =
do Dynamics m Double
y <- Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
forall (m :: * -> *).
Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
smoothI Dynamics m Double
x Dynamics m Double
at (Dynamics m Double
x Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ (Dynamics m Double
1.0 Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
+ Dynamics m Double
i Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m Double
at))
Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics m Double -> Simulation m (Dynamics m Double))
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall a b. (a -> b) -> a -> b
\$ (Dynamics m Double
x Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m Double
y Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
- Dynamics m Double
1.0) Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ Dynamics m 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.
=> Dynamics m a                  -- ^ the difference
-> Dynamics m a                  -- ^ the initial value
-> Simulation m (Dynamics m a)   -- ^ the sum
{-# INLINABLE diffsum #-}
diffsum :: Dynamics m a -> Dynamics m a -> Simulation m (Dynamics m a)
diffsum (Dynamics Point m -> m a
diff) (Dynamics Point m -> m a
i) =
mdo Dynamics m a
y <-
Dynamics m a -> Simulation m (Dynamics m a)
forall (m :: * -> *) e.
Dynamics m e -> Simulation m (Dynamics m e)
MU.memo0Dynamics (Dynamics m a -> Simulation m (Dynamics m a))
-> Dynamics m a -> Simulation m (Dynamics m a)
forall a b. (a -> b) -> a -> b
\$
(Point m -> m a) -> Dynamics m a
forall (m :: * -> *) a. (Point m -> m a) -> Dynamics m a
Dynamics ((Point m -> m a) -> Dynamics m a)
-> (Point m -> m a) -> Dynamics m a
forall a b. (a -> b) -> a -> b
\$ \Point m
p ->
case Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p of
Int
0 -> Point m -> m a
i Point m
p
Int
n -> do
let Dynamics Point m -> m a
m = Dynamics m a
y
sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
ty :: Double
ty = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int
0
py :: Point m
py = Point m
p { pointTime :: Double
pointTime = Double
ty,
pointIteration :: Int
pointIteration = Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1,
pointPhase :: Int
pointPhase = Int
0 }
a
a <- Point m -> m a
m Point m
py
a
b <- Point m -> m a
diff Point m
py
let !v :: a
v = a
a a -> a -> a
forall a. Num a => a -> a -> a
+ a
b
a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return a
v
Dynamics m a -> Simulation m (Dynamics m a)
forall (m :: * -> *) a. Monad m => a -> m a
return Dynamics m a
y

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

--
-- Table Functions
--

-- | Lookup @x@ in a table of pairs @(x, y)@ using linear interpolation.
lookupDynamics :: Monad m => Dynamics m Double -> Array Int (Double, Double) -> Dynamics m Double
{-# INLINABLE lookupDynamics #-}
lookupDynamics :: Dynamics m Double
-> Array Int (Double, Double) -> Dynamics m Double
lookupDynamics (Dynamics Point m -> m Double
m) Array Int (Double, Double)
tbl =
(Point m -> m Double) -> Dynamics m Double
forall (m :: * -> *) a. (Point m -> m a) -> Dynamics m a
Dynamics ((Point m -> m Double) -> Dynamics m Double)
-> (Point m -> m Double) -> Dynamics m Double
forall a b. (a -> b) -> a -> b
\$ \Point m
p ->
do Double
a <- Point m -> m Double
m Point m
p
Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> m Double) -> Double -> m Double
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 :: Monad m => Dynamics m Double -> Array Int (Double, Double) -> Dynamics m Double
{-# INLINABLE lookupStepwiseDynamics #-}
lookupStepwiseDynamics :: Dynamics m Double
-> Array Int (Double, Double) -> Dynamics m Double
lookupStepwiseDynamics (Dynamics Point m -> m Double
m) Array Int (Double, Double)
tbl =
(Point m -> m Double) -> Dynamics m Double
forall (m :: * -> *) a. (Point m -> m a) -> Dynamics m a
Dynamics ((Point m -> m Double) -> Dynamics m Double)
-> (Point m -> m Double) -> Dynamics m Double
forall a b. (a -> b) -> a -> b
\$ \Point m
p ->
do Double
a <- Point m -> m Double
m Point m
p
Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> m Double) -> Double -> m Double
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'.
=> Dynamics m a          -- ^ the value to delay
-> Dynamics m Double     -- ^ the lag time
-> Dynamics m a          -- ^ the delayed value
{-# INLINABLE delay #-}
delay :: Dynamics m a -> Dynamics m Double -> Dynamics m a
delay (Dynamics Point m -> m a
x) (Dynamics Point m -> m Double
d) = Dynamics m a -> Dynamics m a
forall (m :: * -> *) a. Dynamics m a -> Dynamics m a
discreteDynamics (Dynamics m a -> Dynamics m a) -> Dynamics m a -> Dynamics m a
forall a b. (a -> b) -> a -> b
\$ (Point m -> m a) -> Dynamics m a
forall (m :: * -> *) a. (Point m -> m a) -> Dynamics m a
Dynamics Point m -> m a
r
where
r :: Point m -> m a
r Point m
p = do
let t :: Double
t  = Point m -> Double
forall (m :: * -> *). Point m -> Double
pointTime Point m
p
sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
n :: Int
n  = Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p
Double
a <- Point m -> m Double
d Point m
p
let t' :: Double
t' = Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a
n' :: Int
n' = Integer -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Integer -> Int) -> Integer -> Int
forall a b. (a -> b) -> a -> b
\$ Double -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
floor (Double -> Integer) -> Double -> Integer
forall a b. (a -> b) -> a -> b
\$ (Double
t' Double -> Double -> Double
forall a. Num a => a -> a -> a
- Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcStartTime Specs m
sc) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc
y :: m a
y | Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0    = Point m -> m a
x (Point m -> m a) -> Point m -> m a
forall a b. (a -> b) -> a -> b
\$ Point m
p { pointTime :: Double
pointTime = Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcStartTime Specs m
sc,
pointIteration :: Int
pointIteration = Int
0,
pointPhase :: Int
pointPhase = Int
0 }
| Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
n    = Point m -> m a
x (Point m -> m a) -> Point m -> m a
forall a b. (a -> b) -> a -> b
\$ Point m
p { pointTime :: Double
pointTime = Double
t',
pointIteration :: Int
pointIteration = Int
n',
pointPhase :: Int
pointPhase = -Int
1 }
| Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
n    = [Char] -> m a
forall a. HasCallStack => [Char] -> a
error ([Char] -> m a) -> [Char] -> m a
forall a b. (a -> b) -> a -> b
\$
[Char]
"Cannot return the future data: delay. " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++
[Char]
"The lag time cannot be negative."
| Bool
otherwise = [Char] -> m a
forall a. HasCallStack => [Char] -> a
error ([Char] -> m a) -> [Char] -> m a
forall a b. (a -> b) -> a -> b
\$
[Char]
"Cannot return the current data: delay. " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++
[Char]
"The lag time is too small."
m 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'.
=> Dynamics m a                    -- ^ the value to delay
-> Dynamics m Double               -- ^ the lag time
-> Dynamics m a                    -- ^ the initial value
-> Simulation m (Dynamics m a)     -- ^ the delayed value
{-# INLINABLE delayI #-}
delayI :: Dynamics m a
-> Dynamics m Double -> Dynamics m a -> Simulation m (Dynamics m a)
delayI (Dynamics Point m -> m a
x) (Dynamics Point m -> m Double
d) (Dynamics Point m -> m a
i) = Dynamics m a -> Simulation m (Dynamics m a)
forall (m :: * -> *) e.
Dynamics m e -> Simulation m (Dynamics m e)
M.memo0Dynamics (Dynamics m a -> Simulation m (Dynamics m a))
-> Dynamics m a -> Simulation m (Dynamics m a)
forall a b. (a -> b) -> a -> b
\$ (Point m -> m a) -> Dynamics m a
forall (m :: * -> *) a. (Point m -> m a) -> Dynamics m a
Dynamics Point m -> m a
r
where
r :: Point m -> m a
r Point m
p = do
let t :: Double
t  = Point m -> Double
forall (m :: * -> *). Point m -> Double
pointTime Point m
p
sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
n :: Int
n  = Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p
Double
a <- Point m -> m Double
d Point m
p
let t' :: Double
t' = Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a
n' :: Int
n' = Integer -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Integer -> Int) -> Integer -> Int
forall a b. (a -> b) -> a -> b
\$ Double -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
floor (Double -> Integer) -> Double -> Integer
forall a b. (a -> b) -> a -> b
\$ (Double
t' Double -> Double -> Double
forall a. Num a => a -> a -> a
- Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcStartTime Specs m
sc) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc
y :: m a
y | Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0    = Point m -> m a
i (Point m -> m a) -> Point m -> m a
forall a b. (a -> b) -> a -> b
\$ Point m
p { pointTime :: Double
pointTime = Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcStartTime Specs m
sc,
pointIteration :: Int
pointIteration = Int
0,
pointPhase :: Int
pointPhase = Int
0 }
| Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
n    = Point m -> m a
x (Point m -> m a) -> Point m -> m a
forall a b. (a -> b) -> a -> b
\$ Point m
p { pointTime :: Double
pointTime = Double
t',
pointIteration :: Int
pointIteration = Int
n',
pointPhase :: Int
pointPhase = -Int
1 }
| Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
n    = [Char] -> m a
forall a. HasCallStack => [Char] -> a
error ([Char] -> m a) -> [Char] -> m a
forall a b. (a -> b) -> a -> b
\$
[Char]
"Cannot return the future data: delay. " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++
[Char]
"The lag time cannot be negative."
| Bool
otherwise = [Char] -> m a
forall a. HasCallStack => [Char] -> a
error ([Char] -> m a) -> [Char] -> m a
forall a b. (a -> b) -> a -> b
\$
[Char]
"Cannot return the current data: delay. " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++
[Char]
"The lag time is too small."
m a
y

-- | Return the delayed value by the specified positive number of
-- integration time steps used for calculating the lag time.
=> Dynamics m a
-- ^ the value to delay
-> Dynamics m Int
-- ^ the delay as a multiplication of the corresponding number
-- and the integration time step
-> Dynamics m a
-- ^ the delayed value
{-# INLINABLE delayByDT #-}
delayByDT :: Dynamics m a -> Dynamics m Int -> Dynamics m a
delayByDT (Dynamics Point m -> m a
x) (Dynamics Point m -> m Int
d) = Dynamics m a -> Dynamics m a
forall (m :: * -> *) a. Dynamics m a -> Dynamics m a
discreteDynamics (Dynamics m a -> Dynamics m a) -> Dynamics m a -> Dynamics m a
forall a b. (a -> b) -> a -> b
\$ (Point m -> m a) -> Dynamics m a
forall (m :: * -> *) a. (Point m -> m a) -> Dynamics m a
Dynamics Point m -> m a
r
where
r :: Point m -> m a
r Point m
p = do
let sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
n :: Int
n  = Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p
Int
a <- Point m -> m Int
d Point m
p
let p' :: Point m
p' = Point m -> Int -> Point m
forall (m :: * -> *). Point m -> Int -> Point m
delayPoint Point m
p Int
a
n' :: Int
n' = Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p'
y :: m a
y | Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0    = Point m -> m a
x (Point m -> m a) -> Point m -> m a
forall a b. (a -> b) -> a -> b
\$ Point m
p { pointTime :: Double
pointTime = Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcStartTime Specs m
sc,
pointIteration :: Int
pointIteration = Int
0,
pointPhase :: Int
pointPhase = Int
0 }
| Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
n    = Point m -> m a
x Point m
p'
| Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
n    = [Char] -> m a
forall a. HasCallStack => [Char] -> a
error ([Char] -> m a) -> [Char] -> m a
forall a b. (a -> b) -> a -> b
\$
[Char]
"Cannot return the future data: delayByDT. " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++
[Char]
"The lag time cannot be negative."
| Bool
otherwise = [Char] -> m a
forall a. HasCallStack => [Char] -> a
error ([Char] -> m a) -> [Char] -> m a
forall a b. (a -> b) -> a -> b
\$
[Char]
"Cannot return the current data: delayByDT. " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++
[Char]
"The lag time is too small."
m 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.
=> Dynamics m a
-- ^ the value to delay
-> Dynamics m Int
-- ^ the delay as a multiplication of the corresponding number
-- and the integration time step
-> Dynamics m a
-- ^ the initial value
-> Simulation m (Dynamics m a)
-- ^ the delayed value
{-# INLINABLE delayIByDT #-}
delayIByDT :: Dynamics m a
-> Dynamics m Int -> Dynamics m a -> Simulation m (Dynamics m a)
delayIByDT (Dynamics Point m -> m a
x) (Dynamics Point m -> m Int
d) (Dynamics Point m -> m a
i) = Dynamics m a -> Simulation m (Dynamics m a)
forall (m :: * -> *) e.
Dynamics m e -> Simulation m (Dynamics m e)
M.memoDynamics (Dynamics m a -> Simulation m (Dynamics m a))
-> Dynamics m a -> Simulation m (Dynamics m a)
forall a b. (a -> b) -> a -> b
\$ (Point m -> m a) -> Dynamics m a
forall (m :: * -> *) a. (Point m -> m a) -> Dynamics m a
Dynamics Point m -> m a
r
where
r :: Point m -> m a
r Point m
p = do
let sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
n :: Int
n  = Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p
Int
a <- Point m -> m Int
d Point m
p
let p' :: Point m
p' = Point m -> Int -> Point m
forall (m :: * -> *). Point m -> Int -> Point m
delayPoint Point m
p Int
a
n' :: Int
n' = Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p'
y :: m a
y | Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0    = Point m -> m a
i (Point m -> m a) -> Point m -> m a
forall a b. (a -> b) -> a -> b
\$ Point m
p { pointTime :: Double
pointTime = Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcStartTime Specs m
sc,
pointIteration :: Int
pointIteration = Int
0,
pointPhase :: Int
pointPhase = Int
0 }
| Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
n    = Point m -> m a
x Point m
p'
| Int
n' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
n    = [Char] -> m a
forall a. HasCallStack => [Char] -> a
error ([Char] -> m a) -> [Char] -> m a
forall a b. (a -> b) -> a -> b
\$
[Char]
"Cannot return the future data: delayIByDT. " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++
[Char]
"The lag time cannot be negative."
| Bool
otherwise = [Char] -> m a
forall a. HasCallStack => [Char] -> a
error ([Char] -> m a) -> [Char] -> m a
forall a b. (a -> b) -> a -> b
\$
[Char]
"Cannot return the current data: delayIByDT. " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++
[Char]
"The lag time is too small."
m 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
-- @
=> Dynamics m Double                  -- ^ the stream
-> Dynamics m Double                  -- ^ the discount rate
-> Dynamics m Double                  -- ^ the initial value
-> Dynamics m Double                  -- ^ factor
-> Simulation m (Dynamics m Double)   -- ^ the Net Present Value (NPV)
{-# INLINABLE npv #-}
npv :: Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
npv Dynamics m Double
stream Dynamics m Double
rate Dynamics m Double
init Dynamics m Double
factor =
mdo let dt' :: Dynamics m Double
dt' = Parameter m Double -> Dynamics m Double
forall (t :: (* -> *) -> * -> *) (m :: * -> *) a.
ParameterLift t m =>
Parameter m a -> t m a
liftParameter Parameter m Double
forall (m :: * -> *). Monad m => Parameter m Double
dt
Dynamics m Double
df <- Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *).
Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
integ (- Dynamics m Double
df Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m Double
rate) Dynamics m Double
1
Dynamics m Double
accum <- Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *).
Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
integ (Dynamics m Double
stream Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m Double
df) Dynamics m Double
init
Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics m Double -> Simulation m (Dynamics m Double))
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall a b. (a -> b) -> a -> b
\$ (Dynamics m Double
accum Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
+ Dynamics m Double
dt' Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m Double
stream Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m Double
df) Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m 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
-- @
=> Dynamics m Double                  -- ^ the stream
-> Dynamics m Double                  -- ^ the discount rate
-> Dynamics m Double                  -- ^ the initial value
-> Dynamics m Double                  -- ^ factor
-> Simulation m (Dynamics m Double)   -- ^ the Net Present Value End (NPVE)
{-# INLINABLE npve #-}
npve :: Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Dynamics m Double
-> Simulation m (Dynamics m Double)
npve Dynamics m Double
stream Dynamics m Double
rate Dynamics m Double
init Dynamics m Double
factor =
mdo let dt' :: Dynamics m Double
dt' = Parameter m Double -> Dynamics m Double
forall (t :: (* -> *) -> * -> *) (m :: * -> *) a.
ParameterLift t m =>
Parameter m a -> t m a
liftParameter Parameter m Double
forall (m :: * -> *). Monad m => Parameter m Double
dt
Dynamics m Double
df <- Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *).
Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
integ (- Dynamics m Double
df Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m Double
rate Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ (Dynamics m Double
1 Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
+ Dynamics m Double
rate Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m Double
dt')) (Dynamics m Double
1 Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Fractional a => a -> a -> a
/ (Dynamics m Double
1 Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
+ Dynamics m Double
rate Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m Double
dt'))
Dynamics m Double
accum <- Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *).
Dynamics m Double
-> Dynamics m Double -> Simulation m (Dynamics m Double)
integ (Dynamics m Double
stream Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m Double
df) Dynamics m Double
init
Dynamics m Double -> Simulation m (Dynamics m Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Dynamics m Double -> Simulation m (Dynamics m Double))
-> Dynamics m Double -> Simulation m (Dynamics m Double)
forall a b. (a -> b) -> a -> b
\$ (Dynamics m Double
accum Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
+ Dynamics m Double
dt' Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m Double
stream Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m Double
df) Dynamics m Double -> Dynamics m Double -> Dynamics m Double
forall a. Num a => a -> a -> a
* Dynamics m Double
factor

-- | Computation that returns 0 until the step time and then returns the specified height.
=> Dynamics m Double
-- ^ the height
-> Dynamics m Double
-- ^ the step time
-> Dynamics m Double
{-# INLINABLE step #-}
step :: Dynamics m Double -> Dynamics m Double -> Dynamics m Double
step Dynamics m Double
h Dynamics m Double
st =
Dynamics m Double -> Dynamics m Double
forall (m :: * -> *) a. Dynamics m a -> Dynamics m a
discreteDynamics (Dynamics m Double -> Dynamics m Double)
-> Dynamics m Double -> Dynamics m Double
forall a b. (a -> b) -> a -> b
\$
(Point m -> m Double) -> Dynamics m Double
forall (m :: * -> *) a. (Point m -> m a) -> Dynamics m a
Dynamics ((Point m -> m Double) -> Dynamics m Double)
-> (Point m -> m Double) -> Dynamics m Double
forall a b. (a -> b) -> a -> b
\$ \Point m
p ->
do let sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
t :: Double
t  = Point m -> Double
forall (m :: * -> *). Point m -> Double
pointTime Point m
p
Double
st' <- Point m -> Dynamics m Double -> m Double
forall (m :: * -> *) a. Point m -> Dynamics m a -> m a
invokeDynamics Point m
p Dynamics m Double
st
let t' :: Double
t' = Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
if Double
st' Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
t'
then Point m -> Dynamics m Double -> m Double
forall (m :: * -> *) a. Point m -> Dynamics m a -> m a
invokeDynamics Point m
p Dynamics m Double
h
else Double -> m Double
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.
=> Dynamics m Double
-- ^ the time start
-> Dynamics m Double
-- ^ the interval width
-> Dynamics m Double
{-# INLINABLE pulse #-}
pulse :: Dynamics m Double -> Dynamics m Double -> Dynamics m Double
pulse Dynamics m Double
st Dynamics m Double
w =
Dynamics m Double -> Dynamics m Double
forall (m :: * -> *) a. Dynamics m a -> Dynamics m a
discreteDynamics (Dynamics m Double -> Dynamics m Double)
-> Dynamics m Double -> Dynamics m Double
forall a b. (a -> b) -> a -> b
\$
(Point m -> m Double) -> Dynamics m Double
forall (m :: * -> *) a. (Point m -> m a) -> Dynamics m a
Dynamics ((Point m -> m Double) -> Dynamics m Double)
-> (Point m -> m Double) -> Dynamics m Double
forall a b. (a -> b) -> a -> b
\$ \Point m
p ->
do let sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
t :: Double
t  = Point m -> Double
forall (m :: * -> *). Point m -> Double
pointTime Point m
p
Double
st' <- Point m -> Dynamics m Double -> m Double
forall (m :: * -> *) a. Point m -> Dynamics m a -> m a
invokeDynamics Point m
p Dynamics m Double
st
let t' :: Double
t' = Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
if Double
st' Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
t'
then do Double
w' <- Point m -> Dynamics m Double -> m Double
forall (m :: * -> *) a. Point m -> Dynamics m a -> m a
invokeDynamics Point m
p Dynamics m Double
w
Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> m Double) -> Double -> m Double
forall a b. (a -> b) -> a -> b
\$ if Double
t' Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
st' Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
w' then Double
1 else Double
0
else Double -> m Double
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.
=> Dynamics m Double
-- ^ the time start
-> Dynamics m Double
-- ^ the interval width
-> Dynamics m Double
-- ^ the time period
-> Dynamics m Double
{-# INLINABLE pulseP #-}
pulseP :: Dynamics m Double
-> Dynamics m Double -> Dynamics m Double -> Dynamics m Double
pulseP Dynamics m Double
st Dynamics m Double
w Dynamics m Double
period =
Dynamics m Double -> Dynamics m Double
forall (m :: * -> *) a. Dynamics m a -> Dynamics m a
discreteDynamics (Dynamics m Double -> Dynamics m Double)
-> Dynamics m Double -> Dynamics m Double
forall a b. (a -> b) -> a -> b
\$
(Point m -> m Double) -> Dynamics m Double
forall (m :: * -> *) a. (Point m -> m a) -> Dynamics m a
Dynamics ((Point m -> m Double) -> Dynamics m Double)
-> (Point m -> m Double) -> Dynamics m Double
forall a b. (a -> b) -> a -> b
\$ \Point m
p ->
do let sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
t :: Double
t  = Point m -> Double
forall (m :: * -> *). Point m -> Double
pointTime Point m
p
Double
p'  <- Point m -> Dynamics m Double -> m Double
forall (m :: * -> *) a. Point m -> Dynamics m a -> m a
invokeDynamics Point m
p Dynamics m Double
period
Double
st' <- Point m -> Dynamics m Double -> m Double
forall (m :: * -> *) a. Point m -> Dynamics m a -> m a
invokeDynamics Point m
p Dynamics m Double
st
let y' :: Double
y' = if (Double
p' Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) Bool -> Bool -> Bool
&& (Double
t Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
st')
then Integer -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Double -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
floor (Double -> Integer) -> Double -> Integer
forall a b. (a -> b) -> a -> b
\$ (Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
st') Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
p') Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
p'
else Double
0
let st' :: Double
st' = Double
st' Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
y'
let t' :: Double
t' = Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
if Double
st' Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
t'
then do Double
w' <- Point m -> Dynamics m Double -> m Double
forall (m :: * -> *) a. Point m -> Dynamics m a -> m a
invokeDynamics Point m
p Dynamics m Double
w
Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> m Double) -> Double -> m Double
forall a b. (a -> b) -> a -> b
\$ if Double
t' Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
st' Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
w' then Double
1 else Double
0
else Double -> m Double
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.
=> Dynamics m Double
-- ^ the slope parameter
-> Dynamics m Double
-- ^ the time start
-> Dynamics m Double
-- ^ the end time
-> Dynamics m Double
{-# INLINABLE ramp #-}
ramp :: Dynamics m Double
-> Dynamics m Double -> Dynamics m Double -> Dynamics m Double
ramp Dynamics m Double
slope Dynamics m Double
st Dynamics m Double
e =
Dynamics m Double -> Dynamics m Double
forall (m :: * -> *) a. Dynamics m a -> Dynamics m a
discreteDynamics (Dynamics m Double -> Dynamics m Double)
-> Dynamics m Double -> Dynamics m Double
forall a b. (a -> b) -> a -> b
\$
(Point m -> m Double) -> Dynamics m Double
forall (m :: * -> *) a. (Point m -> m a) -> Dynamics m a
Dynamics ((Point m -> m Double) -> Dynamics m Double)
-> (Point m -> m Double) -> Dynamics m Double
forall a b. (a -> b) -> a -> b
\$ \Point m
p ->
do let sc :: Specs m
sc = Point m -> Specs m
forall (m :: * -> *). Point m -> Specs m
pointSpecs Point m
p
t :: Double
t  = Point m -> Double
forall (m :: * -> *). Point m -> Double
pointTime Point m
p
Double
st' <- Point m -> Dynamics m Double -> m Double
forall (m :: * -> *) a. Point m -> Dynamics m a -> m a
invokeDynamics Point m
p Dynamics m Double
st
if Double
st' Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
t
then do Double
slope' <- Point m -> Dynamics m Double -> m Double
forall (m :: * -> *) a. Point m -> Dynamics m a -> m a
invokeDynamics Point m
p Dynamics m Double
slope
Double
e' <- Point m -> Dynamics m Double -> m Double
forall (m :: * -> *) a. Point m -> Dynamics m a -> m a
invokeDynamics Point m
p Dynamics m Double
e
if Double
t Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
e'
then Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> m Double) -> Double -> m Double
forall a b. (a -> b) -> a -> b
\$ Double
slope' Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
st')
else Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> m Double) -> Double -> m Double
forall a b. (a -> b) -> a -> b
\$ Double
slope' Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
e' Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
st')
else Double -> m Double
forall (m :: * -> *) a. Monad m => a -> m a
return Double
0
```