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)
data MuAtom r = MuAtom
{ muNucleus :: !(Complex r)
, muSize :: !Double
, muOrient :: !Double
, muPeriod :: !Int
}
deriving (Read, Show, Eq)
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
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)
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
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]
findAtom' :: [FindAtom r] -> Maybe (MuAtom r)
findAtom' ps = fromAtomSuccess =<< listToMaybe (filter isAtomSuccess ps)
findAtom_ :: (Floating r, NearZero r, Real r) => AngledInternalAddress -> Maybe (MuAtom r)
findAtom_ = findAtom' . findAtom
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
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]
findAddress' :: [FindAddress r] -> Maybe AngledInternalAddress
findAddress' ps = fromAddressSuccess =<< listToMaybe (filter isAddressSuccess ps)
findAddress_ :: (Floating r, NearZero r, Real r, RealFrac r) => MuAtom r -> Maybe AngledInternalAddress
findAddress_ = findAddress' . findAddress
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
locate :: (Floating r, NearZero r, Real r) => Complex r -> Double -> [Locate r]
locate c r = LocateScanTodo : LocateScan : case findPeriod 10000000 (realToFrac r) c of
Nothing -> [LocateFailed]
Just p -> LocateScanDone p : LocateNucleusTodo :
let ok w = magnitude2 w < 16
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
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]
locate' :: [Locate r] -> Maybe (MuAtom r)
locate' ps = fromLocateSuccess =<< listToMaybe (filter isLocateSuccess ps)
locate_ :: (Floating r, NearZero r, Real r) => Complex r -> Double -> Maybe (MuAtom r)
locate_ c r = locate' (locate c r)