Clifford Attractor

gistfile1.hs

import Data.Array.Repa ( Array, DIM2, DIM3, Z(..), (:.)(..) )
import qualified Data.Array.Repa          as R
import qualified Data.Array.Repa.IO.DevIL as D
import Data.Word  ( Word8   )
import Data.Fixed ( divMod' )
import Debug.Trace

import qualified Data.Map as M
import qualified Data.Vector.Unboxed as UV

traceS a = trace $ show a

type R     = Float
type R2    = (R, R)
type Angle = R

pixels :: Int
pixels = 640
pixels_half = div pixels 2
scale :: R
scale = 160

a = -1.4
b = 1.5
c = 1.0
d = 0.7
type R3 = (R,R,R)
start :: R3 
start = (0.156452879309654, 0.466939155012369, 0.144256137311459)
clifford_iter :: R -> R3 -> R3 
clifford_iter dt (x,y,z) = (dx, dy, dz)
  where dx = dt * (sin (a*y) - c * cos (a*x))
        dy = dt * (sin (b*x) - d * cos (b*y))
        dz = 0
clifford_num :: R -> R3 -> Int -> R3
clifford_num _ v 0 = v
clifford_num dt v n = clifford_num dt (clifford_iter dt v) (n-1)
clifford_expand :: R -> R3 -> Int -> [R3]
clifford_expand _ v 0 = [v]
clifford_expand dt v n = v : clifford_expand dt (clifford_iter dt v) (n-1)
project :: R3 -> Int
project (x,y,z) = 
  if vy >= pixels || vy < 0 then error $ "vy out of range (" ++ show x ++ " ," ++ show y ++ " ," ++ show z ++ ")" else
  vx * pixels + vy
  where vx = truncate (x*scale) + pixels_half -- + truncate (z * 0.1 * scale)
        vy = truncate (y*scale) + pixels_half -- + truncate (z * 0.1 * scale)
updateFromPositions :: UV.Vector R -> [R3] -> UV.Vector R
updateFromPositions arr ps =
  let alist = M.toList $ foldr (\p map->M.insertWith (+) (project p) 0.1 map) (M.empty) ps in
  arr UV.// alist

vectorToImage :: UV.Vector R -> Array DIM3 Word8
vectorToImage vec = R.traverse arr8 (:. 4) chans where
    arr  = R.fromVector (Z :. pixels :. pixels) vec
    arr8 = R.map (floor . (*255) . min 1 . max 0) arr
    chans _ (Z :. _ :. _ :. 3) = 255  -- alpha channel
    chans a (Z :. x :. y :. _) = a (Z :. x :. y)

main = do
  let trajectory = clifford_expand 1.0 start 1000000
  let whiteVector = UV.fromList $ (take (pixels*pixels) $ cycle [0])
  let vect = updateFromPositions whiteVector trajectory
  let image = vectorToImage vect
  D.runIL $ D.writeImage "outfrac.png" image

References:
as to using repa-devil library:
http://mainisusuallyafunction.blogspot.com/2011/10/quasicrystals-as-sums-of-waves-in-plane.html
about Clifford attractor:
http://paulbourke.net/fractals/clifford/