{-# LANGUAGE BangPatterns, RecordWildCards #-}
{-
gmndl -- Mandelbrot Set explorer
Copyright (C) 2010 Claude Heiland-Allen
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)