diff --git a/MuAtom.hs b/MuAtom.hs new file mode 100644 --- /dev/null +++ b/MuAtom.hs @@ -0,0 +1,115 @@ +{-# LANGUAGE BangPatterns, RecordWildCards #-} + +{- + + gmndl -- Mandelbrot Set explorer + Copyright (C) 2010 Claude Heiland-Allen + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +-} + +module MuAtom (muAtom) where + +import Data.Complex (Complex((:+)), cis, magnitude, mkPolar, phase) +import Data.Ratio (denominator) +import Numeric.GSL.Root (RootMethod(Hybrids)) +import qualified Numeric.GSL.Root as GSL + +type N = Integer +type Q = Rational +type R = Double +type C = Complex R + +-- a Mandelbrot Set mu-atom + +data Atom = Atom{ nucleus :: !C, period :: !N, root :: !C, cardioid :: !Bool } + deriving Show + +continent :: Atom +continent = Atom 0 1 1 True + +-- finding bond points + +f :: N -> C -> C -> C +f !n !z !c + | n == 0 = z + | otherwise = let w = f (n - 1) z c in w * w + c + +df :: N -> C -> C -> C +df !n !z !c = 2 ^ n * product [ f i z c | i <- [0 .. n - 1] ] + +bondIter :: N -> C -> [R] -> [R] +bondIter !n !b [x0, x1, x2, x3] = + let z = x0:+x1 + c = x2:+x3 + y0 :+ y1 = f n z c - z + y2 :+ y3 = df n z c - b + in [y0, y1, y2, y3] +bondIter _ _ _ = error "MuAtom.bondIter: internal error" + +-- finding nucleus + +l :: N -> C -> C +l !n !c + | n == 0 = 0 + | otherwise = let z = l (n - 1) c in z * z + c + +nucleusIter :: N -> [R] -> [R] +nucleusIter !n [x0, x1] = + let c = x0:+x1 + y0 :+ y1 = l n c + in [y0, y1] +nucleusIter _ _ = error "MuAtom.nucleusIter: internal error" + +-- finding descendants + +muChild :: Atom -> Q -> Atom +muChild Atom{..} address = + let -- some properties of the parent and its relation to the child + size = magnitude (root - nucleus) + angle = 2 * pi * realToFrac address + bondAngle = cis angle + child = period * denominator address + solve = GSL.root Hybrids 1e-10 10000 + -- perturb from the stable nucleus to help ensure convergence to the bond point + _initial@(ir :+ ii) = nucleus + mkPolar (size / 2) (phase (root - nucleus) + angle) + ([_, _, br, bi], _) = solve (bondIter period bondAngle) [ ir, ii, ir, ii ] + bondPoint = br :+ bi + -- estimate where the nucleus will be + radiusEstimate + | cardioid = size / m2 * sin (pi * realToFrac address) + | otherwise = size / m2 + where m2 = fromIntegral (denominator address) ^ (2 :: N) + deltaEstimate = bondPoint - nucleus + _guess@(gr :+ gi) = bondPoint + mkPolar radiusEstimate (phase deltaEstimate) + -- refine the nucleus estimate + ([cr, ci], _) = solve (nucleusIter child) [gr, gi] + childNucleus = cr :+ ci + in Atom childNucleus child bondPoint False + +muChildren :: Atom -> [Q] -> [Atom] +muChildren a [] = [a] +muChildren a (q:qs) = let b = muChild a q in a : muChildren b qs + +-- interface to the outside world + +muAtom :: [Rational] -> (Double, Double, Double, Integer) +muAtom qs = + let Atom{..} = last $ muChildren continent qs + r :+ i = nucleus + s = magnitude (nucleus - root) + p = period + in (r, i, s, p) diff --git a/gmndl.cabal b/gmndl.cabal --- a/gmndl.cabal +++ b/gmndl.cabal @@ -1,14 +1,20 @@ Name: gmndl -Version: 0.2 +Version: 0.3 Synopsis: Mandelbrot Set explorer using GTK Description: - A Mandelbrot Set explorer, with (currently) only one - control (left-click to center and zoom in). Multple - render threads use higher precision maths at higher - zoom levels. Suggested usage: + A Mandelbrot Set explorer. Multiple render threads + use higher precision maths at higher zoom levels. + Suggested usage: . @gmndl +RTS -N -qa -RTS --width=640 --height=480@ + . + Left-click to zoom in, right-click to zoom out. The + status bar shows where you are, and the entry field + takes a space-separated list of fractions strictly + between 0 and 1, try for example: + . + @1\/2 2\/3 1\/4 3\/5@ Cabal-version: >=1.4 License: GPL-2 @@ -20,10 +26,12 @@ Executable gmndl Main-is: gmndl.hs + Other-modules: MuAtom Build-depends: base >= 4 && < 5, array >= 0.3 && < 0.4, gtk >= 0.11 && < 0.13, gtkglext >= 0.11 && < 0.13, + hmatrix >= 0.10 && < 0.11, mtl, OpenGL >= 2.4 && < 2.5, priority-queue >= 0.2.1 && < 0.3, diff --git a/gmndl.hs b/gmndl.hs --- a/gmndl.hs +++ b/gmndl.hs @@ -1,4 +1,26 @@ {-# LANGUAGE BangPatterns, ForeignFunctionInterface #-} + +{- + + gmndl -- Mandelbrot Set explorer + Copyright (C) 2010 Claude Heiland-Allen + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +-} + module Main (main) where -- some simple helpers @@ -6,7 +28,7 @@ import Data.List (isPrefixOf) -- concurrent renderer with capability-specific scheduling -import Control.Concurrent (killThread) +import Control.Concurrent (forkIO, killThread) import GHC.Conc (forkOnIO, numCapabilities) -- the dependency on mtl is just for this! @@ -46,6 +68,9 @@ -- ugly! but the default realToFrac :: (C)Double -> (C)Double is slooow import Unsafe.Coerce (unsafeCoerce) +-- mu-atom properties +import MuAtom (muAtom) + -- some type aliases to shorten things type B = Word8 type N = Int @@ -330,7 +355,8 @@ widgetSetSizeRequest canvas width height -- allocate some image bytes and clear with white imgdata <- mallocBytes $ width * height * 3 - _ <- memset (castPtr imgdata) 255 (fromIntegral $ height * width * 3) + let clear = memset (castPtr imgdata) 255 (fromIntegral $ height * width * 3) >> return () + clear let output x y r g b = do let p = imgdata `plusPtr` ((y * width + x) * 3) pokeByteOff p 0 r @@ -339,15 +365,17 @@ window <- windowNew eventb <- eventBoxNew vbox <- vBoxNew False 0 - hbox <- hBoxNew False 0 + status <- vBoxNew False 0 statusRe <- labelNew Nothing statusIm <- labelNew Nothing statusZo <- labelNew Nothing + ratios <- entryNew boxPackStart vbox eventb PackGrow 0 - boxPackStart vbox hbox PackGrow 0 - boxPackStart hbox statusRe PackGrow 0 - boxPackStart hbox statusIm PackGrow 0 - boxPackStart hbox statusZo PackGrow 0 + boxPackStart vbox status PackGrow 0 + boxPackStart vbox ratios PackGrow 0 + boxPackStart status statusRe PackGrow 0 + boxPackStart status statusIm PackGrow 0 + boxPackStart status statusZo PackGrow 0 let updateStatus re im zo = do -- this updates the status bar labelSetText statusRe (reshow $ show re) labelSetText statusIm (reshow $ show im) @@ -372,17 +400,40 @@ updateStatus re' im' zoom -- when the mouse button is pressed, center and zoom in _ <- eventb `on` buttonPressEvent $ {-# SCC "cbEv" #-} tryEvent $ do - LeftButton <- eventButton + b <- eventButton (x, y) <- eventCoordinates liftIO $ do (c, zoom, stop) <- readIORef sR stop - _ <- memset (castPtr imgdata) 255 (fromIntegral $ height * width * 3) + clear let c'@(re' :+ im') = c + ((convert x :+ convert (-y)) - (fromIntegral width :+ fromIntegral (-height)) * (0.5 :+ 0)) * ((1/2^^zoom) :+ 0) - zoom' = zoom + 1 + zoom' = zoom + delta + delta | b == LeftButton = 1 + | b == RightButton = -1 + | otherwise = 0 stop' <- renderer' rng output c' zoom' - writeIORef sR (c', zoom', stop') + writeIORef sR $! (c', zoom', stop') updateStatus re' im' zoom' + -- when pressing return in the ratios list, zoom to that mu-atom + _ <- ratios `onEntryActivate` do + s <- entryGetText ratios + case rationalize s of + Nothing -> return () + Just qs -> do + _ <- ratios `widgetSetSensitive` False + _ <- forkIO $ do + let (cr, ci, radius, _period) = muAtom qs + zoom' = floor $ (logBase 2 . fromIntegral $ width `min` height) - (logBase 2 radius) - 2 + c'@(cr' :+ ci') = convert cr :+ convert ci + cr `seq` ci `seq` zoom' `seq` postGUISync $ do + (_, _, stop) <- readIORef sR + stop + clear + stop' <- renderer' rng output c' zoom' + writeIORef sR $! (c', zoom', stop') + _ <- ratios `widgetSetSensitive` True + updateStatus cr' ci' zoom' + return () -- need to set up OpenGL stuff in callback just because that's how... _ <- onRealize canvas $ {-# SCC "cbRz" #-} withGLDrawingArea canvas $ \_ -> do GL.clearColor $= (GL.Color4 0.0 0.0 0.0 0.0) @@ -488,3 +539,14 @@ [] -> "0." ++ y x' -> x' ++ "." ++ y in reverse . dropWhile (=='.') . dropWhile (=='0') . reverse . (if sign then ('-':) else id) $ digits + +-- convert a string of space separated fractions into rationals +-- moreover ensure that they are all in 0 Maybe [Rational] +rationalize s = + let f '/' = '%' + f c = c + r t = case reads t of + [(q, "")] | 0 < q && q < 1 -> Just q + _ -> Nothing + in mapM r . words . map f $ s