{-# LANGUAGE TypeFamilies #-}

-- |
-- Module     : Simulation.Aivika.Trans.Internal.Specs
-- 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
--
-- It defines the simulation specs and related stuff.
module Simulation.Aivika.Trans.Internal.Specs
       (Specs(..),
        Method(..),
        Run(..),
        Point(..),
        basicTime,
        integIterationBnds,
        integIterationHiBnd,
        integIterationLoBnd,
        integPhaseBnds,
        integPhaseHiBnd,
        integPhaseLoBnd,
        integTimes,
        integPoints,
        integPointsStartingFrom,
        integStartPoint,
        integStopPoint,
        simulationStopPoint,
        timeGrid,
        pointAt,
        delayPoint) where

import Simulation.Aivika.Trans.Internal.Types

-- | Returns the integration iterations starting from zero.
integIterations :: Specs m -> [Int]
integIterations :: Specs m -> [Int]
integIterations Specs m
sc = [Int
i1 .. Int
i2] where
  i1 :: Int
i1 = Specs m -> Int
forall (m :: * -> *). Specs m -> Int
integIterationLoBnd Specs m
sc
  i2 :: Int
i2 = Specs m -> Int
forall (m :: * -> *). Specs m -> Int
integIterationHiBnd Specs m
sc

-- | Returns the first and last integration iterations.
integIterationBnds :: Specs m -> (Int, Int)
integIterationBnds :: Specs m -> (Int, Int)
integIterationBnds Specs m
sc = (Int
i1, Int
i2) where
  i1 :: Int
i1 = Specs m -> Int
forall (m :: * -> *). Specs m -> Int
integIterationLoBnd Specs m
sc
  i2 :: Int
i2 = Specs m -> Int
forall (m :: * -> *). Specs m -> Int
integIterationHiBnd Specs m
sc

-- | Returns the first integration iteration, i.e. zero.
integIterationLoBnd :: Specs m -> Int
integIterationLoBnd :: Specs m -> Int
integIterationLoBnd Specs m
sc = Int
0

-- | Returns the last integration iteration.
integIterationHiBnd :: Specs m -> Int
integIterationHiBnd :: Specs m -> Int
integIterationHiBnd Specs m
sc =
  let n :: Int
n = Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
round ((Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcStopTime Specs m
sc 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)
  in if Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0
     then
       [Char] -> Int
forall a. HasCallStack => [Char] -> a
error ([Char] -> Int) -> [Char] -> Int
forall a b. (a -> b) -> a -> b
$
       [Char]
"Either the simulation specs are incorrect, " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++
       [Char]
"or a step time is too small, because of which " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++
       [Char]
"a floating point overflow occurred on 32-bit Haskell implementation."
     else Int
n

-- | Returns the phases for the specified simulation specs starting from zero.
integPhases :: Specs m -> [Int]
integPhases :: Specs m -> [Int]
integPhases Specs m
sc = 
  case Specs m -> Method
forall (m :: * -> *). Specs m -> Method
spcMethod Specs m
sc of
    Method
Euler -> [Int
0]
    Method
RungeKutta2 -> [Int
0, Int
1]
    Method
RungeKutta4 -> [Int
0, Int
1, Int
2, Int
3]
    Method
RungeKutta4b -> [Int
0, Int
1, Int
2, Int
3]

-- | Returns the first and last integration phases.
integPhaseBnds :: Specs m -> (Int, Int)
integPhaseBnds :: Specs m -> (Int, Int)
integPhaseBnds Specs m
sc = 
  case Specs m -> Method
forall (m :: * -> *). Specs m -> Method
spcMethod Specs m
sc of
    Method
Euler -> (Int
0, Int
0)
    Method
RungeKutta2 -> (Int
0, Int
1)
    Method
RungeKutta4 -> (Int
0, Int
3)
    Method
RungeKutta4b -> (Int
0, Int
3)

-- | Returns the first integration phase, i.e. zero.
integPhaseLoBnd :: Specs m -> Int
integPhaseLoBnd :: Specs m -> Int
integPhaseLoBnd Specs m
sc = Int
0
                  
-- | Returns the last integration phase, 0 for Euler's method, 1 for RK2 and 3 for RK4.
integPhaseHiBnd :: Specs m -> Int
integPhaseHiBnd :: Specs m -> Int
integPhaseHiBnd Specs m
sc = 
  case Specs m -> Method
forall (m :: * -> *). Specs m -> Method
spcMethod Specs m
sc of
    Method
Euler -> Int
0
    Method
RungeKutta2 -> Int
1
    Method
RungeKutta4 -> Int
3
    Method
RungeKutta4b -> Int
3

-- | Returns a simulation time for the integration point specified by 
-- the specs, iteration and phase.
basicTime :: Specs m -> Int -> Int -> Double
{-# INLINE basicTime #-}
basicTime :: Specs m -> Int -> Int -> Double
basicTime Specs m
sc Int
n Int
ph =
  if Int
ph Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 then 
    [Char] -> Double
forall a. HasCallStack => [Char] -> a
error [Char]
"Incorrect phase: basicTime"
  else
    Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcStartTime Specs m
sc Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
n' 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
+ Method -> Int -> Double
delta (Specs m -> Method
forall (m :: * -> *). Specs m -> Method
spcMethod Specs m
sc) Int
ph 
      where n' :: Double
n' = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
            delta :: Method -> Int -> Double
delta Method
Euler       Int
0 = Double
0
            delta Method
RungeKutta2 Int
0 = Double
0
            delta Method
RungeKutta2 Int
1 = Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc
            delta Method
RungeKutta4 Int
0 = Double
0
            delta Method
RungeKutta4 Int
1 = Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
            delta Method
RungeKutta4 Int
2 = Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
            delta Method
RungeKutta4 Int
3 = Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc
            delta Method
RungeKutta4b Int
0 = Double
0
            delta Method
RungeKutta4b Int
1 = Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3
            delta Method
RungeKutta4b Int
2 = Double
2 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
            delta Method
RungeKutta4b Int
3 = Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc

-- | Return the integration time values.
integTimes :: Specs m -> [Double]
integTimes :: Specs m -> [Double]
integTimes Specs m
sc = (Int -> Double) -> [Int] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Int -> Double
t [Int
nl .. Int
nu]
  where (Int
nl, Int
nu) = Specs m -> (Int, Int)
forall (m :: * -> *). Specs m -> (Int, Int)
integIterationBnds Specs m
sc
        t :: Int -> Double
t Int
n = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc Int
n Int
0

-- | Return the integration time points.
integPoints :: Run m -> [Point m]
integPoints :: Run m -> [Point m]
integPoints Run m
r = [Point m]
points
  where sc :: Specs m
sc = Run m -> Specs m
forall (m :: * -> *). Run m -> Specs m
runSpecs Run m
r
        (Int
nl, Int
nu) = Specs m -> (Int, Int)
forall (m :: * -> *). Specs m -> (Int, Int)
integIterationBnds Specs m
sc
        points :: [Point m]
points   = (Int -> Point m) -> [Int] -> [Point m]
forall a b. (a -> b) -> [a] -> [b]
map Int -> Point m
point [Int
nl .. Int
nu]
        point :: Int -> Point m
point Int
n  = Point :: forall (m :: * -> *).
Specs m -> Run m -> Double -> Int -> Int -> Point m
Point { pointSpecs :: Specs m
pointSpecs = Specs m
sc,
                           pointRun :: Run m
pointRun = Run m
r,
                           pointTime :: Double
pointTime = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc Int
n Int
0,
                           pointIteration :: Int
pointIteration = Int
n,
                           pointPhase :: Int
pointPhase = Int
0 }

-- | Return the start time point.
integStartPoint :: Run m -> Point m
integStartPoint :: Run m -> Point m
integStartPoint Run m
r = Int -> Point m
point Int
nl
  where sc :: Specs m
sc = Run m -> Specs m
forall (m :: * -> *). Run m -> Specs m
runSpecs Run m
r
        (Int
nl, Int
nu) = Specs m -> (Int, Int)
forall (m :: * -> *). Specs m -> (Int, Int)
integIterationBnds Specs m
sc
        point :: Int -> Point m
point Int
n  = Point :: forall (m :: * -> *).
Specs m -> Run m -> Double -> Int -> Int -> Point m
Point { pointSpecs :: Specs m
pointSpecs = Specs m
sc,
                           pointRun :: Run m
pointRun = Run m
r,
                           pointTime :: Double
pointTime = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc Int
n Int
0,
                           pointIteration :: Int
pointIteration = Int
n,
                           pointPhase :: Int
pointPhase = Int
0 }

-- | Return the integration stop time point.
integStopPoint :: Run m -> Point m
integStopPoint :: Run m -> Point m
integStopPoint Run m
r = Int -> Point m
point Int
nu
  where sc :: Specs m
sc = Run m -> Specs m
forall (m :: * -> *). Run m -> Specs m
runSpecs Run m
r
        (Int
nl, Int
nu) = Specs m -> (Int, Int)
forall (m :: * -> *). Specs m -> (Int, Int)
integIterationBnds Specs m
sc
        point :: Int -> Point m
point Int
n  = Point :: forall (m :: * -> *).
Specs m -> Run m -> Double -> Int -> Int -> Point m
Point { pointSpecs :: Specs m
pointSpecs = Specs m
sc,
                           pointRun :: Run m
pointRun = Run m
r,
                           pointTime :: Double
pointTime = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc Int
n Int
0,
                           pointIteration :: Int
pointIteration = Int
n,
                           pointPhase :: Int
pointPhase = Int
0 }

-- | Return the simulation stop time point.
simulationStopPoint :: Run m -> Point m
simulationStopPoint :: Run m -> Point m
simulationStopPoint Run m
r = Run m -> Double -> Point m
forall (m :: * -> *). Run m -> Double -> Point m
pointAt Run m
r (Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcStopTime (Specs m -> Double) -> Specs m -> Double
forall a b. (a -> b) -> a -> b
$ Run m -> Specs m
forall (m :: * -> *). Run m -> Specs m
runSpecs Run m
r)

-- | Return the point at the specified time.
pointAt :: Run m -> Double -> Point m
{-# INLINABLE pointAt #-}
pointAt :: Run m -> Double -> Point m
pointAt Run m
r Double
t = Point m
p
  where sc :: Specs m
sc = Run m -> Specs m
forall (m :: * -> *). Run m -> Specs m
runSpecs Run m
r
        t0 :: Double
t0 = Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcStartTime Specs m
sc
        dt :: Double
dt = Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc
        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
t Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
t0) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
dt)
        p :: Point m
p = Point :: forall (m :: * -> *).
Specs m -> Run m -> Double -> Int -> Int -> Point m
Point { pointSpecs :: Specs m
pointSpecs = Specs m
sc,
                    pointRun :: Run m
pointRun = Run m
r,
                    pointTime :: Double
pointTime = Double
t,
                    pointIteration :: Int
pointIteration = Int
n,
                    pointPhase :: Int
pointPhase = -Int
1 }

-- | Return the integration time points startin from the specified iteration.
integPointsStartingFrom :: Point m -> [Point m]
integPointsStartingFrom :: Point m -> [Point m]
integPointsStartingFrom Point m
p = [Point m]
points
  where r :: Run m
r  = Point m -> Run m
forall (m :: * -> *). Point m -> Run m
pointRun Point m
p
        sc :: Specs m
sc = Run m -> Specs m
forall (m :: * -> *). Run m -> Specs m
runSpecs Run m
r
        (Int
nl, Int
nu) = Specs m -> (Int, Int)
forall (m :: * -> *). Specs m -> (Int, Int)
integIterationBnds Specs m
sc
        n0 :: Int
n0       = if Point m -> Int
forall (m :: * -> *). Point m -> Int
pointPhase Point m
p Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0
                   then Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p
                   else Point m -> Int
forall (m :: * -> *). Point m -> Int
pointIteration Point m
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1
        points :: [Point m]
points   = (Int -> Point m) -> [Int] -> [Point m]
forall a b. (a -> b) -> [a] -> [b]
map Int -> Point m
point [Int
n0 .. Int
nu]
        point :: Int -> Point m
point Int
n  = Point :: forall (m :: * -> *).
Specs m -> Run m -> Double -> Int -> Int -> Point m
Point { pointSpecs :: Specs m
pointSpecs = Specs m
sc,
                           pointRun :: Run m
pointRun = Run m
r,
                           pointTime :: Double
pointTime = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc Int
n Int
0,
                           pointIteration :: Int
pointIteration = Int
n,
                           pointPhase :: Int
pointPhase = Int
0 }

-- | Return the indexed time values in the grid by specified size.
timeGrid :: Specs m -> Int -> [(Int, Double)]
timeGrid :: Specs m -> Int -> [(Int, Double)]
timeGrid Specs m
sc Int
n =
  let t0 :: Double
t0 = Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcStartTime Specs m
sc
      t2 :: Double
t2 = Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcStopTime Specs m
sc
      n' :: Int
n' = Int -> Int -> Int
forall a. Ord a => a -> a -> a
max (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int
1
      dt :: Double
dt = (Double
t2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
t0) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n')
      f :: Int -> (Int, Double)
f Int
i | Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0    = (Int
i, Double
t0)
          | Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
n'   = (Int
i, Double
t2)
          | Bool
otherwise = (Int
i, Double
t0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
i) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
dt)
  in (Int -> (Int, Double)) -> [Int] -> [(Int, Double)]
forall a b. (a -> b) -> [a] -> [b]
map Int -> (Int, Double)
f [Int
0 .. Int
n']

-- | Delay the point by the specified positive number of iterations.
delayPoint :: Point m -> Int -> Point m
delayPoint :: Point m -> Int -> Point m
delayPoint Point m
p Int
dn
  | Int
dn Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
0   = [Char] -> Point m
forall a. HasCallStack => [Char] -> a
error [Char]
"Expected the positive number of iterations: delayPoint"
  | Bool
otherwise =
    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
        ph :: Int
ph = Point m -> Int
forall (m :: * -> *). Point m -> Int
pointPhase Point m
p
    in if Int
ph Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0
       then let t' :: Double
t' = Point m -> Double
forall (m :: * -> *). Point m -> Double
pointTime Point m
p Double -> Double -> Double
forall a. Num a => a -> a -> a
- Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
dn Double -> Double -> Double
forall a. Num a => a -> a -> a
* Specs m -> Double
forall (m :: * -> *). Specs m -> Double
spcDT Specs m
sc
                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
            in Point m
p { pointTime :: Double
pointTime = Double
t',
                   pointIteration :: Int
pointIteration = Int
n',
                   pointPhase :: Int
pointPhase = -Int
1 }
       else let n' :: Int
n' = Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
dn
                t' :: Double
t' = Specs m -> Int -> Int -> Double
forall (m :: * -> *). Specs m -> Int -> Int -> Double
basicTime Specs m
sc Int
n' Int
ph
            in Point m
p { pointTime :: Double
pointTime = Double
t',
                   pointIteration :: Int
pointIteration = Int
n',
                   pointPhase :: Int
pointPhase = Int
ph }