Plotting the Mandelbrot Set

In this article, I'll show you how to draw the Mandelbrot set with a short Haskell program.

Prerequisites

You'll need a Haskell compiler/interpreter. For speed you'll probably want to use ghc, rather than hugs or ghci which will be very slow.

Ubuntu users should install the package ghc6.

The Code

We'll reuse the Image module from the ray tracer article. This will let us write PPM image files.

module Main where
import Complex
import Image

type CD = Complex Double
type MandelbrotState = (CD, CD)

-- The Mandlebrot sequence: x_{n+1} = x_n^2 + c_n, c_{n+1} = c_n.
mandelbrot_step :: MandelbrotState -> MandelbrotState
mandelbrot_step    (x, c)          =  (x * x + c, c)

-- Starting state for given position.
mandelbrot_initial :: CD -> MandelbrotState
mandelbrot_initial    z   = (z,z)

-- Terminate if the magnitude of x is bigger than two.
mandelbrot_termination_condition :: MandelbrotState -> Bool
mandelbrot_termination_condition = (> 2.0) . magnitude . fst

--f_index gives the least value of n s.t. termination_cond(f^n(initial))
--or max_n if no such value less than max_n exists
f_index :: (a-> a) -> (a -> Bool) -> a  -> Integer -> Integer
f_index    f          term_cond      state max_n   =  f_index' f term_cond
                                                                state max_n 0

f_index' :: (a-> a) -> (a -> Bool) -> a     -> Integer -> Integer -> Integer
f_index'    f          term_cond      state    max_n      n  
    | term_cond state = n
    | n == max_n      = n
    | otherwise       = f_index' f term_cond (f state) max_n (n+1)


-- We'll simplify the code for drawing by hiding all this from it - we'll just
-- give a simple complex -> integer function.
type DrawableFunction = CD -> Integer


-- stepper_to_drawable will convert to this function type.
stepper_to_drawable   :: (a-> a) -> (CD -> a) -> (a -> Bool) -> Integer
                         -> DrawableFunction
                         
stepper_to_drawable step init term max_n =  f
    where f z = ((f_index step term (init z) max_n) * 255) `div` max_n

-- int_to_colour will take an integer between 0 and 255 to a shade of grey
int_to_colour x = [x,x,x]

-- plot_function generates a list with one value per pixel
plot_function :: DrawableFunction -> Double -> Double -> Double -> Double
                 -> Integer -> Integer -> [Integer]

plot_function   f x_min x_max y_min y_max x_res y_res =
    concat [ [ f (x :+ y)
                          | x <- map fromPixelX [0..x_res - 1] ]
                          | y <- map fromPixelY [0..y_res - 1] ]
                         
    where fromPixelX x = x_min + (fromInteger x * x_scale)
          fromPixelY y = y_min + (fromInteger y * y_scale)
          x_scale = (x_max - x_min) / fromInteger (x_res - 1)
          y_scale = (y_max - y_min) / fromInteger (y_res - 1)


mandelbrot_f :: DrawableFunction
mandelbrot_f = stepper_to_drawable mandelbrot_step mandelbrot_initial
                                   mandelbrot_termination_condition 255
-- Test for mandlebrot_f
test_mandlebrot_f = mandlebrot_f 0 == 0 &&
                    mandlebrof f 4 == 255
                       
-- Finally, we bring it all together to draw the mandlebrot set
main = writeFile "mandelbrot.ppm" (create_ppm (x_res, y_res, pixels))
    where x_res  = 4096
          y_res  = 4096
          pixels = map int_to_colour (plot_function mandelbrot_f x_min x_max
                                                    y_min y_max x_res y_res)
          x_min  = -2.0
          x_max  =  2.0
          y_min  = -2.0
          y_max  =  2.0
         

Here's the kind of output you'll get (constrast upped and scaled down in GIMP).

Mandelbrot fractal

100% crop (constrast improved in GIMP):

Mandelbrot fractal (cropped)

Get the code.