Problem of Apollonius

include "NSLog.incl"

begin record Circle
CGPoint     center
double      radius
CFStringRef locator
end record

local fn CircleToString( c as Circle ) as CFStringRef
end fn = fn StringWithFormat( @"%@ Circle( x = %0.3f, y = %0.3f, radius = %0.3f )", c.locator, c.center.x, c.center.y, c.radius )

local fn SolveApollonius( c1 as Circle, c2 as Circle, c3 as Circle, s1 as Double, s2 as Double, s3 as Double ) as Circle
'~'1
Circle result

double x1 = c1.center.x
double y1 = c1.center.y
double r1 = c1.radius

double x2 = c2.center.x
double y2 = c2.center.y
double r2 = c2.radius

double x3 = c3.center.x
double y3 = c3.center.y
double r3 = c3.radius

double v11 = 2*x2 - 2*x1
double v12 = 2*y2 - 2*y1
double v13 = x1*x1 - x2*x2 + y1*y1 - y2*y2 - r1*r1 + r2*r2
double v14 = 2*s2*r2 - 2*s1*r1

double v21 = 2*x3 - 2*x2
double v22 = 2*y3 - 2*y2
double v23 = x2*x2 - x3*x3 + y2*y2 - y3*y3 - r2*r2 + r3*r3
double v24 = 2*s3*r3 - 2*s2*r2

double w12 = v12/v11
double w13 = v13/v11
double w14 = v14/v11

double w22 = v22/v21-w12
double w23 = v23/v21-w13
double w24 = v24/v21-w14

double P = -w23/w22
double Q = w24/w22
double M = -w12*P-w13
double N = w14 - w12*Q

double a = N*N + Q*Q - 1
double b = 2*M*N - 2*N*x1 + 2*P*Q - 2*Q*y1 + 2*s1*r1
double c = x1*x1 + M*M - 2*M*x1 + P*P + y1*y1 - 2*P*y1 - r1*r1

double D = b*b-4*a*c

double rs = (-b - sqr(D)) / (2*a)
double xs = M + N * rs
double ys = P + Q * rs

result.center.x = xs
result.center.y = ys
result.radius   = rs

if ( s1 < 1 )
result.locator = @"Internal Tangent:"
else
result.locator = @"External Tangent:"
end if	
end fn = result

Circle c1, c2, c3, c

c1.center.x = 0.0 : c1.center.y = 0.0 : c1.radius = 1.0
c2.center.x = 4.0 : c2.center.y = 0.0 : c2.radius = 1.0
c3.center.x = 2.0 : c3.center.y = 4.0 : c3.radius = 2.0

// External tangent
c = fn SolveApollonius( c1, c2, c3,  1,  1,  1 )
NSLog( @"%@", fn CircleToString( c ) )

// Internal tangent
c = fn SolveApollonius( c1, c2, c3, -1, -1, -1 )
NSLog( @"%@", fn CircleToString( c ) )

HandleEvents