diff --git a/Fractal/RUFF/Mandelbrot/Atom.hs b/Fractal/RUFF/Mandelbrot/Atom.hs new file mode 100644 --- /dev/null +++ b/Fractal/RUFF/Mandelbrot/Atom.hs @@ -0,0 +1,241 @@ +{-# LANGUAGE BangPatterns #-} +{- | +Module : Fractal.RUFF.Mandelbrot.Atom +Copyright : (c) Claude Heiland-Allen 2011 +License : BSD3 + +Maintainer : claudiusmaximus@goto10.org +Stability : unstable +Portability : portable + +Mu-atom coordinate and address algorithms. +-} +module Fractal.RUFF.Mandelbrot.Atom + ( MuAtom(..) + , FindAtom(..), findAtom, findAtom', findAtom_ + , FindAddress(..), findAddress, findAddress', findAddress_ + , Locate(..), locate, locate', locate_ + ) where + +import Control.Arrow ((***)) +import Data.Maybe (listToMaybe) +import Data.Ratio ((%)) +import Data.Vec (NearZero, nearZero) + +import Fractal.RUFF.Mandelbrot.Address (AngledInternalAddress, Angle, splitAddress, addressPeriod, externalAngles, angledInternalAddress) +import Fractal.RUFF.Mandelbrot.Nucleus (findNucleus, findBond, findPeriod) +import Fractal.RUFF.Mandelbrot.Ray (externalRay, externalRayOut) +import Fractal.RUFF.Types.Complex (Complex, magnitude, magnitude2, phase, mkPolar) + +-- | Mu-atom properties. +data MuAtom r = MuAtom + { muNucleus :: !(Complex r) + , muSize :: !Double + , muOrient :: !Double + , muPeriod :: !Integer + } + deriving (Read, Show, Eq) + +-- | Progress updates for 'findAtom'. +data FindAtom r + = AtomSplitTodo + | AtomSplitDone AngledInternalAddress [Angle] + | AtomAnglesTodo + | AtomAnglesDone !Rational !Rational + | AtomRayTodo + | AtomRay !Integer + | AtomRayDone !(Complex r) + | AtomNucleusTodo + | AtomNucleus !Integer + | AtomNucleusDone !(Complex r) + | AtomBondTodo + | AtomBond !Integer + | AtomBondDone !(Complex r) + | AtomSuccess !(MuAtom r) + | AtomFailed + deriving (Read, Show, Eq) + +isAtomSuccess :: FindAtom r -> Bool +isAtomSuccess (AtomSuccess _) = True +isAtomSuccess _ = False + +fromAtomSuccess :: FindAtom r -> Maybe (MuAtom r) +fromAtomSuccess (AtomSuccess m) = Just m +fromAtomSuccess _ = Nothing + +-- | Try to find an atom, providing progress updates. +findAtom :: (Floating r, NearZero r, Real r) => AngledInternalAddress -> [FindAtom r] +findAtom addr = AtomSplitTodo : + let (!iaddr, !caddr) = splitAddress addr + !p = addressPeriod iaddr + in AtomSplitDone iaddr caddr : AtomAnglesTodo : case externalAngles iaddr of + Nothing -> [AtomFailed] + Just (!lo, !hi) -> AtomAnglesDone lo hi : AtomRayTodo : + let sharpness = 8 + er = 65536 + accuracy = 1e-10 + ok w = magnitude2 w < 2 * er ^ (2::Int) -- NaN -> False + rayl = externalRay accuracy sharpness er lo + rayh = externalRay accuracy sharpness er hi + ray' = takeWhile (uncurry (&&) . (ok *** ok) . snd) $ [ 1 .. ] `zip` (rayl `zip` rayh) + rgo [] _ = [AtomFailed] + rgo [_] _ = [AtomFailed] + rgo ((i, (xl, xh)):m@((_, (yl, yh)):_)) f + | i < fromIntegral sharpness * (p + 16) = AtomRay i : rgo m f + | dl + dh > dx + dy = AtomRay i : rgo m f + | otherwise = AtomRayDone x : f x + where + x = 0.5 * (xl + xh) + dl = magnitude2 (xl - yl) + dh = magnitude2 (xh - yh) + dx = magnitude2 (xl - xh) + dy = magnitude2 (yl - yh) + in rgo ray' $ \rayend -> AtomNucleusTodo : + let nuc = findNucleus p rayend + nuc' = takeWhile (ok . snd) $ [ 1 .. ] `zip` nuc + ngo [] _ = [AtomFailed] + ngo [_] _ = [AtomFailed] + ngo ((i, x):m@((_, y):_)) f + | not (nearZero (x - y)) = AtomNucleus i : ngo m f + | otherwise = AtomNucleusDone x : f x + in ngo nuc' $ \nucleus -> AtomBondTodo : + let bnd = findBond p nucleus 0.5 + bnd' = takeWhile (ok . snd) $ [ 1.. ] `zip` bnd + bgo [] _ = [AtomFailed] + bgo [_] _ = [AtomFailed] + bgo ((i, x):m@((_, y):_)) f + | not (nearZero (x - y)) = AtomBond i : bgo m f + | otherwise = AtomBondDone x : f x + in bgo bnd' $ \bond -> + let delta = bond - nucleus + size = realToFrac $ magnitude delta / 0.75 -- FIXME check for island-hood + orient = realToFrac $ phase delta + atom = MuAtom{ muNucleus = nucleus, muSize = size, muOrient = orient, muPeriod = p } + in if 10 > size && size > 0 then [AtomSuccess atom] else [AtomFailed] + +-- | Find the first success in the progress list. +findAtom' :: [FindAtom r] -> Maybe (MuAtom r) +findAtom' ps = fromAtomSuccess =<< listToMaybe (filter isAtomSuccess ps) + +-- | Find an atom from its address. +findAtom_ :: (Floating r, NearZero r, Real r) => AngledInternalAddress -> Maybe (MuAtom r) +findAtom_ = findAtom' . findAtom + +-- | Progress updates for 'findAddress'. +data FindAddress r + = AddressCuspTodo + | AddressCuspDone !(Complex r) + | AddressDwellTodo + | AddressDwell !Integer + | AddressDwellDone !Integer + | AddressRayOutTodo + | AddressRayOut !Double + | AddressRayOutDone !(Complex r) + | AddressExternalTodo + | AddressExternalDone !Rational + | AddressAddressTodo + | AddressSuccess AngledInternalAddress + | AddressFailed + deriving (Read, Show, Eq) + +isAddressSuccess :: FindAddress r -> Bool +isAddressSuccess (AddressSuccess _) = True +isAddressSuccess _ = False + +fromAddressSuccess :: FindAddress r -> Maybe AngledInternalAddress +fromAddressSuccess (AddressSuccess a) = Just a +fromAddressSuccess _ = Nothing + +-- | Try to find an address, providing progress updates. +findAddress :: (Floating r, NearZero r, Real r, RealFrac r) => MuAtom r -> [FindAddress r] +findAddress mu = AddressCuspTodo : + let cusp = muNucleus mu - mkPolar (realToFrac (muSize mu)) (realToFrac ((muOrient mu))) + er = 65536 + er2 = er * er + in AddressCuspDone cusp : AddressDwellTodo : + let dgo z n f = AddressDwell n : if magnitude2 z > er2 then f n else dgo (z * z + cusp) (n + 1) f + in dgo 0 0 $ \n -> AddressDwellDone n : AddressRayOutTodo : + let rgo ((i,!_):izs@(_:_)) f = AddressRayOut (fromIntegral i / (fromIntegral sharpness * fromIntegral n)) : rgo izs f + rgo [(_,!z)] f | magnitude2 z > er2 = AddressRayOutDone z : f z + rgo _ _ = [AddressFailed] + accuracy = 1e-16 + sharpness = 16 + epsilon0 = realToFrac (muSize mu) * accuracy + in rgo ([(1 :: Integer) ..] `zip` externalRayOut (fromIntegral n + 100) epsilon0 accuracy sharpness er cusp) $ \rend -> AddressExternalTodo : + let den = 2 ^ muPeriod mu - 1 + num' = fromIntegral den * warp (phase rend / (2 * pi)) + num = round num' + warp t + | t > 0 = t + | otherwise = t + 1 + angle = num % den + in AddressExternalDone angle : AddressAddressTodo : case angledInternalAddress angle of + Nothing -> [AddressFailed] + Just addr -> if addressPeriod addr /= muPeriod mu then [AddressFailed] else [AddressSuccess addr] + +-- | Find the first success in the progress list. +findAddress' :: [FindAddress r] -> Maybe AngledInternalAddress +findAddress' ps = fromAddressSuccess =<< listToMaybe (filter isAddressSuccess ps) + +-- | Find an address for a mu-atom. +findAddress_ :: (Floating r, NearZero r, Real r, RealFrac r) => MuAtom r -> Maybe AngledInternalAddress +findAddress_ = findAddress' . findAddress + +-- | Progress updates for 'locate'. +data Locate r + = LocateScanTodo + | LocateScan + | LocateScanDone !Integer + | LocateNucleusTodo + | LocateNucleus !Integer + | LocateNucleusDone !(Complex r) + | LocateBondTodo + | LocateBond !Integer + | LocateBondDone !(Complex r) + | LocateSuccess !(MuAtom r) + | LocateFailed + deriving (Read, Show, Eq) + +isLocateSuccess :: Locate r -> Bool +isLocateSuccess (LocateSuccess _) = True +isLocateSuccess _ = False + +fromLocateSuccess :: Locate r -> Maybe (MuAtom r) +fromLocateSuccess (LocateSuccess m) = Just m +fromLocateSuccess _ = Nothing + +-- | Try to find an atom close to a coordinate, providing progress updates. +locate :: (Floating r, NearZero r, Real r) => Complex r {- ^ center -} -> Double {- ^ radius -} -> [Locate r] +locate c r = LocateScanTodo : LocateScan : case findPeriod 10000000 (realToFrac r) c of -- FIXME hardcoded + Nothing -> [LocateFailed] + Just p -> LocateScanDone p : LocateNucleusTodo : + let ok w = magnitude2 w < 16 -- NaN -> False + nuc = findNucleus p c + nuc' = takeWhile (ok . snd) $ [ 1 .. ] `zip` nuc + ngo [] _ = [LocateFailed] + ngo [_] _ = [LocateFailed] + ngo ((i, x):m@((_, y):_)) f + | not (nearZero (x - y)) = LocateNucleus i : ngo m f + | otherwise = LocateNucleusDone x : f x + in ngo nuc' $ \nucleus ->LocateBondTodo : + let bnd = findBond p nucleus 0.5 + bnd' = takeWhile (ok . snd) $ [ 1.. ] `zip` bnd + bgo [] _ = [LocateFailed] + bgo [_] _ = [LocateFailed] + bgo ((i, x):m@((_, y):_)) f + | not (nearZero (x - y)) = LocateBond i : bgo m f + | otherwise = LocateBondDone x : f x + in bgo bnd' $ \bond -> + let delta = bond - nucleus + size = realToFrac $ magnitude delta / 0.75 -- FIXME only valid for islands + orient = realToFrac $ phase delta + atom = MuAtom{ muNucleus = nucleus, muSize = size, muOrient = orient, muPeriod = p } + in if 10 > size && size > 0 then [LocateSuccess atom] else [LocateFailed] + +-- | Find the first success in the progress list. +locate' :: [Locate r] -> Maybe (MuAtom r) +locate' ps = fromLocateSuccess =<< listToMaybe (filter isLocateSuccess ps) + +-- | Find an atom close to a coordinate. +locate_ :: (Floating r, NearZero r, Real r) => Complex r -> Double -> Maybe (MuAtom r) +locate_ c r = locate' (locate c r) diff --git a/Fractal/RUFF/Types/Complex.hs b/Fractal/RUFF/Types/Complex.hs --- a/Fractal/RUFF/Types/Complex.hs +++ b/Fractal/RUFF/Types/Complex.hs @@ -19,6 +19,7 @@ import Data.Data (Data) import Data.Typeable (Typeable) +import Data.Vec (NearZero(nearZero)) -- | Complex number type without the 'RealFloat' constraint. data Complex r = !r :+ !r @@ -53,6 +54,9 @@ asinh z = log (z + sqrt (z*z + 1)) acosh z = log (z + sqrt (z*z - 1)) atanh z = 1/2 * log ((1 + z) / (1 - z)) + +instance NearZero r => NearZero (Complex r) where + nearZero (r :+ i) = nearZero r && nearZero i -- | Extract the real part. realPart :: Complex r -> r diff --git a/ruff.cabal b/ruff.cabal --- a/ruff.cabal +++ b/ruff.cabal @@ -1,5 +1,5 @@ Name: ruff -Version: 0.2 +Version: 0.3 Synopsis: relatively useful fractal functions Description: A library for analysis and exploration of fractals, providing @@ -16,10 +16,11 @@ Category: Math Build-type: Simple -Cabal-version: >=1.2 +Cabal-version: >=1.6 Library Exposed-modules: Fractal.RUFF.Mandelbrot.Address + Fractal.RUFF.Mandelbrot.Atom Fractal.RUFF.Mandelbrot.Image Fractal.RUFF.Mandelbrot.Iterate Fractal.RUFF.Mandelbrot.Nucleus @@ -29,6 +30,16 @@ Build-depends: base >= 3 && < 6, array >= 0.3 && < 0.4, mtl >= 2 && < 3, - parsec >= 3.1 && < 3.2 + parsec >= 3.1 && < 3.2, + Vec >= 0.9 && < 1 GHC-Options: -Wall -O2 GHC-Prof-Options: -prof -auto-all -caf-all + +source-repository head + type: git + location: git://gitorious.org/ruff/ruff.git + +source-repository this + type: git + location: git://gitorious.org/ruff/ruff.git + tag: v0.3