--
-- mandelbrot.hs                                     Steffen Priebe WS 2000/2001
--
-- Bei groesseren Aufloesungen: a.out > bild.ppm
-- Anzeige der Datei bild.ppm unter Linux mit "xview bild.ppm"
--------------------------------------------------------------------------------

import Complex
import System(getArgs)

usage = "Usage: <program> <dimension> <which frame (a to c)>"

main = do
        args <- getArgs
        if (length args /= 2)
            then error usage 
            else do 
                  let dimx   = (read (args!!0))::Int
                      sw     = 10
                      (lo,ru)= frame (head (args!!1))
                      b      = bild sw lo ru dimx
                  putStr b

-- Ergebnis: ("links oben","rechts unten")
-- hier koennen weitere (moeglichst quadratische) Ausschnitte definiert werden:
frame :: Char -> (Complex Double, Complex Double)
frame c = (lo, ru)
	where lo = case c of 
		       'a' -> (-0.75104) :+ (0.11441)
	               'b' -> (-2.0) :+ (1.5)
		       'c' -> (-1.4) :+ (0.1)
                       _   -> error (c:" : invalid frame selection")
	      ru = case c of 
		       'a' -> (-0.7408) :+ (0.10511)
	               'b' -> (2.0) :+ (-1.5)
		       'c' -> (-1.2) :+ (-0.1)
{-  where dimx = 200
--        lo   = (-0.75104) :+ (0.10511)    -- Beispiel 1
--        ru   = (-0.74080) :+ (0.11441)
        lo   = (-2.0) :+ (-1.5)         -- Beispiel 2
        ru   =  (2.0) :+  (1.5)
        sw   = 10.0
        b    = bild sw lo ru dimx
-}

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

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]]

iter :: Double -> Complex Double -> Int -> Complex Double -> Int
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)