-- |
-- Module      :  Mandelbrot
-- Copyright   :  (c) Philipps Universitaet Marburg 2009-2010
-- License     :  BSD-style (see the file LICENSE)
-- 
-- Maintainer  :  eden@mathematik.uni-marburg.de
-- Stability   :  beta
-- Portability :  not portable
--
-- The following Haskell module implements Mandelbrot fractals with Eden.
--
-- Depends on the Eden Compiler.
--
-- Eden Project



--------------------------------------------------------------------------------
--
-- Mandelbrot fractal computation
--
-- Arguments: 1. Example-coordinates (0 or 1)   
--            2. number of pixels per image row   
--            3. Skeleton selection (0=farm (unshuffle),
--                                   1=farm (splitIntoN), 
--                                   2=Workpool (erroneous!)
--
-- Output:    PPM-Format to Stdout, only if 5th argument is "-out"
--            otherwise output is suppressed, return value: "Done"
--
--------------------------------------------------------------------------------



import System.Environment
import Control.Seq
import Control.Parallel (pseq)
import Control.Parallel.Eden

import Data.Complex
import Data.List

import ParMap

usage = "Mandelbrot fractal computation\n" ++
      " Arguments: 1. Example-coordinates (0 or 1)\n" ++   
      "            2. number of pixels per image row\n" ++   
      "            3. Skeleton selection (incomplete!)\n" ++
      "\n" ++
      " Output:    PPM-Format to Stdout, only if 4th argument is \"-out\"\n" ++
      "            otherwise output is suppressed, return value: \"Done\"\n"

-- Parallel --------------------------------------------------------------------

-- instance NFData a => NFData (Complex a)
instance (RealFloat a, Trans a)  => Trans  (Complex a)

-- Main ------------------------------------------------------------------------

main = do
  args <- getArgs
  if length args < 3
     then putStrLn $ "Arguments missing...\n" ++ usage 
     else do
      let (ul, lr) = examples!!(read (args!!0))
      let dimx     = read (args!!1)
      let np       = noPe
      let skel     = if length args > 2 then read (args!!2) else 2
      let b        = pic 5.0 ul lr dimx np skel
      if length args > 3 && (args!!3) == "-out"
     	then putStr b
     	else rnf b `seq` putStrLn "Done (no output)"

-- examples, given by upper left and lower right corner
examples = [
            ((-0.75104) :+ (0.10511), (-0.74080) :+ (0.11441)),
            ((-2.5)     :+ (-1.5),    (1.5)      :+ (1.5))
           ]

-- Compute image (pic) -------------------------------------------------

pic :: Double -> Complex Double -> Complex Double -> Integer -> Int -> Int ->
       String
pic threshold ul lr dimx np s = header ++ concat (parMap computeRow tasks)
  where parMap = (cycle [map, parMapEden, farm0, parMapWP] )!!s
        --------------------------------------------------------
        tasks :: [[Complex Double]]
        tasks =  lines

        computeRow :: [Complex Double] -> String
        computeRow t = concatMap (rgb . (iter 255 threshold (0.0 :+ 0.0) 0)) t

        header         = "P3\n"++(show (dimx+1))++" "++(show dimy)++"\n255\n"
        (dimy, lines)  = coord ul lr dimx

-- rgb mapping for pixel values
rgb :: Integer -> String
rgb i = (show ri) ++ " " ++ (show gi) ++ " " ++ (show bi) ++ "\n" 
      where ri = (i*26) `mod` 255
	    gi = (i*2) `mod` 255
	    bi = (i*35) `mod` 255
	  {-gi | i < 170 = (i `mod` 85) * 3
	       | otherwise = 255
	    bi | i < 170 = 0
	       | otherwise = (i `mod` 170)*3-}

-- Coordinates and iteration -------------------------------------------

coord :: Complex Double -> Complex Double ->
         Integer ->
         (Integer, [[Complex Double]])
coord (x1 :+ y1) (x2 :+ y2) dimx = (dimy, ks)
  where breadth   = abs (x2 - x1)
        height    = abs (y2 - y1)
        stepx     = breadth / ((fromInteger dimx)::Double)
        stepy     = height / ((fromInteger dimy)::Double)
        sx2       = stepx / 2.0
        sy2       = stepy / 2.0
        dimy     = round ((((fromInteger dimx)::Double)*height)/breadth)
        ks       = [ [(x+sx2) :+ (y+sy2) | x <- [x1,x1+stepx..x2]]
                     | y <- [y1,y1+stepy..y2]
                   ]

iter :: Integer -> Double -> Complex Double -> Integer -> Complex Double -> 
        Integer
iter limit threshold x it c | it > limit                  = limit
                            | (magnitude x) >= threshold = it
                            | otherwise                 
                                = iter limit threshold (x*x+c) (it+1) c