Tuesday, December 22, 2015

lyapunov fractals

It took me some experimentation to figure out how to color this derivation of the logistic map; I'm still not quite sure how the hues should scale as you zoom. But the bi-tonal method shown below works well enough to produce the image at right.
{-# LANGUAGE NoMonomorphismRestriction #-}

import Control.Monad (liftM2)
import Data.Bits (shiftL)
import Graphics.UI.SDL as SDL
 
(xres, yres, cast) = (400, 400, fromIntegral)

main = withInit [InitVideo] $ do
  w <- setVideoMode xres yres 32 [NoFrame]
  enableEvent SDLMouseMotion False
  setCaption "Lyapunov" "Lyapunov"
  render w 1024 (3.2, 3.2) (3.8, 3.8)

render w i (x,y) (s,t) = do
  draw w (zip ps $ area x s y t) i
  run w
 where
  ps = liftM2 (,) [0.. xres] [0.. yres]

draw w cs i = do
  s <- createRGBSurface [] 1 1 32 0 0 0 0
  mapM_ (draw s) cs
  SDL.flip w
 where
  draw s (p,q) = do 
    let rect (x,y) = Just $ Rect x y 1 1
    fillRect s Nothing $ rgb $ lyapunov i q
    blitSurface s Nothing w $ rect p

run w = do
  delay 32
  e <- pollEvent
  case e of
   KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
   _                              -> run w

-----------------------------------------------------

rgb n | n < 0 = Pixel $ f n
      | True  = Pixel $ f n `shiftL` 16
 where
  f x = round (min 255 $ abs n * 255)

lyapunov i p = (1 / cast i) * sum rs
 where
  rs = filter (not . isInfinite) $ f xs
  f  = map (log . abs) . zipWith ($) (str p)
  xs = map ((1 -) . (* 2)) $ verhulst i p

verhulst i p = take i $ go 0.5 $ str p
 where
  go x (f:fs) = x : go (f (x * (1-x))) fs

str (x,y) = map (*) $ concat $ repeat [x,x,y,x,y]

area x x' y y' = liftM2 (,) [s, s+a.. t] [u, u+b.. v]
 where
  (a,b) = ((t - s) / m, (v - u) / n)
  (s,t) = (min x x', max x x')
  (u,v) = (min y y', max y y')
  (m,n) = (cast xres, cast yres)

No comments:

Post a Comment