{-# LANGUAGE Rank2Types
           , FlexibleContexts
           , ViewPatterns
  #-}
{-# OPTIONS_GHC -fno-warn-missing-methods #-}

-----------------------------------------------------------------------------
-- |
-- Module      :  Diagrams.TwoD.Apollonian
-- Copyright   :  (c) 2011 Brent Yorgey
-- License     :  BSD-style (see LICENSE)
-- Maintainer  :  byorgey@cis.upenn.edu
--
-- Generation of Apollonian gaskets.  Any three mutually tangent
-- circles uniquely determine exactly two others which are mutually
-- tangent to all three.  This process can be repeated, generating a
-- fractal circle packing.
--
-- See J. Lagarias, C. Mallows, and A. Wilks, \"Beyond the Descartes
-- circle theorem\", /Amer. Math. Monthly/ 109 (2002), 338--361.
-- <http://arxiv.org/abs/math/0101066>.
--
-----------------------------------------------------------------------------

module Diagrams.TwoD.Apollonian
       ( -- * Circles

         Circle(..), mkCircle, center, radius

         -- * Descartes' Theorem
       , descartes, other, initialConfig

         -- * Apollonian gasket generation

       , apollonian

         -- * Diagram generation

       , drawCircle
       , drawGasket
       , apollonianGasket

       ) where

import Data.Complex
import Data.Foldable (foldMap)

import Diagrams.Prelude hiding (radius, center)

import Control.Arrow (second)

------------------------------------------------------------
--  Circles
------------------------------------------------------------

-- | Representation for circles that lets us quickly compute an
--   Apollonian gasket.
data Circle = Circle { bend :: Double
                       -- ^ The bend is the reciprocal of signed
                       --   radius: a negative radius means the
                       --   outside and inside of the circle are
                       --   switched.  The bends of any four mutually
                       --   tangent circles satisfy Descartes'
                       --   Theorem.
                     , cb   :: Complex Double
                       -- ^ /Product/ of bend and center represented
                       --   as a complex number.  Amazingly, these
                       --   products also satisfy the equation of
                       --   Descartes' Theorem.
                     }
  deriving (Eq, Show)

-- | Create a @Circle@ given a signed radius and a location for its center.
mkCircle :: Double -- ^ signed radius
         -> P2     -- ^ center
         -> Circle
mkCircle r (unp2 -> (x,y)) = Circle (1/r) (b*x :+ b*y)
  where b = 1/r

-- | Get the center of a circle.
center :: Circle -> P2
center (Circle b (cbx :+ cby)) = p2 (cbx / b, cby / b)

-- | Get the (unsigned) radius of a circle.
radius :: Circle -> Double
radius = abs . recip . bend

liftF :: (forall a. Floating a => a -> a) -> Circle -> Circle
liftF f (Circle b c) = Circle (f b) (f c)

liftF2 :: (forall a. Floating a => a -> a -> a) -> Circle -> Circle -> Circle
liftF2 f (Circle b1 cb1) (Circle b2 cb2) = Circle (f b1 b2) (f cb1 cb2)

instance Num Circle where
  (+) = liftF2 (+)
  (-) = liftF2 (-)
  (*) = liftF2 (*)
  negate = liftF negate
  abs = liftF abs
  fromInteger n = Circle (fromInteger n) (fromInteger n)

instance Fractional Circle where
  (/) = liftF2 (/)
  recip = liftF recip

-- | The @Num@, @Fractional@, and @Floating@ instances for @Circle@
--   (all simply lifted elementwise over @Circle@'s fields) let us use
--   Descartes' Theorem directly on circles.
instance Floating Circle where
  sqrt = liftF sqrt

------------------------------------------------------------
--  Descartes' Theorem
------------------------------------------------------------

-- | Descartes' Theorem states that if @b1@, @b2@, @b3@ and @b4@ are
--   the bends of four mutually tangent circles, then
--
--   @
--     b1^2 + b2^2 + b3^2 + b4^2 = 1/2 * (b1 + b2 + b3 + b4)^2.
--   @
--
--   Surprisingly, if we replace each of the @bi@ with the /product/
--   of @bi@ and the center of the corresponding circle (represented
--   as a complex number), the equation continues to hold! (See the
--   paper referenced at the top of the module.)
--
--   @descartes [b1,b2,b3]@ solves for @b4@, returning both solutions.
--   Notably, @descartes@ works for any instance of @Floating@, which
--   includes both @Double@ (for bends), @Complex Double@ (for
--   bend/center product), and @Circle@ (for both at once).
descartes :: Floating a => [a] -> [a]
descartes [b1,b2,b3] = [r + s, -r + s]
  where r = 2 * sqrt (b1*b2 + b1*b3 + b2*b3)
        s = b1+b2+b3
descartes _ = error "descartes must be called on a list of length 3"

-- | If we have /four/ mutually tangent circles we can choose one of
--   them to replace; the remaining three determine exactly one other
--   circle which is mutually tangent.  However, in this situation
--   there is no need to apply 'descartes' again, since the two
--   solutions @b4@ and @b4'@ satisfy
--
--   @
--     b4 + b4' = 2 * (b1 + b2 + b3)
--   @
--
--   Hence, to replace @b4@ with its dual, we need only sum the other
--   three, multiply by two, and subtract @b4@.  Again, this works for
--   bends as well as bend/center products.
other :: Num a => [a] -> a -> a
other xs x = 2 * sum xs - x

-- | Generate an initial configuration of four mutually tangent
--   circles, given just the signed bends of three of them.
initialConfig :: Double -> Double -> Double -> [Circle]
initialConfig b1 b2 b3 = cs ++ [c4]
  where cs     = [Circle b1 0, Circle b2 ((b2/b1 + 1) :+ 0), Circle b3 cb3]
        a      = 1/b1 + 1/b2
        b      = 1/b1 + 1/b3
        c      = 1/b2 + 1/b3
        x      = (b*b + a*a - c*c)/(2*a)
        y      = sqrt (b*b - x*x)
        cb3    = b3*x :+ b3*y
        [c4,_] = descartes cs

------------------------------------------------------------
--  Gasket generation
------------------------------------------------------------

select :: [a] -> [(a, [a])]
select [] = []
select (x:xs) = (x,xs) : (map . second) (x:) (select xs)

-- | Given a threshold radius and a list of /four/ mutually tangent
--   circles, generate the Apollonian gasket containing those circles.
--   Stop the recursion when encountering a circle with an (unsigned)
--   radius smaller than the threshold.
apollonian :: Double -> [Circle] -> [Circle]
apollonian thresh cs
  =  cs
  ++ (concat . map (\(c,cs') -> apollonian' thresh (other cs' c) cs') . select $ cs)

apollonian' :: Double -> Circle -> [Circle] -> [Circle]
apollonian' thresh cur others
  | radius cur < thresh = []
  | otherwise = cur
              : (concat $
                   map (\(c, cs') -> apollonian' thresh
                                       (other (cur:cs') c)
                                       (cur:cs')
                       )
                       (select others)
                )

------------------------------------------------------------
--  Diagram generation
------------------------------------------------------------

-- | Draw a circle.
drawCircle :: (Renderable (Path R2) b) => Double -> Circle -> Diagram b R2
drawCircle w c = circle (radius c) # moveTo (center c)
                                   # lw w # fcA transparent

-- | Draw a generated gasket, using a line width 0.003 times the
--   radius of the largest circle.
drawGasket :: (Renderable (Path R2) b) => [Circle] -> Diagram b R2
drawGasket cs = foldMap (drawCircle w) cs
  where w = (*0.003) . maximum . map radius $ cs

-- | Draw an Apollonian gasket: the first argument is the threshold;
--   the recursion will stop upon reaching circles with radii less than
--   it. The next three arguments are bends of three circles.
apollonianGasket :: (Renderable (Path R2) b)
                 => Double -> Double -> Double -> Double -> Diagram b R2
apollonianGasket thresh b1 b2 b3 = drawGasket . apollonian thresh $ (initialConfig b1 b2 b3)