gruff-0.2: src/MuAtom.hs

{-# LANGUAGE BangPatterns #-}
module MuAtom(MuAtom(..), MuProgress(..), muFromAddress, MuProgress'(..), muToAddress, MuProgress''(..), muLocate) where

import Control.Arrow ((***))
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)

import Number (R)

data MuAtom = MuAtom
  { muNucleus :: !(Complex R)
  , muSize    :: !Double
  , muOrient  :: !Double
  , muPeriod  :: !Integer
  }
  deriving (Read, Show, Eq)

data MuProgress
  = MuSplitTodo
  | MuSplitDone AngledInternalAddress [Angle]
  | MuAnglesTodo
  | MuAnglesDone !Rational !Rational
  | MuRayTodo
  | MuRay !Integer
  | MuRayDone !(Complex R)
  | MuNucleusTodo
  | MuNucleus !Integer
  | MuNucleusDone !(Complex R)
  | MuBondTodo
  | MuBond !Integer
  | MuBondDone !(Complex R)
  | MuSuccess !MuAtom
  | MuFailed
  deriving (Read, Show, Eq)

muFromAddress :: AngledInternalAddress -> [MuProgress]
muFromAddress addr = MuSplitTodo :
  let (!iaddr, !caddr) = splitAddress addr
      !p = addressPeriod iaddr
  in  MuSplitDone iaddr caddr : MuAnglesTodo : case externalAngles iaddr of
    Nothing -> [MuFailed]
    Just (!lo, !hi) -> MuAnglesDone lo hi : MuRayTodo :
      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 []  _ = [MuFailed]
          rgo [_] _ = [MuFailed]
          rgo ((i, (xl, xh)):m@((_, (yl, yh)):_)) f
            | i < fromIntegral sharpness * (p + 16) = MuRay i : rgo m f
            | dl + dh > dx + dy = MuRay i : rgo m f
            | otherwise = MuRayDone 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 -> MuNucleusTodo :
        let nuc = findNucleus p rayend
            nuc' = takeWhile (ok . snd) $ [ 1 .. ] `zip` nuc
            ngo []  _ = [MuFailed]
            ngo [_] _ = [MuFailed]
            ngo ((i, x):m@((_, y):_)) f
              | not (nearZero (x - y)) = MuNucleus i : ngo m f
              | otherwise = MuNucleusDone x : f x
        in  ngo nuc' $ \nucleus -> MuBondTodo :
          let bnd = findBond p nucleus 0.5
              bnd' = takeWhile (ok . snd) $ [ 1.. ] `zip` bnd
              bgo []  _ = [MuFailed]
              bgo [_] _ = [MuFailed]
              bgo ((i, x):m@((_, y):_)) f
                | not (nearZero (x - y)) = MuBond i : bgo m f
                | otherwise = MuBondDone 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 [MuSuccess atom] else [MuFailed]

data MuProgress'
  = MuCuspTodo
  | MuCuspDone !(Complex R)
  | MuDwellTodo
  | MuDwell !Integer
  | MuDwellDone !Integer
  | MuRayOutTodo
  | MuRayOut !Double
  | MuRayOutDone !(Complex R)
  | MuExternalTodo
  | MuExternalDone !Rational
  | MuAddressTodo
  | MuSuccess' AngledInternalAddress
  | MuFailed'
  deriving (Read, Show, Eq)

muToAddress :: MuAtom -> [MuProgress']
muToAddress mu = MuCuspTodo :
  let cusp = muNucleus mu - mkPolar (realToFrac (muSize mu)) (realToFrac ((muOrient mu)))
      er = 65536
      er2 = er * er
  in  MuCuspDone cusp : MuDwellTodo :
    let dgo z n f = MuDwell n : if magnitude2 z > er2 then f n else dgo (z * z + cusp) (n + 1) f
    in  dgo 0 0 $ \n -> MuDwellDone n : MuRayOutTodo :
      let rgo ((i,!_):izs@(_:_)) f = MuRayOut (fromIntegral i / (fromIntegral sharpness * fromIntegral n)) : rgo izs f
          rgo [(_,!z)] f | magnitude2 z > er2 = MuRayOutDone z : f z
          rgo _ _ = [MuFailed']
          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 -> MuExternalTodo :
        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  MuExternalDone angle : MuAddressTodo : case angledInternalAddress angle of
              Nothing -> [MuFailed']
              Just addr -> if addressPeriod addr /= muPeriod mu then [MuFailed'] else [MuSuccess' addr]

data MuProgress''
  = MuScanTodo
  | MuScan
  | MuScanDone !Integer
  | MuNucleusTodo'
  | MuNucleus' !Integer
  | MuNucleusDone' !(Complex R)
  | MuBondTodo'
  | MuBond' !Integer
  | MuBondDone' !(Complex R)
  | MuSuccess'' !MuAtom
  | MuFailed''
  deriving (Read, Show, Eq)

muLocate :: Complex R -> Double -> [MuProgress'']
muLocate c r = MuScanTodo : MuScan : case findPeriod 10000000 (realToFrac r) c of
  Nothing -> [MuFailed'']
  Just p -> MuScanDone p : MuNucleusTodo' :
    let ok w = magnitude2 w < 16 -- NaN -> False
        nuc = findNucleus p c
        nuc' = takeWhile (ok . snd) $ [ 1 .. ] `zip` nuc
        ngo []  _ = [MuFailed'']
        ngo [_] _ = [MuFailed'']
        ngo ((i, x):m@((_, y):_)) f
          | not (nearZero (x - y)) = MuNucleus' i : ngo m f
          | otherwise = MuNucleusDone' x : f x
    in  ngo nuc' $ \nucleus -> MuBondTodo' :
      let bnd = findBond p nucleus 0.5
          bnd' = takeWhile (ok . snd) $ [ 1.. ] `zip` bnd
          bgo []  _ = [MuFailed'']
          bgo [_] _ = [MuFailed'']
          bgo ((i, x):m@((_, y):_)) f
            | not (nearZero (x - y)) = MuBond' i : bgo m f
            | otherwise = MuBondDone' 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 [MuSuccess'' atom] else [MuFailed'']