{-# LANGUAGE BangPatterns, DeriveDataTypeable #-}
{- |
Module      :  Fractal.RUFF.Mandelbrot.Ray
Copyright   :  (c) Claude Heiland-Allen 2011
License     :  BSD3

Maintainer  :  claudiusmaximus@goto10.org
Stability   :  unstable
Portability :  portable

External angles define external rays which can be traced back from
the circle at infinity to the boundary of the Mandelbrot Set.

The algorithm is based on Tomoki Kawahira's paper
/An algorithm to draw external rays of the Mandelbrot set/
<http://www.math.nagoya-u.ac.jp/~kawahira/programs/mandel-exray.pdf>.
-}

module Fractal.RUFF.Mandelbrot.Ray (externalRay) where

import Fractal.RUFF.Types.Complex (Complex(), magnitude, mkPolar)
import Fractal.RUFF.Mandelbrot.Address (Angle, double)

-- | Compute the external ray for an external angle with a given
--   accuracy, sharpness and starting radius.  For example:
--
-- > externalRay 1e-10 8 (2**24) (1/3)
--
externalRay :: (Ord r, Floating r) => r -> Int -> r -> Angle -> [Complex r]
externalRay epsilon sharpness radius angle = map fst . iterate step $ (mkPolar radius (2 * pi * fromRational angle), (0, 0))
  where
    -- step :: (NearZero r, Floating r) => (Complex r, (Int, Int)) -> (Complex r, (Int, Int))
    step (!c, (!k0, !j0))
      | j > sharpness = step (c, (k0 + 1, 0))
      | otherwise = (n c, (k0, j0 + 1))
      where
        k = k0 + 1
        j = j0 + 1
        m = (k - 1) * sharpness + j
        r = radius ** ((1/2) ** (fromIntegral m / fromIntegral sharpness))
        t = mkPolar (r ** (2 ** fromIntegral k0)) (2 * pi * fromRational (iterate double angle !! k0))
        n !z = let d = (cc - t) / dd in if not (magnitude d > epsilon) then z else n (z - d)
          where
            (cc, dd) = ncnd k
            ncnd 1 = (z, 1)
            ncnd i = let (!nc, !nd) = ncnd (i - 1) in (nc * nc + z, 2 * nc * nd + 1)