--
-- Explicit Euler method for 1D linear diffusion equation
-- (using lists)
--

module Euler1D ( stepEuler, evalEuler, makeu0 ) where

-- impose zero flux condition
zeroflux mu (boundary:inner:xs) = boundary+mu*2*(inner-boundary)

-- one step of integration
stepEuler mu u0 = (applyBC . (diffused mu)) $ u0
    where
          diffused mu (left:x:[]) = []    -- ignore outer points
          diffused mu (left:x:right:xs) = -- integrate inner points
                   let y = (x+mu*(left+right-2*x))
                       in y `seq` y : diffused mu (x:right:xs)
          applyBC inner = [lbc] ++ inner ++ [rbc] -- boundary conditions
               where
                     lbc = zeroflux mu ((head u0):(take 1 inner))           -- left boundary
                     rbc = zeroflux mu ((last u0):(take 1 $ reverse inner)) -- right boundary

-- initial condition
makeu0 n = [ f x n | x <- [0..n]]

f x n = ((^2) . sin . (pi*) . xi) x
    where xi x = fromIntegral x / fromIntegral n

-- integrate until time t with time step dt
evalEuler :: (Floating a, Ord a) => a -> a -> a -> [a] -> [a]
evalEuler dt t mu u
  | t <= 0.0 = u
  | otherwise =
      let prevstep = stepEuler mu u
      in  evalEuler (min dt (t-dt)) (t-dt) mu $! prevstep