Chimera model with two different delays

Chimera.hs

-- Chimera and dissipative solitons demonstrations in a system
-- with two strongly different delays \tau_2 = 100*\tau_1 from arXiv:1712.03283
--
-- How to execute:
-- 1. Save the code to Chimera.hs
-- 2. Install stack (command line interface is marked by $):
--   $ wget -qO- https://get.haskellstack.org/ | sh
-- OR
--   $ curl -sSL https://get.haskellstack.org/ | sh
-- 3. Run:
--   $ stack --resolver lts-10.6 --install-ghc runghc --package free-vector-spaces-0.1.4.0 --package dde-0.3.0 Chimera.hs
--
-- Troubleshooting:
-- 1. Make sure gcc compiler is installed by typing
--   $ gcc
-- To install gcc on Mac OS X:
--   $ xcode-select --intall
-- To install gcc on Ubuntu:
--   $ sudo apt-get install build-essential

-- Parameters to experiment with
--
-- 1. Dissipative solitons
--
--   beta = 1.6
--   gamma = 0.83
--   phi0 = -0.67
--   m = 50
--   epsilon = 0.01
--   delta = 0.009
--
-- NB: Transient time has to be at least
-- 20 long delays \tau_2,
-- i.e. 20 * 100 * 2000 = 4,000,000 iterations.
--
-- 2. Incoherent core
--
--   beta = 1.1
--   gamma = 0.5
--   phi0 = -0.5
--   m = 50
--   epsilon = 0.01
--   delta = 0.009
--
-- NB: Transient time has to be at least
-- 50 long delays \tau_2.
--
-- 3. Coherent core
--
--   beta = 1.1
--   gamma = 0.5
--   phi0 = -0.05
--   m = 50
--   epsilon = 0.01
--   delta = 0.009
--
-- NB: Transient time has to be at least
-- 50 long delays \tau_2.
--
-- Note: to compute for longer transient times, consider discarding
-- intermediate data.

import           Linear ( V2 (..) )  -- From `linear` package
import qualified Data.Vector.Storable as V  -- From `vector` package
import qualified Numeric.DDE as DDE  -- From `dde-0.3.0` package

-- Model equation
rhs phi0 = DDE.RHS derivative
  where
    derivative ((V2 x y), (DDE.Hist histSnapshots), _) = V2 x' y'
      where
        -- DDE Eq. (4) from arXiv:1712.03283
        x' = (-x - delta * y + (1 - gamma) * f x_tau1 + gamma * f x_tau2) / epsilon
        y' = x

        f = airy phi0

        -- Delay terms where tau2 / tau1 = 100
        (V2 x_tau1 _):(V2 x_tau2 _):_ = histSnapshots

        -- Constants
        epsilon = 0.01
        gamma = 0.5
        delta = 0.009

-- Airy function definition
airy phi0 x = beta / (1 + m * (sin (x + phi0))^2)
  where
    m = 50
    beta = 1.6

model phi0 hStep len1 len2 totalIter = (state1, V.map (\(V2 x y) -> x) trace)
  where
    -- Initial conditions:
    -- dynamical state and delay history.
    state0 = V2 0.0 0.0
    hist0 = V.fromList $ map (\n -> let x = sin(2 * pi * fromIntegral n / 1000)
                                    in V2 x 0.0) [1..len2]

    -- Input is ignored in `rhs`
    inp = DDE.Input $ V.replicate (totalIter + 1) 0

    delaysInSamples = [len1, len2]

    (state1, trace) = DDE.integ DDE.rk4 state0 hist0 delaysInSamples hStep (rhs phi0) inp

-- Control parameter
phi0 = -0.45

main = do
  let hStep = 0.0005  -- Integration step
      tauD1 = 1.0  -- Short delay time
      len1 = round $ tauD1 / hStep  -- Samples per short delay
      len2 = 100 * len1  -- Long delay (in samples)
      delays = 20  -- 20 long delays
      total = delays * len2

  let (state1, trace) = model phi0 hStep len1 len2 total

  mapM_ print $ V.toList trace