ruff-0.4: Fractal/RUFF/Mandelbrot/Atom.hs
{-# LANGUAGE BangPatterns #-}
{- |
Module : Fractal.RUFF.Mandelbrot.Atom
Copyright : (c) Claude Heiland-Allen 2011,2015
License : BSD3
Maintainer : claude@mathr.co.uk
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.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)
import Fractal.RUFF.Types.Ratio ((%), fromQ)
-- | Mu-atom properties.
data MuAtom r = MuAtom
{ muNucleus :: !(Complex r)
, muSize :: !Double
, muOrient :: !Double
, muPeriod :: !Int
}
deriving (Read, Show, Eq)
-- | Progress updates for 'findAtom'.
data FindAtom r
= AtomSplitTodo
| AtomSplitDone AngledInternalAddress [Angle]
| AtomAnglesTodo
| AtomAnglesDone !Angle !Angle
| AtomRayTodo
| AtomRay !Int
| AtomRayDone !(Complex r)
| AtomNucleusTodo
| AtomNucleus !Int
| AtomNucleusDone !(Complex r)
| AtomBondTodo
| AtomBond !Int
| 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 (fromQ lo)
rayh = externalRay accuracy sharpness er (fromQ 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 !Int
| AddressDwellDone !Int
| AddressRayOutTodo
| AddressRayOut !Double
| AddressRayOutDone !(Complex r)
| AddressExternalTodo
| AddressExternalDone !Angle
| 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 :: Int) ..] `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 !Int
| LocateNucleusTodo
| LocateNucleus !Int
| LocateNucleusDone !(Complex r)
| LocateBondTodo
| LocateBond !Int
| 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)