Apollonian Gasket

A very interesting and quite pretty mathematical object for those who like such things. There are many other ways to colour the circles...

Now that I've cracked the formula I might do the full circle version which is even more pleasing to the eye.

' A version of the 'Apollonian Gasket'
'
' 20060418 Paul Fuller
'
' Use as you please
'
' References:
'
' - astronomy.swin.edu.au/~pbourke/fractals/apollony/	General description of the Apollonian Gasket
' - mathworld.wolfram.com/SoddyCircles.html		Radii of inner and outer Soddy circles
' - sonoma.edu/users/w/wilsonst/papers/geometry/circles/	Intersection of two circles
'

SuperStrict

Type T_Circle

	Field X:Float
	Field Y:Float
	Field R:Float
	Field CR:Int
	Field CG:Int
	Field CB:Int
	
	Function CreateCircleXY:T_Circle(x:Float, y:Float)
	
		Return CreateCircle(x, y, 0)
		
	End Function
	
	Function CreateCircle:T_Circle(x:Float, y:Float, r:Float)
		Local c:T_Circle = New T_Circle
		
		c.X = x
		c.Y = y
		c.R = r
		
		c.CR = 50
		c.CG = 50
		c.CB = 50
	
		Return c
	End Function
	
	Method Draw()
		
		SetColor(CR, CG, CB)
		DrawOval(X - R - MidX, Y - R- MidY, R * 2, R * 2)
		
	End Method

	Method SetRGB(r:Int, g:Int, b:Int)
	
		CR = r
		CG = g
		CB = b
		
	End Method
		
End Type

' Global offset to position the interesting part in the centre of the screen
Global MidX:Float
Global MidY:Float

' Mainline
Initialise()

Local C1:T_Circle = T_Circle.CreateCircleXY(0, 0)		' Controlled via the Mouse
Local C2:T_Circle = T_Circle.CreateCircleXY(0, 0)
Local C3:T_Circle = T_Circle.CreateCircleXY(1280, -500)

C1.SetRGB(0, 255, 0)
C2.SetRGB(0, 0, 255)
C3.SetRGB(255, 0, 0)

While Not KeyHit(KEY_ESCAPE)
	
	Cls
	
	C1.X = MouseX() + MidX
	C1.Y = MouseY() + MidY
	
	Adjust(C1, C2, C3)
	Divide(1, C1, C2, C3)
	
	C1.Draw()
	C2.Draw()
	C3.Draw()
	
	Flip True
Wend

End

Function Initialise()

	Graphics 1280, 1024, 16

	Cls
	
End Function

' Adjust the radii of three circles so that they are pairwaise tangent
Function Adjust(A:T_Circle, B:T_Circle, C:T_Circle)

	Local dA:Float = Distance(B.X, B.Y, C.X, C.Y)
	Local dB:Float = Distance(A.X, A.Y, C.X, C.Y)
	Local dC:Float = Distance(A.X, A.Y, B.X, B.Y)
	
	A.R = (dB + dC - dA) / 2
	B.R = (dA + dC - dB) / 2
	C.R = (dA + dB - dC) / 2
	
End Function

' Recursively divide three pairwise tangent circles by adding the inner Soddy circle
Function Divide(level:Int, C1:T_Circle, C2:T_Circle, C3:T_Circle)

	' Take the radii of kissing circles
	Local r1:Float = C1.R
	Local r2:Float = C2.R
	Local r3:Float = C3.R
	
	' Calculate radii of inner and outer Soddy circles (outer is not currently used)
	Local top:Float = r1 * r2 * r3
	Local part1:Float = r1 * r2 + r1 * r3 + r2 * r3
	Local part2:Float = 2 * Sqr(top * (r1 + r2 + r3))
	 
	Local innerR:Float = top / (part1 + part2)
	Local outerR:Float = top / (part1 - part2)
	
	' Bail-out if circle is too small to show
	If innerR < 1
		Return
	EndIf
	
	' Calculate coordinates of inner Soddy circle
	'
	' - Calculate the two points of intersection of C2 and C3 extended by the radius of the inner Soddy circle
	' - See which point is closer to C1
	'
	
	Local x1:Float = C2.X
	Local y1:Float = C2.Y
	Local a:Float = C2.R + innerR
	
	Local x2:Float = C3.X
	Local y2:Float = C3.Y
	Local b:Float = C3.R + innerR
	
	Local dx:Float = x2 - x1
	Local dy:Float = y2 - y1
	Local ds:Float = dx * dx + dy * dy	' Distance squared
	
	Local term1:Float = Sqr(((a + b) * (a + b) - ds) * (ds - (b - a) * (b - a))) / (2 * ds)
	Local term2:Float = (a * a - b * b) / (2 * ds)
	
	' Mid-point
	Local mx:Float = (x2 + x1) / 2 + (x2 - x1) * term2
	Local my:Float = (y2 + y1) / 2 + (y2 - y1) * term2
	
	' Offset from the mid-point
	Local ox:Float = (y2 - y1) * term1
	Local oy:Float = (x2 - x1) * term1
	
	Local d1:Float = Distance(mx + ox, my - oy, C1.X, C1.Y)
	Local d2:Float = Distance(mx - ox, my + oy, C1.X, C1.Y)
	
	' See which is closer and define that as the inner Soddy circle centre
	Local innerS:T_Circle
	If d1 <= d2
		innerS = T_Circle.CreateCircle(mx + ox, my - oy, innerR)
	Else
		innerS = T_Circle.CreateCircle(mx - ox, my + oy, innerR)
	EndIf

	' Take the first inner Soddy circle centre as the centre for drawing everything else
	If level = 1
		MidX = innerS.X - 640
		MidY = innerS.Y - 512
	EndIf

	' Determine a new colour by combining the three parent circle colours in some way		
	Local f:Float = 0.95 ^ (level - 1)
	Local cR:Int = 255 - (C1.CR + C2.CR + C3.CR) / 3 * f
	Local cG:Int = 255 - (C1.CG + C2.CG + C3.CG) / 3 * f
	Local cB:Int = 255 - (C1.CB + C2.CB + C3.CB) / 3 * f

	innerS.SetRGB(cR, cG, cB)
	innerS.Draw()
	
	' Recurse
	Divide(level + 1, C1, C2, innerS)
	Divide(level + 1, C1, C3, innerS)
	Divide(level + 1, C2, C3, innerS)
	
End Function

' Return the distance between (x1, y1) and (x2, y2)
Function Distance:Float(x1:Float, y1:Float, x2:Float, y2:Float)
	Local xd:Float = x1 - x2
	Local yd:Float = y1 - y2
	
	Return Sqr(xd * xd + yd * yd)
End Function