-- |
-- 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. prefetch (used only for workpool skeletons)
--            4. Skeleton selection (0=farm (splitIntoN), 
--                                   1=offlineFarm (splitIntoN), 
--                                   2=workpool,
--                                   3=farm (unshuffle)
--                                   4=offlineFarm (unshuffle)
--
-- 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 Control.Parallel.Eden.Map
import Control.Parallel.Eden.Workpool
import Control.Parallel.Eden.Auxiliary 

import Data.Complex
import Data.List


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

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

           

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

main = do
  args <- getArgs
  if length args < 4 then putStrLn "Arguments missing..." else return ()
  let (lo, ru) = examples!!(read (args!!0))
  let dimx     = read (args!!1)
  let np       = noPe
  let pf       = read (args!!2)
  let skel     = if length args > 3 then read (args!!3) else 2
  let b        = bild 10.0 lo ru dimx np pf skel
  if length args > 4 && (args!!4) == "-out"
     then writeFile "out.ppm" b -- putStr b
     else rnf b `seq` putStrLn "Done (no output)"

examples = [
            ((-0.75104) :+ (0.10511), (-0.74080) :+ (0.11441)),
            ((-2.5)     :+ (-1.5),    (1.5)      :+ (1.5))
           ]

-- Compute image (bild) -------------------------------------------------

bild :: Double -> Complex Double -> Complex Double -> Integer -> Int -> Int -> Int -> 
        String
bild schwellwert lo ru dimx np pf s = 
  header ++ (  concat result  )
--
  where 
--  result :: [(Int, String)]
----------------------------------------------
	result = (skels!!s) tasks
	skels = cycle [farm (splitIntoN (np-1)) concat wf
		 , offlineFarm (np-1) (splitIntoN (np-1)) concat wf
		 , workpoolSorted (np-1) pf wf
		 , farm (unshuffle (np-1)) shuffle wf
		 , offlineFarm (np-1) (unshuffle (np-1)) shuffle wf
		 ]
----------------------------------------------
--        result'        = map snd result
--'
        tasks :: [[Complex Double]]
        tasks =  lines

        wf :: [Complex Double] -> String
        wf t      = concatMap (rgb . (iter schwellwert (0.0 :+ 0.0) 0)) t 

        header         = "P3\n"++(show (dimx+1))++" "++(show dimy)++"\n255\n"
        (dimy, lines)  = koord lo ru dimx
        rgb i          = i' ++ " " ++ i' ++ " " ++ i' ++ "\n"
          where i' = show i

-- Koordinaten aufspannen, Iteration -------------------------------------------

koord :: Complex Double -> Complex Double ->
         Integer ->
         (Integer, [[Complex Double]])
koord (x1 :+ y1) (x2 :+ y2) dimx = (dimy, ks)
  where breite   = abs (x2 - x1)
        hoehe    = abs (y2 - y1)
        schrittx = breite / ((fromInteger dimx)::Double)
        schritty = hoehe / ((fromInteger dimy)::Double)
        sx2      = schrittx / 2.0
        sy2      = schritty / 2.0
        dimy     = round ((((fromInteger dimx)::Double)*hoehe)/breite)
        ks       = [ [(x+sx2) :+ (y+sy2) | x <- [x1,x1+schrittx..x2]]
                     | y <- [y1,y1+schritty..y2]
                   ]

iter :: Double -> Complex Double -> Integer -> Complex Double -> Integer
iter schwellwert x it c | it > 255                  = 255
                        | (betrag x) >= schwellwert = it
                        | otherwise                 = iter schwellwert x' it' c
                                                      where it' = it + 1
                                                            x'  = x*x + c
betrag :: Complex Double -> Double
betrag (x :+ y) = sqrt (x*x + y*y)