module OsciDiffEq where

{-
ghci -fglasgow-exts -fno-implicit-prelude ../numericprelude/MyPrelude.hs ../numericprelude/NumExtras.lhs ../numericprelude/VectorSpace.lhs ../numericprelude/PreludeBase.lhs ../numericprelude/NumericPrelude.lhs OsciDiffEq.hs

import MyPrelude
but then (+) ends up in a loop Exception :-(
-}

import Number.Complex((+:),phase)

infixl 6 .+
infixr 7 *>

integrate :: Num a => a -> [a] -> [a]
integrate = scanl (+)

(.+) :: Num a => [a] -> [a] -> [a]
(.+) = zipWith (+)

(*>) :: Num a => a -> [a] -> [a]
(*>) v = map (v*)

wave :: Num a => (a,a) -> (a,a) -> [a]
wave (k0,c0) (k1,c1) =
   let y'  = integrate c1 y''
       y   = integrate c0 y'
       y'' = map negate (k0 *> y  .+  k1 *> y')
   in  y

waveExample :: [Double]
waveExample = wave (0.07, 1) (0.08, 0)


waveSqr :: Num a => (a,a,a) -> (a,a) -> (a,a) -> [a]
waveSqr (a00,a01,a11) (k0,c0) (k1,c1) =
   let mul = zipWith (*)
       y'  = integrate c1 y''
       y   = integrate c0 y'
       y'' = map negate (foldl1 (.+)
               (zipWith (*>) [k0, k1, a00, a01, a11]
                             [y, y', mul y y, mul y y', mul y' y']))
   in  y

{- the square term destabilizes the solution -}
waveSqrExample :: [Double]
waveSqrExample = waveSqr (0.04,0,0) (0.07, 1) (0.08, 0)


waveSin :: Floating a => (a,a) -> (a,a) -> (a,a) -> [a]
waveSin (a0,a1) (k0,c0) (k1,c1) =
   let y'  = integrate c1 y''
       y   = integrate c0 y'
       y'' = map negate (foldl1 (.+)
               (zipWith (*>) [k0, k1, a0, a1]
                             [y, y', map sin y, map sin y']))
   in  y

{- the square term destabilizes the solution -}
waveSinExample :: [Double]
waveSinExample = waveSin (0.1,0) (0.07, 10) (0.08, 0)


wavePhase :: RealFloat a => a -> (a,a) -> (a,a) -> [a]
wavePhase (a0) (k0,c0) (k1,c1) =
   let y'  = integrate c1 y''
       y   = integrate c0 y'
       y'' = map negate (foldl1 (.+)
               (zipWith (*>) [k0, k1, a0]
                             [y, y', zipWith (\r i -> phase (r +: i)) y y']))
   in  y

{- the square term destabilizes the solution -}
wavePhaseExample :: [Double]
wavePhaseExample = wavePhase (0.005) (0.07, 1) (0.08, 0)