1 -- MANDELBROT SET GENERATOR
    2 -- DAVID HANLEY 11/03/1992
    3 -- BUILDS A MAND-TREE, which contains, at its leaves, a colour.
    4 -- Tree is then traversed and points along with their colour are output
    5 
    6 module Main where
    7 
    8 a `par` b = b
    9 --1.3 a `seq` b = b
   10 
   11 -- MandTree - either contains a NS node with two subtrees
   12 --                   contains a EW node with two subtrees
   13 --                   a leaf holding a colour
   14 data MandTree = NS   MandTree MandTree
   15                |EW   MandTree MandTree
   16                |Leaf Colour
   17 
   18 -- Colour - colour of a point or rectangle
   19 type Colour = Int
   20 
   21 -- Point - 2-Tuple of integers, (Int1,Int2),
   22 --         where Int1 - x-coord of point,
   23 --               Int2 - y-coord of point.
   24 type Point = (Int,Int)
   25 
   26 -- MandTuple
   27 type MandTuple = (Int,Int,Int,Int,Colour)
   28 
   29 -- size: Size of window - Square 400x400
   30 size :: Int
   31 size =  200
   32 
   33 -- build_tree - Constructs mandtree
   34 build_tree :: Point -> Point -> MandTree
   35 build_tree p1@(x1,y1)
   36            p2@(x2,y2) 
   37            = 
   38             if rec_col /= -1 then   -- All points in currnet rectangle are same colour
   39               Leaf rec_col
   40             else
   41               if split == "NS" then
   42                  btns1 `par` (btns2 `seq` (NS btns1 btns2))
   43               else
   44                  btew1 `par` (btew2 `seq` (EW btew1 btew2))
   45             where
   46                  rec_col = check_perim p1 p2
   47                  split   = if (x2-x1) >= (y2-y1) then  -- x-axis longer, NS split
   48                               "NS" 
   49                            else                        -- y-axis longer, EW split
   50                               "EW"
   51                  nsp1    = p1
   52                  nsp2    = (split_x, y2)
   53                  nsp3    = (split_x+1, y1)
   54                  nsp4    = p2
   55                  ewp1    = p1
   56                  ewp2    = (x2, split_y)
   57                  ewp3    = (x1, split_y+1)
   58                  ewp4    = p2
   59                  split_x = (x2+x1) `div` 2
   60                  split_y = (y2+y1) `div` 2
   61                  btns1   = build_tree nsp1 nsp2
   62                  btns2   = build_tree nsp3 nsp4
   63                  btew1   = build_tree ewp1 ewp2
   64                  btew2   = build_tree ewp3 ewp4
   65 
   66 
   67 check_perim :: Point -> Point -> Colour
   68 check_perim p1@(x1,y1)
   69             p2@(x2,y2)
   70             = if (equalp p1 p2) then  -- single point, just return its colour
   71                  point_colour p1
   72               else if corners_diff then  -- check corners of current rectangle aren't same
   73                  -1   
   74               else
   75                 check_sides
   76               where
   77               corners_diff = if col1 == col2 then
   78                                if col1 == col3 then
   79                                  if col1 == col4 then
   80                                    False
   81                                  else
   82                                    True
   83                                else
   84                                  True
   85                              else
   86                                True
   87 
   88               col1         = point_colour p1
   89               col2         = point_colour (x2,y1)
   90               col3         = point_colour (x2,y2)
   91               col4         = point_colour (x1,y2)
   92 
   93               check_sides  = if (check_line (x1+1,y1) right) then
   94                               if (check_line (x2,y1+1) down) then
   95                                if (check_line (x2-1,y2) left) then
   96                                 if (check_line (x1,y2-1) up) then
   97                                   col1
   98                                 else
   99                                   -1
  100                                else
  101                                  -1
  102                               else
  103                                -1
  104                              else
  105                               -1
  106 
  107                              where
  108                              check_line pc@(xc,yc) pd@(xd,yd)
  109                                =
  110                                if finished then
  111                                  True
  112                                else if (point_colour pc) /= col1 then
  113                                  False
  114                                else 
  115                                  check_line (xc+xd, yc+yd) pd
  116                                where
  117                                     finished  = if (equalp pd right) then
  118                                                   (xc >= x2)
  119                                                 else if (equalp pd down) then
  120                                                   (yc <= y2)
  121                                                 else if (equalp pd left) then
  122                                                   (xc <= x1)
  123                                                 else
  124                                                   (yc >= y1)
  125 
  126 
  127 -- Evaluate the color index of a point
  128 -- This is the algoritm described on page 121 of "The Beauty of Fractals"
  129 -- The code is commented with the step numbers from the algorithm.
  130 
  131 -- point_colour - Calculates the dwell value of a point.
  132 point_colour :: Point -> Colour
  133 point_colour (x, y)
  134              = check_radius (np x) (nq y) 0 0.0 0.0 -- step1
  135 
  136 
  137 -- check_radius - Calculates the escape radius of a point
  138 check_radius :: Float -> Float -> Int -> Float -> Float -> Colour
  139 check_radius p q k x y = if kp == num_cols then
  140                            0                             -- step 3.ii
  141                          else
  142                            if r > (fromIntegral m ) then
  143                             kp -- step 3.i
  144                            else
  145                              check_radius p q (kp) xn yn -- step 3.iii
  146                                   where xn = new_x x y p -- step 2
  147                                         yn = new_y x y q -- step 2
  148                                         r = radius xn yn -- step 3
  149                                         kp = k + 1       -- step 2
  150 
  151 
  152 -- M Set Properties.
  153  
  154 pmn :: Float -- Min p value.
  155 pmn =  -2.25
  156  
  157 pmx :: Float -- Max p value.
  158 pmx =   0.75
  159  
  160 qmn :: Float -- Min q value.
  161 qmn = -1.5
  162  
  163 qmx :: Float -- Max q value
  164 qmx = 1.5
  165 
  166 m :: Int      -- The escape radius, M.
  167 m =  20
  168  
  169 --- Misc functions.
  170  
  171 equalp :: Point -> Point -> Bool
  172 equalp (x1, y1) (x2, y2) = ((x1 == x2) && (y1 == y2))
  173  
  174  
  175 -- Set calculation functions.
  176  
  177 num_cols :: Int -- The number of colors; num_cols+1 = length (the_colors).
  178 num_cols = 26
  179  
  180 delta_p :: Float      -- The change in p per pixel.
  181 delta_p =  (pmx - pmn) / fromIntegral (size - 1)
  182 
  183 delta_q :: Float      -- The change in q per pixel.
  184 delta_q =  (qmx - qmn) / fromIntegral (size - 1)
  185  
  186 np :: Int -> Float     -- Calculate a new value of p.
  187 np x = pmn + (fromIntegral x) * delta_p 
  188  
  189 nq :: Int -> Float     -- Calculate a new value of q.
  190 nq y = qmn + (fromIntegral y) * delta_q
  191 
  192 radius :: Float -> Float -> Float       -- Current radius of apoint (x,y).
  193 radius x y = x*x + y*y
  194  
  195 new_x :: Float -> Float -> Float -> Float       -- iterate new x value.
  196 new_x x y p = x*x - y*y + p
  197  
  198 new_y :: Float -> Float -> Float -> Float       -- iterate new y value.
  199 new_y x y q = 2.0 * x * y + q
  200 
  201 
  202 -- Misc. functions for traversing perimeter of rectangle.
  203 
  204 up,down,left,right :: Point
  205 up    = ( 0,-1)  
  206 down  = ( 0, 1)
  207 left  = (-1, 0)
  208 right = ( 1, 0)
  209 
  210 -- finite: forces evaluation of a tree
  211 finite            :: MandTree -> Bool
  212 finite (Leaf c)   =  (c == c)
  213 finite (NS t1 t2) =  (finite t1 && finite t2)
  214 finite (EW t1 t2) =  (finite t1 && finite t2)
  215 
  216 main =  if finite(build_tree (0,0) (size,size `div` 2)) then
  217             print "Success"
  218           else
  219             print "Fail"
  220 
  221 
  222