
Task

Implement a solution to the Problem of Apollonius   (description on Wikipedia)   which is the problem of finding the circle that is tangent to three specified circles   (colored black in the diagram below to the right).

There is an   algebraic solution   which is pretty straightforward.



The solutions to the example in the code are shown in the diagram   (below and to the right).

The red circle is "internally tangent" to all three black circles,   and the green circle is "externally tangent" to all three black circles.




 
Dim As String circle1 = "  0.000,  0.000,  1.000"
Dim As String circle2 = "  4.000,  0.000,  1.000"
Dim As String circle3 = "  2.000,  4.000,  2.000"
 
Sub ApolloniusSolver(c1 As String, c2 As String, c3 As String, s1 As Single, s2 As Single, s3 As Single)
    Dim As Single x1, x2, x3, y1, y2, y3, r1, r2, r3
    Dim As Single v11, v12, v13, v14, v21, v22, v23, v24, w12, w13, w14
    Dim As Single w22, w23, w24,P, Q, M, N, a, b, c, D
    Dim As Single Radius, XPos, YPos
 
    x1 = Val(Mid(c1, 3, 1)): y1 = Val(Mid(c1, 11, 1)): r1 = Val(Mid(c1, 19, 1))
    x2 = Val(Mid(c2, 3, 1)): y2 = Val(Mid(c2, 11, 1)): r2 = Val(Mid(c2, 19, 1))
    x3 = Val(Mid(c3, 3, 1)): y3 = Val(Mid(c3, 11, 1)): r3 = Val(Mid(c3, 19, 1))
 
    v11 = 2 * x2 - 2 * x1
    v12 = 2 * y2 - 2* y1
    v13 = x1 * x1 - x2 * x2 + y1 * y1 - y2 * y2 - r1 * r1 + r2 * r2
    v14 = 2 * s2 * r2 - 2 * s1 * r1
 
    v21 = 2 * x3 - 2 * x2
    v22 = 2 * y3 - 2 * y2
    v23 = x2 * x2 - x3 * x3 + y2 * y2 - y3 * y3 - r2 * r2 + r3 * r3
    v24 = 2 * s3 * r3 - 2 * s2 * r2
 
    w12 = v12 / v11
    w13 = v13 / v11
    w14 = v14 / v11
 
    w22 = v22 / v21 - w12
    w23 = v23 / v21 - w13
    w24 = v24 / v21 - w14
 
    P   = 0 - w23 / w22
    Q   =     w24 / w22
    M   = 0 - w12 * P - w13
    N   =     w14 - w12 *  Q
 
    a   = N * N + Q * Q - 1
    b   = 2 * M * N - 2 * N * x1 + 2 * P * Q - 2 * Q * y1 + 2 * s1 * r1
    c   = x1 * x1 + M * M -2 * M * x1 + P * P + y1 * y1 - 2 * P * y1 - r1 * r1
 
    D   = b *  b - 4 * a * c
 
    Radius = (0 - b - Sqr(D)) / (2 * a)
    XPos   = M + N * Radius
    YPos   = P + Q * Radius
 
    Print Using " ##.###, ##.###, ##.###"; XPos; YPos; Radius
End Sub
 
Print "  x_pos   y_pos   radius"
Print circle1
Print circle2
Print circle3
Print
Print "R1: " : ApolloniusSolver(circle1, circle2, circle3,  1,  1,  1)
Print "R2: " : ApolloniusSolver(circle1, circle2, circle3, -1, -1, -1)
Sleep
 