Following on from the previous post in this series, today’s post completes the implementation to create a full Apollonian gasket in AutoCAD using F#.
As a comment on the original Common LISP implementation, someone had contributed a more complete version which allowed me to complete today’s F# version.
Here’s the additional F# file for the project (which I’ll be providing in full at the end of the series):
module CirclePackingFullFs
open System.Numerics;
// Use Descartes' theorem to calculate the radius/position
// of the 4th circle
// k4 = k1 + k2 + k3 +/- sqrt(k1k2 + k2k3 + k3k1)
let solve a b c =
[(+);(-)] |>
List.map
(fun f -> (f (a + b + c) (2. * (sqrt (a*b + b*c + a*c)))))
let multComplex a (c : Complex) =
Complex(a * c.Real, a * c.Imaginary)
let solveComplex (a : Complex) (b : Complex) (c : Complex) =
[(+);(-)] |>
List.map
(fun f ->
(f (a + b + c)
(multComplex 2.
(Complex.Sqrt
(Complex.Multiply(a,b) +
Complex.Multiply(b,c) +
Complex.Multiply(a,c))))))
// Return a list of all solution circles (positive and negative)
let apollonianCircle (x1,y1,k1) (x2,y2,k2) (x3,y3,k3) =
let p1 = Complex(x1,y1)
let p2 = Complex(x2,y2)
let p3 = Complex(x3,y3)
let kz1 = multComplex k1 p1
let kz2 = multComplex k2 p2
let kz3 = multComplex k3 p3
let curSols = solve k1 k2 k3
let posSols = solveComplex kz1 kz2 kz3
List.map2
(fun cur pos ->
let r4 = 1. / cur
let pos4 = multComplex r4 pos
pos4.Real, pos4.Imaginary, cur)
curSols posSols
// Only return certain circles - for use with "hard" gaps
let apollonianCircle2 (x1,y1,k1) (x2,y2,k2) (x3,y3,k3) =
let p1 = Complex(x1,y1)
let p2 = Complex(x2,y2)
let p3 = Complex(x3,y3)
let kz1 = multComplex k1 p1
let kz2 = multComplex k2 p2
let kz3 = multComplex k3 p3
let cur = List.head (solve k1 k2 k3)
let pos = List.nth (solveComplex kz1 kz2 kz3) 1
let r4 = 1. / cur
let pos4 = multComplex r4 pos
pos4.Real, pos4.Imaginary, cur
let rec easyGap c1 c2 c3 steps =
match steps with
| 0 -> []
| _ ->
let c4 = List.head(apollonianCircle c1 c2 c3)
[(c4,steps-1)] @
easyGap c1 c2 c4 (steps-1) @
easyGap c2 c3 c4 (steps-1) @
easyGap c3 c1 c4 (steps-1)
let hardGap c1 c2 c3 steps =
match steps with
| 0 -> []
| _ ->
let g1 = apollonianCircle2 c1 c2 c3
[(g1,steps-1)] @
easyGap c2 c3 g1 (steps-1) @
easyGap c2 c1 g1 (steps-1) @
if steps = 1 then
[]
else
let g2 = apollonianCircle2 c1 c3 g1
[(g2,steps-2)] @
easyGap c3 g1 g2 (steps-2) @
easyGap c1 g1 g2 (steps-2) @
if steps = 2 then
[]
else
let g3 = apollonianCircle2 c1 c3 g2
[(g3,steps-3)] @
easyGap c1 g2 g3 (steps-3) @
easyGap c1 c3 g3 (steps-3) @
if steps = 3 then
[]
else
let g4 = apollonianCircle2 c3 g2 g3
[(g4,steps-4)] @
easyGap c3 g2 g4 (steps-4) @
easyGap g2 g3 g4 (steps-4)
type Packer() =
static member ApollonianGasket outerRad steps =
let rt3 = sqrt(3.)
let size = outerRad * 2.
let a = 1. + 2. / rt3
let innerRad = outerRad / a
let innerCur = 1. / innerRad
let h = innerRad * rt3
let c0 = outerRad, outerRad, -1. / outerRad
let c1 = outerRad, size - innerRad, innerCur
let c2 = outerRad + innerRad, size - (h + innerRad), innerCur
let c3 = outerRad - innerRad, size - (h + innerRad), innerCur
let c4 = List.head(apollonianCircle c1 c2 c3)
let nextStep = steps-1
// Main circles
// (Assumption is that the outer circle already exists)
[(c1,steps);(c2,steps);(c3,steps);(c4,nextStep)] @
// Inner gaps
easyGap c1 c2 c4 nextStep @
easyGap c2 c3 c4 nextStep @
easyGap c3 c1 c4 nextStep @
// Peripheral gaps
easyGap c1 c2 c0 nextStep @
easyGap c2 c3 c0 nextStep @
hardGap c3 c1 c0 nextStep
The hardGap function does a fair amount of manual filling of the top-left gap… I suspect this could also have been addressed by mirroring or rotating the results from one of the easyGap calls, but there you go. Ultimately the code performs well enough.
The corresponding C# loader file clearly also needs updating to define an additional command calling the F# code:
using Autodesk.AutoCAD.ApplicationServices;
using Autodesk.AutoCAD.DatabaseServices;
using Autodesk.AutoCAD.EditorInput;
using Autodesk.AutoCAD.Geometry;
using Autodesk.AutoCAD.Runtime;
namespace CirclePackingLoader
{
public class Commands
{
[CommandMethod("AGFC")]
public static void ApollonianGasket()
{
PackApollonianCircles(true);
}
[CommandMethod("AGFCI")]
public static void InnerApollonianGasket()
{
PackApollonianCircles(false);
}
public static void PackApollonianCircles(bool full)
{
Document doc =
Application.DocumentManager.MdiActiveDocument;
Database db = doc.Database;
Editor ed = doc.Editor;
PromptEntityOptions peo =
new PromptEntityOptions("\nSelect circle: ");
peo.AllowNone = true;
peo.SetRejectMessage("\nMust be a circle.");
peo.AddAllowedClass(typeof(Circle), false);
PromptEntityResult per = ed.GetEntity(peo);
if (per.Status != PromptStatus.OK &&
per.Status != PromptStatus.None
)
return;
Transaction tr =
doc.TransactionManager.StartTransaction();
using (tr)
{
BlockTableRecord btr =
(BlockTableRecord)tr.GetObject(
db.CurrentSpaceId, OpenMode.ForWrite
);
// If the user selected a circle, use it
Circle cir;
if (per.Status != PromptStatus.None && per.ObjectId != null)
{
cir =
tr.GetObject(per.ObjectId, OpenMode.ForRead) as Circle;
}
else
{
// Otherwise create a new one at the default location
cir = new Circle(Point3d.Origin, Vector3d.ZAxis, 500);
btr.AppendEntity(cir);
tr.AddNewlyCreatedDBObject(cir, true);
}
int numSteps = 11;
var res =
full ?
CirclePackingFullFs.Packer.ApollonianGasket(
cir.Radius, numSteps
)
:
CirclePackingFs.Packer.ApollonianGasket(
cir.Radius, numSteps
);
Vector3d offset =
cir.Center - new Point3d(cir.Radius, cir.Radius, 0.0);
if (res.Length > 0)
{
// We get an F# list of tuples containing our circle
// definitions (no need for a special class)
foreach (var tup in res)
{
// Our circles are defined in terms of position (x,y)
// and curvature (the 3rd item in the nested tuple)
double curvature = System.Math.Abs(tup.Item1.Item3);
if (1.0 / curvature > 0.0)
{
// x and y are Items 1 and 2 (in the nested tuple)
// respectively
Circle c =
new Circle(
new Point3d(
tup.Item1.Item1, tup.Item1.Item2, 0.0
) + offset,
Vector3d.ZAxis,
1.0 / curvature
);
// Color index will be based on the "level" of each
// circle (item 2 in the top-level tuple)
c.ColorIndex = (numSteps - tup.Item2) + 1;
btr.AppendEntity(c);
tr.AddNewlyCreatedDBObject(c, true);
}
}
}
tr.Commit();
}
}
}
}
When we run our new AGFC command, we now see a full Apollonian gasket generated:
Over the next few posts we’ll make the transition across to 3D, looking at sphere packing.
Subscribe via RSS