gmndl-0.3: MuAtom.hs

{-# LANGUAGE BangPatterns, RecordWildCards #-}

{-

    gmndl -- Mandelbrot Set explorer
    Copyright (C) 2010  Claude Heiland-Allen <claudiusmaximus@goto10.org>

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License along
    with this program; if not, write to the Free Software Foundation, Inc.,
    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

-}

module MuAtom (muAtom) where

import Data.Complex (Complex((:+)), cis, magnitude, mkPolar, phase)
import Data.Ratio (denominator)
import Numeric.GSL.Root (RootMethod(Hybrids))
import qualified Numeric.GSL.Root as GSL

type N = Integer
type Q = Rational
type R = Double
type C = Complex R

-- a Mandelbrot Set mu-atom

data Atom = Atom{ nucleus :: !C, period :: !N, root :: !C, cardioid :: !Bool }
  deriving Show

continent :: Atom
continent = Atom 0 1 1 True

-- finding bond points

f :: N -> C -> C -> C
f !n !z !c
  | n == 0 = z
  | otherwise = let w = f (n - 1) z c in w * w + c

df :: N -> C -> C -> C
df !n !z !c = 2 ^ n * product [ f i z c | i <- [0 .. n - 1] ]

bondIter :: N -> C -> [R] -> [R]
bondIter !n !b [x0, x1, x2, x3] =
  let z = x0:+x1
      c = x2:+x3
      y0 :+ y1 =  f n z c - z
      y2 :+ y3 = df n z c - b
  in  [y0, y1, y2, y3]
bondIter _ _ _ = error "MuAtom.bondIter: internal error"

-- finding nucleus

l :: N -> C -> C
l !n !c
  | n == 0 = 0
  | otherwise = let z = l (n - 1) c in z * z + c

nucleusIter :: N -> [R] -> [R]
nucleusIter !n [x0, x1] =
  let c = x0:+x1
      y0 :+ y1 = l n c
  in  [y0, y1]
nucleusIter _ _ = error "MuAtom.nucleusIter: internal error"

-- finding descendants

muChild :: Atom -> Q -> Atom
muChild Atom{..} address =
  let -- some properties of the parent and its relation to the child
      size = magnitude (root - nucleus)
      angle = 2 * pi * realToFrac address
      bondAngle = cis angle
      child = period * denominator address
      solve = GSL.root Hybrids 1e-10 10000
      -- perturb from the stable nucleus to help ensure convergence to the bond point
      _initial@(ir :+ ii) = nucleus + mkPolar (size / 2) (phase (root - nucleus) + angle)
      ([_, _, br, bi], _) = solve (bondIter period bondAngle) [ ir, ii, ir, ii ]
      bondPoint = br :+ bi
      -- estimate where the nucleus will be
      radiusEstimate
        | cardioid  = size / m2 * sin (pi * realToFrac address)
        | otherwise = size / m2
        where m2 = fromIntegral (denominator address) ^ (2 :: N)
      deltaEstimate = bondPoint - nucleus
      _guess@(gr :+ gi) = bondPoint + mkPolar radiusEstimate (phase deltaEstimate)
      -- refine the nucleus estimate
      ([cr, ci], _) = solve (nucleusIter child) [gr, gi]
      childNucleus = cr :+ ci
  in  Atom childNucleus child bondPoint False

muChildren :: Atom -> [Q] -> [Atom]
muChildren a [] = [a]
muChildren a (q:qs) = let b = muChild a q in a : muChildren b qs

-- interface to the outside world

muAtom :: [Rational] -> (Double, Double, Double, Integer)
muAtom qs =
  let Atom{..} = last $ muChildren continent qs
      r :+ i = nucleus
      s = magnitude (nucleus - root)
      p = period
  in (r, i, s, p)