{-
  Visualisierung von Julia-Mengen
  Jost Berthold, Philipps-Universitaet Marburg, SS 2006 

 (Grundlage: mandelbrot.hs, Steffen Priebe WS 2000/2001)

 Bei groesseren Aufloesungen Ausgabe umleiten: 
   #> programm > bild.ppm
 Anzeige der Datei (z.B. bild.ppm) unter Linux mit 
   #> kview bild.ppm
       ODER
   #> xview bild.ppm
--------------------------------------------------------------------------- -}

import Complex
import System(getArgs, getProgName)

import Parallel
import Strategies
import U8 -- Listenstrategien, Aufgabe 8.x

usage name = unlines $ 
	[""
	,"Visualisierung fuer Julia-Mengen"
	,""
	,"Verwendung: #>" ++ name ++ " <Aufloesung(x)> <Ausschnitt (a bis d / frei)>"
	,""
	]

main = do
	name <- getProgName
        args <- getArgs
        if (length args < 2)
            then putStrLn (usage name)
            else let dimx    = (read (args!!0))::Int
                     (lo,ru) = frame (tail args)
		     ppm     = juliaSet juliastart threshold lo ru dimx
                 in putStr ppm

-- Konstanten (sofern nicht Programmparameter)
threshold :: Double
threshold = 10.0
juliastart :: Complex Double
juliastart = 0.35 :+ 0.35

-- Ausschnitt: ("links oben","rechts unten")
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))
			-- hier ggf. weitere Ausschnitte definieren...
		       ]
	      a     = (read (args!!1))::Double
	      b     = (read (args!!2))::Double
	      range = (read (args!!3))::Double
	      free1 = a :+ b
	      free2 = free1 + (range :+ (-range))


---------- Mandelbrotmenge:  --------------------------

mandelbrot :: Double -> Complex Double -> Complex Double -> Int -> String
mandelbrot schwellwert lo ru dimx = header ++ (unlines pixels)
    where pixels = map (rgb . (iter schwellwert (0.0 :+ 0.0) 0)) koords
          header         = "P3\n"++(show dimx)++" "++(show dimy)++"\n255\n" 
          (dimy, koords) = koord lo ru dimx 

---------- (einfache Art) Juliamengen:  --------------------------

juliaSet :: Complex Double -> Double -> Complex Double -> Complex Double -> Int -> String
juliaSet cWert schwellwert lo ru dimx = header ++ (unlines pixels)
  where 
    -- Iteration for Julia-Set:
    iterJulia      = \z -> iter schwellwert z 0 cWert
    pixels = map rgb (map iterJulia koords)
    -------------
    (dimy, koords) = koord lo ru dimx
    header = "P3\n"++(show dimx)++" "++(show dimy)++"\n255\n" 

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

-- Pixel faerben
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

-- Koordinatenmenge fuer den Ausschnitt bestimmen
koord :: Complex Double -> Complex Double ->
         Int ->
         (Int, [Complex Double])
koord (x1 :+ y1) (x2 :+ y2) dimx = (dimy, ks)
  where breite   = abs (x2 - x1)
        hoehe    = abs (y1 - y2)     -- y1 > y2 ("oben" > "unten")
        schrittx = breite / (fromIntegral dimx)
        schritty = hoehe / dimy'
        sx2      = schrittx / 2
        sy2      = schritty / 2
        dimy     = round dimy'
        dimy'    = ((fromIntegral dimx)::Double)*hoehe/breite
        ks       = [ (x+sx2) :+ (y-sy2) | y <- [y1,y1-schritty..y2+schritty],
                                          x <- [x1,x1+schrittx..x2-schrittx]]

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

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