module Fractal.RUFF.Mandelbrot.Ray (externalRay) where
import Fractal.RUFF.Types.Complex (Complex(), magnitude, mkPolar)
import Fractal.RUFF.Mandelbrot.Address (Angle, double)
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 (!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)