gistfile1.hs

linSolveStep :: Int -> Grid -> Grid -> Double -> Double -> IO ()
linSolveStep b x x0 a c = forEachPixel x
                            (\(i,j) -> 
                             do
                               (up,down,left,right) <- neighbours x (i,j)
                               x0v <- readVal x0 (i,j)
                               writeVal x (i,j) ((x0v + a*(up + down + left + right)) / c))
                          >> setBnd b x

linSolve :: Int -> Grid -> Grid -> Double -> Double -> IO()
linSolve b x x0 a c = forM_ [1..20] (\_ -> linSolveStep b x x0 a c)

diffuse :: Int -> Grid -> Grid -> Double -> Double -> IO()
diffuse b x@(Grid n _) x0 diff dt = linSolve b x x0 a (1 + 4*a) where
    a = dt * diff * (fromIntegral n * fromIntegral n)