// F# adaptation of Kean Walmsley's Apollonian gasket code - 2012-02-06
// http://through-the-interface.typepad.com/through_the_interface/2012/02/circle-packing-in-autocad-creating-an-apollonian-gasket-using-net.html
// with Gaston Nunez's comment (#1) taken into account, see also
// http://en.wikibooks.org/wiki/Fractals/Apollonian_fractals#Soddy.27s_algorithm
module CirclePacking
open Autodesk.AutoCAD.EditorInput
open Autodesk.AutoCAD.DatabaseServices
open Autodesk.AutoCAD.Geometry
open Autodesk.AutoCAD.Runtime
type acadApp = Autodesk.AutoCAD.ApplicationServices.Application
type MyCircle = {
Center : Point2d
Radius : float
ColorIndex : int } with
static member Create(x, y, radius) = {
Center = new Point2d(x, y)
Radius = radius
ColorIndex = 5 }
type TangentType = OutsideAll | OutsideTwo | Outer
let findApollonianCircle c1 c2 c3 tanType =
let (k1, k2, k3) =
match tanType with
| OutsideAll -> 1. / c1.Radius, 1. / c2.Radius, 1. / c3.Radius
| OutsideTwo -> 1. / c1.Radius, 1. / c2.Radius, -1. / c3.Radius
| Outer -> -1. / c1.Radius, -1. / c2.Radius, -1. / c3.Radius
let k4 = k1 + k2 + k3 + 2. * sqrt(k1 * k2 + k1 * k3 + k2 * k3)
let z1 = Complex.Create(c1.Center.X, c1.Center.Y)
let z2 = Complex.Create(c2.Center.X, c2.Center.Y)
let z3 = Complex.Create(c3.Center.X, c3.Center.Y)
let z4 = k1 * z1 + k2 * z2 + k3 * z3
let z4 = z4 + 2. * sqrt(k1 * k2 * z1 * z2 + k1 * k3 * z1 * z3 + k2 * k3 * z2 * z3)
let z4 = z4 / Complex.Create(k4, 0.)
MyCircle.Create(z4.r, z4.i, 1. / k4)
// Draw a circle tangent to these three circles and that is
// outside all three
let rec findCircleOutsideAll level cir0 cir1 cir2 = seq{
let newCir = findApollonianCircle cir0 cir1 cir2 OutsideAll
yield newCir
if level > 1 then
yield! findCircleOutsideAll (level - 1) cir0 cir1 newCir
yield! findCircleOutsideAll (level - 1) cir0 cir2 newCir
yield! findCircleOutsideAll (level - 1) cir1 cir2 newCir }
// Draw a circle tangent to these three circles and that is
// outside two of them
let rec findCircleOutsideTwo level cir0 cir1 contains = seq{
let newCir = findApollonianCircle cir0 cir1 contains OutsideTwo
yield newCir
if level > 1 then
yield! findCircleOutsideTwo (level - 1) newCir cir0 contains
yield! findCircleOutsideTwo (level - 1) newCir cir1 contains
yield! findCircleOutsideAll (level - 1) cir0 cir1 newCir }
// Find the Apollonian packing
let findApollonianPacking level =
// Create the three central tangent circles
let radius = 1. / (2. + 4. / sqrt 3.)
// The centroid of the central circles lies at point (1, 1)
let (x0, y0) = 1.0, 0.5 + radius
let (x1, y1) = x0 - radius, y0 + radius * sqrt 3.
let (x2, y2) = x0 + radius, y1
let cir0 = MyCircle.Create(x0, y0, radius)
let cir1 = MyCircle.Create(x1, y1, radius)
let cir2 = MyCircle.Create(x2, y2, radius)
// Find the circle that contains them all, we know already
// its center (1, 1) and diameter (1)
let outer = findApollonianCircle cir0 cir1 cir2 Outer
seq{
yield outer
yield cir0
yield cir1
yield cir2
if level > 1 then
// Find the central circle
yield! findCircleOutsideAll (level - 1) cir0 cir1 cir2
// Find circles tangent to the big circle
yield! findCircleOutsideTwo (level - 1) cir0 cir1 outer
yield! findCircleOutsideTwo (level - 1) cir1 cir2 outer
yield! findCircleOutsideTwo (level - 1) cir2 cir0 outer }
let drawCircles cirId level =
let db = HostApplicationServices.WorkingDatabase
use tr = db.TransactionManager.StartTransaction()
let btr =
tr.GetObject(
db.CurrentSpaceId, OpenMode.ForWrite )
:?> BlockTableRecord
let cir =
if cirId = ObjectId.Null then
// Create a new circlee at the default location
let cir =
new Circle(Point3d.Origin, Vector3d.ZAxis, 500.)
btr.AppendEntity cir |> ignore
tr.AddNewlyCreatedDBObject(cir, true)
cir
else
// Otherwise use the selected circle
tr.GetObject(cirId, OpenMode.ForRead) :?> Circle
let mat =
Matrix3d.Displacement(cir.Center.GetAsVector()) *
Matrix3d.Scaling(cir.Diameter, Point3d.Origin) *
Matrix3d.WorldToPlane(cir.Normal).Inverse() *
Matrix3d.Displacement(new Vector3d(-1., -1., 0.))
// Add each of the returned circles to the current space
let cnt =
findApollonianPacking level
let c =
new Circle(
new Point3d(myC.Center.X, myC.Center.Y, 0.),
Vector3d.ZAxis, myC.Radius,
ColorIndex = myC.ColorIndex )
c.TransformBy mat
btr.AppendEntity c |> ignore
tr.AddNewlyCreatedDBObject(c, true)
n + 1 ) 0
tr.Commit()
cnt
[<CommandMethod "myAG">]
let apollonianGasket() =
let ed = acadApp.DocumentManager.MdiActiveDocument.Editor
let peo =
new PromptEntityOptions(
"\nSelect circle: ",
AllowNone = true )
peo.SetRejectMessage "Must be a circle."
peo.AddAllowedClass(typeof<Circle>, false)
let per = ed.GetEntity peo
if per.Status = PromptStatus.OK ||
per.Status = PromptStatus.None then
let pio =
new PromptIntegerOptions(
"\nEnter maximum level: ",
AllowNone = true,
LowerLimit = 1, // 2*3^0+2 = 4 circles
UpperLimit = 13, // 2*3^12+2 = 1062884 circles
DefaultValue = 8 )
let pir = ed.GetInteger pio
if pir.Status = PromptStatus.OK ||
pir.Status = PromptStatus.None then
drawCircles
( if per.Status = PromptStatus.OK then per.ObjectId
else ObjectId.Null )
( if pir.Status = PromptStatus.OK then pir.Value
else pio.DefaultValue )
|> fun cnt -> ed.WriteMessage("\n{0} Circles. ", cnt)