module Main 
where
{-
  Visualisation of Julia sets

 With high resolutions output should be redirected:  
   #> program > image.ppm
 View the file (e.g. image.ppm) under Linux with 
   #> eog image.ppm

--------------------------------------------------------------------------- -}

import Data.Complex
import System.Environment(getArgs, getProgName)

import Control.Parallel.Eden
import Control.Parallel.Eden.Auxiliary
import Control.Parallel.Eden.Map
import Control.Parallel.Eden.Workpool
import Data.Time.Clock (diffUTCTime, getCurrentTime)

usage name = unlines $ 
	[""
	,"Visualisation for Julia sets"
	,""
	,"Usage: #>" ++ name ++ " <Version No> <Output (0/1)> <resolution(x)> <section (a to d / free)>"
	,""
        ,"Version 0: Sequential"
	]

main = do
        name <- getProgName
        args <- getArgs
        time1 <- getCurrentTime
        if (length args < 4)
            then putStrLn (usage name)
            else let version  = ((read.head) args)::Int
                     verbose  = 0 /= (read (args!!1))
                     dimx     = (read (args!!2))::Int
                     (lo,ru)  = frame (drop 3 args)
                     ppm = case version of
                             0 -> juliaSet juliastart threshold lo ru dimx
                             _ -> juliaSet juliastart threshold lo ru dimx
                     pixels = map rgb $ concat ppm                 -- add for coloring
                     header = "P3\n"++(show dimx)++" "++(show (length ppm))++"\n255\n"

                 in if verbose 
                       then putStr $ header ++ (unlines pixels)    --Output
                       else do
                             rnf ppm `seq` putStrLn "done"
                             time2 <- getCurrentTime
                             putStrLn $ "Version "++show version++" Runtime:"++ (show (diffUTCTime time2 time1))


-- Constants (as far as not program parameters)
threshold :: Double
threshold = 10.0
juliastart :: Complex Double
juliastart = 0.35 :+ 0.35

-- Section: ("upper left","down right")
frame :: [String] -> (Complex Double, Complex Double)
-- frame ("free":args) = ...
frame args = case lookup (head (head args)) flist of 
                    Just f  -> f
                    Nothing -> if length args > 4 
                                     then (free1,free2) 
                                     else frame ["c"]
        where flist = zip ['a'..]  
                       [((-0.75):+(0.11),(-0.74) :+ (0.10))
                       ,((-2.0) :+ (1.5) ,(2.0)  :+ (-1.5))
                       ,((-0.3) :+ (0.5) ,(0.7)  :+ (-0.5))
                       ,((-0.7) :+ (0.0) ,(-0.2) :+ (-1.0))
                        -- define further sections here ...
                       ]
              a     = (read (args!!1))::Double
              b     = (read (args!!2))::Double
              range = (read (args!!3))::Double
              free1 = a :+ b
              free2 = free1 + (range :+ (-range))

---------- (simple kind) Julia sets:  --------------------------

type JuliaSet = Complex Double -> Double -> Complex Double -> Complex Double -> Int -> [[Int]] 


--Julia Set
juliaSet :: JuliaSet
juliaSet cValue threshold ul lr dimx = pixels 


  where 
    -- Iteration for Julia-Set:
    iterJulia = \z -> iter threshold z 0 cValue
    pixels    =  map (map iterJulia . unbox) $ map Box coords
    -------------
    (dimy, coords) = coord ul lr dimx

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


-- color pixel 
rgb :: Int -> String
rgb i = unwords [show ri, show gi, show bi]
      where ri = (i*26) `mod` 255
            gi = (i*2) `mod` 255
            bi = (i*35) `mod` 255

-- determine coordinates for section 
coord :: Complex Double -> Complex Double ->
         Int ->
         (Int, [[Complex Double]])
coord (x1 :+ y1) (x2 :+ y2) dimx = (dimy, ks)
  where width    = abs (x2 - x1)
        height   = abs (y1 - y2)     -- y1 > y2 ("upper" > "lower")
        stepx    = width / (fromIntegral dimx)
        stepy    = height / dimy'
        sx2      = stepx / 2
        sy2      = stepy / 2
        dimy     = round dimy'
        dimy'    = ((fromIntegral dimx)::Double)*height/width
        ks       = map (\i -> mkxs (y1 -  (fromIntegral i) * stepy)) [0..dimy]
        mkxs y   = map (\i -> let x = x1 +  (fromIntegral i) * stepx
                              in  ((x+sx2) :+ (y-sy2))
                       ) [0..dimx] 
--        ks       = [    
--                    [ (x+sx2) :+ (y-sy2) | x <- [x1,x1+stepx..x2-stepx]
--                    ]
--                                         | y <- [y1,y1-stepy..y2+stepy]
--                   ]

-- Iteration fuerfor julia/mandelbrot
iter :: Double -> Complex Double -> Int -> Complex Double -> Int
iter threshold x it c 
  | it > 255                  = 255
  | (absValue x) >= threshold = it
  | absValue (x' - x) < 0.001 = 255
  | otherwise                 = iter threshold x' (it+1) c
                    where x'  = x*x + c

absValue :: Complex Double -> Double
absValue (x :+ y) = sqrt (x*x + y*y)

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

--Wrap type for unevaluated communication
newtype Box a = Box {unbox :: a}
instance NFData a =>  NFData (Box a) 
instance Trans a => Trans (Box a)