
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.




link graphics
 
record circle(x,y,r)
global scale,xoffset,yoffset,yadjust
 
procedure main()
 
WOpen("size=400,400") | stop("Unable to open Window")
scale := 28
xoffset := WAttrib("width") / 2
yoffset := ( yadjust := WAttrib("height")) / 2
 
 
WC(c1 := circle(0,0,1),"black")
WC(c2 := circle(4,0,1),"black")
WC(c3 := circle(2,4,2),"black")
WC(c4 := Apollonius(c1,c2,c3,1,1,1),"green")    #/ Expects "Circle[x=2.00,y=2.10,r=3.90]" (green circle in image)
WC(c5 := Apollonius(c1,c2,c3,-1,-1,-1),"red")   #/ Expects "Circle[x=2.00,y=0.83,r=1.17]" (red circle in image)
 
 
WAttrib("fg=blue")
DrawLine( 0*scale+xoffset, yadjust-(-1*scale+yoffset),  0*scale+xoffset, yadjust-(4*scale+yoffset) )
DrawLine( -1*scale+xoffset, yadjust-(0*scale+yoffset),  4*scale+xoffset, yadjust-(0*scale+yoffset) )
WDone()
end
 
procedure WC(c,fg)  # write and plot circle
WAttrib("fg="||fg) 
DrawCircle(c.x*scale+xoffset, yadjust-(c.y*scale+yoffset), c.r*scale)
return write("Circle(x,y,r) := (",c.x,", ",c.y,", ",c.r,")")
end
 
procedure Apollonius(c1,c2,c3,s1,s2,s3)  # solve Apollonius
 
  v11 := 2.*(c2.x - c1.x)
  v12 := 2.*(c2.y - c1.y)
  v13 := c1.x^2 - c2.x^2 + c1.y^2 - c2.y^2 - c1.r^2 + c2.r^2
  v14 := 2.*(s2*c2.r - s1*c1.r)
 
  v21 := 2.*(c3.x - c2.x)
  v22 := 2.*(c3.y - c2.y)
  v23 := c2.x^2 - c3.x^2 + c2.y^2 - c3.y^2 - c2.r^2 + c3.r^2
  v24 := 2.*(s3*c3.r - s2*c2.r)
 
  w12 := v12/v11
  w13 := v13/v11
  w14 := v14/v11
 
  w22 := v22/v21-w12
  w23 := v23/v21-w13
  w24 := v24/v21-w14
 
  P := -w23/w22
  Q := w24/w22
  M := -w12*P-w13
  N := w14 - w12*Q
 
  a := N*N + Q*Q - 1
  b := 2*M*N - 2*N*c1.x + 2*P*Q - 2*Q*c1.y + 2*s1*c1.r
  c := c1.x*c1.x + M*M - 2*M*c1.x + P*P + c1.y*c1.y - 2*P*c1.y - c1.r*c1.r
 
  #// Find a root of a quadratic equation. This requires the circle centers not to be e.g. colinear
  D := b*b-4*a*c
  rs := (-b-sqrt(D))/(2*a)
  xs := M + N * rs
  ys := P + Q * rs
  return circle(xs,ys,rs)
end