{-# 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)