gmndl 0.2 → 0.3
raw patch · 3 files changed
+201/−16 lines, 3 filesdep +hmatrix
Dependencies added: hmatrix
Files
- MuAtom.hs +115/−0
- gmndl.cabal +13/−5
- gmndl.hs +73/−11
+ MuAtom.hs view
@@ -0,0 +1,115 @@+{-# LANGUAGE BangPatterns, RecordWildCards #-}++{-++ gmndl -- Mandelbrot Set explorer+ Copyright (C) 2010 Claude Heiland-Allen <claudiusmaximus@goto10.org>++ 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)
gmndl.cabal view
@@ -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,
gmndl.hs view
@@ -1,4 +1,26 @@ {-# LANGUAGE BangPatterns, ForeignFunctionInterface #-}++{-++ gmndl -- Mandelbrot Set explorer+ Copyright (C) 2010 Claude Heiland-Allen <claudiusmaximus@goto10.org>++ 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<q<1+rationalize :: String -> 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