-- File: tkz_elements_circle.lua -- Copyright (c) 2026 Alain Matthes -- SPDX-License-Identifier: LPPL-1.3c -- Maintainer: Alain Matthes --------------------------------------------------------------------------- -- circles --------------------------------------------------------------------------- circle = {} circle.__index = circle function circle:new(c, t) -- c --> center t --> through local type = "circle" local ct = line:new(c, t) local opp = antipode_(c, t) local radius = point.abs(c - t) local south = c - point(0, radius) local east = c + point(radius, 0) local north = c + point(0, radius) local west = c - point(radius, 0) local perimeter = tkz.tau * radius local area = 0.5 * tkz.tau * radius * radius local cir = { center = c, through = t, ct = ct, opp = opp, radius = radius, south = south, east = east, north = north, west = west, type = type, perimeter = perimeter, area = area, } setmetatable(cir, self) return cir end setmetatable(circle, { __call = function(cls, ...) return cls:new(...) end, }) function circle:get(i) if i == nil then return self.center, self.through elseif i == 1 then return self.center elseif i == 2 then return self.through else return nil end end -- other definition function circle:radius(center, radius) return circle:new(center, center + point(radius, 0)) end function circle:diameter(za, zb) return circle:new(midpoint_(za, zb), zb) end ----------------------- -- Result -> boolean ----------------------- function circle:on_circle(p, EPS) return on_circle_(self.center, self.through, p, EPS) end function circle:in_disk(p, EPS) return in_disk_(self.center, self.through, p, EPS) end function circle:in_disk_strict(p, EPS) return in_disk_strict_(self.center, self.through, p, EPS) end function circle:out_disk_strict(p, EPS) return out_disk_strict_(self.center, self.through, p, EPS) end circle.in_out = circle.on_circle circle.in_out_disk = circle.in_disk circle.in_out_disk_strict = circle.in_disk_strict function circle:position(obj, EPS) EPS = EPS or tkz.epsilon local mt = getmetatable(obj) if mt == point then return circle_point_position_(self, obj, EPS) elseif mt == line then return line_circle_position_(obj, self, EPS) elseif mt == circle then return circles_position_(self.center, self.radius, obj.center, obj.radius, EPS) end error("tkz-elements: circle:position: a point, a line or a circle is expected", 2) end -- legacy alias (keep for backward compatibility) circle.point_position = circle.position -- Disk (closed): IN / OUT (boundary counted as IN) function circle:position_disk(p, EPS) EPS = EPS or tkz.epsilon local d = point.abs(p - self.center) local r = point.abs(self.through - self.center) if d <= r + EPS then return "IN" else return "OUT" end end -- Boolean shortcuts (same default tolerance) function circle:is_disjoint(L) return self:line_position(L) == "DISJOINT" end function circle:is_tangent(L, EPS) return self:line_position(L, EPS) == "TANGENT" end function circle:is_secant(L, EPS) return self:line_position(L, EPS) == "SECANT" end function circle:line_position(L, EPS) EPS = EPS or tkz.epsilon return line_position_(self.center, self.radius, L.pa, L.pb, EPS) end -- Relative position between two lines and a circle- function circle:lines_position(L1, L2, mode) local a, b = L1:get() local c, d = L2:get() local o = self.center local r = self.radius if mode == "parallel" or is_parallel_(a, b, c, d) then return lines_parallel_circle_(o, r, a, b, c, d) else return lines_secant_circle_(o, r, a, b, c, d) end end ------------------------ -- string -------------- ------------------------ function circle:circles_position(C, EPS) return circles_position_(self.center, self.radius, C.center, C.radius, EPS) end ----------------------- -- Result -> real ----------------------- function circle:power(pt) local d = point.abs(self.center - pt) return d * d - self.radius * self.radius end ----------------------- -- Result -> point ----------------------- function circle:antipode(pt) return 2 * self.center - pt end function circle:midarc(z1, z2) local phi = 0.5 * get_angle_normalize_(self.center, z1, z2) return rotation_(self.center, phi, z1) end function circle:point(t, P, EPS) local phi = 2 * t * math.pi local A = P or self.through EPS = EPS or tkz.epsilon -- the starting point must be on the circle if self:position(A, EPS) ~= "ON" then error("tkz-elements: circle:point: the point P is not on the circle", 2) end return rotation_(self.center, phi, A) end function circle:random(inside) inside = (inside == "inside") if inside then local x = self:point(math.random()) local L = line:new(self.center, x) return L:random() else return self:point(math.random()) end end function circle:internal_similitude(C) return internal_similitude_(self.center, self.radius, C.center, C.radius) end function circle:external_similitude(C) if math.abs(self.radius - C.radius) < tkz.epsilon then error("tkz-elements: circle:external_similitude: the radii are equal, " .. "the external center of similitude is at infinity", 2) end local S = external_similitude_(self.center, self.radius, C.center, C.radius) return S end function circle:similitude(C, mode) mode = mode or "internal" local s = (mode == "external") and -1 or 1 if s == -1 and math.abs(self.radius - C.radius) < tkz.epsilon then error("tkz-elements: circle:similitude: the radii are equal, " .. "the external center of similitude is at infinity", 2) end local c1, r1 = self.center, self.radius local c2, r2 = C.center, C.radius return barycenter_({ c2, r1 }, { c1, s * r2 }) end function circle:radical_center(C1, C2) if C2 == nil then if self.radius > C1.radius then return radical_center_(self.center, self.through, C1.center, C1.through) else return radical_center_(C1.center, C1.through, self.center, self.through) end else return radical_center3_(self, C1, C2) end end -- Pole of the line L with respect to the circle self -- Returns the point P (or nil, "INFINITE" if L passes through O) function circle:pole(L, EPS) EPS = EPS or tkz.epsilon local O = self.center local R = self.radius -- 1) H = projection of O on L local H = L:projection(O) if not H then return nil, "NO_PROJECTION" -- pathological case end -- 2) If L passes through O -> pole at infinity (direction orthogonal to L) if length_(O, H) < EPS then return nil, "INFINITE" -- the pole is at infinity end -- 3) P = inverse of H in the circle (center O, radius R) -- OP * OH = R^2, P on the ray OH local P = self:inversion(H) return P end ----------------------- -- Result -> line ----------------------- function circle:tangent_at(pt) return line:new(rotation_(pt, math.pi / 2, self.center), rotation_(pt, -math.pi / 2, self.center)) end function circle:tangent_from(pt) local t1, t2 = tangent_from_(self.center, self.through, pt) return line:new(pt, t1), line:new(pt, t2) end function circle:tangent_parallel(L) local LP = L:orthogonal_from(self.center) local X, Y = intersection(LP, self) return self:tangent_at(X), self:tangent_at(Y) end function circle:radical_axis(C) local t1, t2 if self.radius > C.radius then t1, t2 = radical_axis_(self.center, self.through, C.center, C.through) else t1, t2 = radical_axis_(C.center, C.through, self.center, self.through) end return line:new(t1, t2) end function circle:polar(p) local q = self:inversion(p) local qa = (p - q):orthogonal(1):at(q) local qb = (q - p):orthogonal(1):at(q) return line:new(qa, qb) end -- common_tangent(C, mode): common tangents of two circles, as lines. -- mode = "external" (default), "internal" or "both" (externals first). -- An explicit error is raised when no tangent of the requested family exists. function circle:common_tangent(C, mode, EPS) mode = mode or "external" if mode ~= "external" and mode ~= "internal" and mode ~= "both" and mode ~= "all" then error("tkz-elements: circle:common_tangent: mode must be \"external\", \"internal\" or \"both\", got " .. tostring(mode), 2) end local list, pos = common_tangents_(self, C, mode, EPS) if #list == 0 then local why = { IDENTICAL = "the two circles are identical", CONCENTRIC = "the two circles are concentric", DISJOINT_INT = "one circle is inside the other", SECANT = "the circles intersect", TANGENT_INT = "the circles are internally tangent", } error("tkz-elements: circle:common_tangent: no " .. (mode == "both" and "" or mode .. " ") .. "common tangent (" .. (why[pos] or pos) .. ")", 2) end local lines = {} for i, t in ipairs(list) do lines[i] = line:new(t[1], t[2]) end return table.unpack(lines) end -- Same as common_tangent, without error: nothing is returned when no -- tangent exists (used by the Apollonius solvers). function circle:common_tangent_or_nil_(C, mode, EPS) local list = common_tangents_(self, C, mode, EPS) local lines = {} for i, t in ipairs(list) do lines[i] = line:new(t[1], t[2]) end return table.unpack(lines) end ----------------------- -- Result -> circle ----------------------- --[[ Constructs a circle orthogonal to the current circle (`self`) and passing through a given external point `pt`. The resulting circle: - has its center at the given point `pt`, - is orthogonal to `self`, - and uses the distance from `pt` to a point of tangency as its radius. Geometric condition: - The given point `pt` must lie **outside** the current circle. - This is verified using `in_out_disk`. If the condition is not met, the function raises a TeX error. Returns: A new circle (center = `pt`, radius = distance to the tangent point). --]] function circle:orthogonal_from(pt) -- test if pt is outside the circle if self:in_out_disk(pt) then error("tkz-elements: circle:orthogonal_from: the point must lie outside the circle", 2) else -- Calculate tangents from point 'pt'. local t1, t2 = tangent_from_(self.center, self.through, pt) -- Return a circle with the center in 'pt' and one of the tangents as the radius return circle:new(pt, t1) end end -- orthogonal_through(pta, ptb [, X] [, EPS]) -- Circle orthogonal to self and passing through pta and ptb. -- If pta and ptb are inverse with respect to self, every circle through them -- is orthogonal to self: the circle through X is returned if a third point X -- is given, otherwise the circle with diameter [pta, ptb]. function circle:orthogonal_through(pta, ptb, X, EPS) if type(X) == "number" then -- old call: orthogonal_through(pta, ptb, EPS) EPS, X = X, nil end EPS = EPS or tkz.epsilon if X ~= nil and (type(X) ~= "table" or X.type ~= "point") then error("tkz-elements: circle:orthogonal_through: the third argument must be a point", 2) end if point.abs(pta - ptb) <= EPS then error("tkz-elements: circle:orthogonal_through: the two points must be distinct", 2) end -- Case where the three points are collinear if is_linear_(self.center, pta, ptb) then -- the two points are inverses of each other: a family of solutions if are_inverses_(self.center, self.through, pta, ptb, EPS) then if X then if is_linear_(pta, ptb, X) then error("tkz-elements: circle:orthogonal_through: the third point is on the line " .. "through the two points", 2) end return circle:new(circum_center_(pta, ptb, X), pta) end return circle:new(midpoint_(pta, ptb), pta) else error("tkz-elements: circle:orthogonal_through: the points are collinear with the " .. "center but not inverse: no orthogonal circle", 2) end else local o = orthogonal_through_(self.center, self.through, pta, ptb, EPS) if not o then error("tkz-elements: circle:orthogonal_through: no orthogonal circle through these points", 2) end return circle:new(o, pta) end end -- alias retained circle.orthogonal_circle_through = circle.orthogonal_through function circle:radical_circle(C1, C2) local rc = self:radical_center(C1, C2) if C2 == nil then return self:orthogonal_from(rc) else return C1:orthogonal_from(rc) end end function circle:midcircle(obj, EPS) if type(obj) == "table" and obj.type == "circle" then -- Returns the midcircle of 'self' and 'C' return midcircle_cc_(self.center, self.through, obj.center, obj.through, EPS) elseif type(obj) == "table" and obj.type == "line" then return midcircle_cl_(self.center, self.through, obj.pa, obj.pb, EPS) else error("tkz-elements: circle:midcircle: a circle or a line is expected", 2) end end ----------------------------------------------------------- -- Result -> drawing ----------------------------------------------------------- function circle:draw() -- Coordinates of the center and radius of the circle local x, y = self.center:get() local r = self.radius -- Format of the LaTeX command that draws the circle local frmt = "\\draw (%0.3f,%0.3f) circle [radius=%0.3f];" -- Print the LaTeX command with the values of x, y and r tex.sprint(string.format(frmt, x, y, r)) end function circle:inversion_L(L) -- Check whether the center of the circle is on the line L if L:in_out(self.center) then return L -- Returns the line L unchanged if the center is on it else -- Projection of the center on the line L local p = L:projection(self.center) -- Inverse of the projected point in the circle (self.center, self.through) local q = inversion_(self.center, self.through, p) -- Returns the circle with center the midpoint of 'self.center' and 'q', through 'q' return circle:new(midpoint_(self.center, q), q) end end function circle:inversion_C(C) local p, q, x, y, X, Y if C:in_out(self.center) then p = C:antipode(self.center) q = inversion_(self.center, self.through, p) x = ortho_from_(q, self.center, p) y = ortho_from_(q, p, self.center) return line:new(x, y) else x, y = intersection_lc_(self.center, C.center, C.center, C.through) X = inversion_(self.center, self.through, x) Y = inversion_(self.center, self.through, y) return circle:new(midpoint_(X, Y), X) end end -- Inversion in the circle self. Each argument may be a point, a line -- or a circle; a line or a circle may be sent to a circle or a line. -- z.P = C:inversion(z.M) ; z.P, C.new = C:inversion(z.M, L.AB) function circle:inversion(...) local n = select("#", ...) local t = {} for i = 1, n do local obj = select(i, ...) local ty = type(obj) == "table" and obj.type or type(obj) if ty == "point" then if point.abs(obj - self.center) <= tkz.epsilon then error("tkz-elements: circle:inversion: the inversion is undefined at the center", 2) end t[i] = inversion_(self.center, self.through, obj) elseif ty == "line" then t[i] = self:inversion_L(obj) elseif ty == "circle" then t[i] = self:inversion_C(obj) else error("tkz-elements: circle:inversion: expected a point, a line or a circle, got a " .. tostring(ty), 2) end end return table.unpack(t, 1, n) end -- Alias simple : inversion de rapport < 0 function circle:inversion_neg(obj) return self.center:symmetry(self:inversion(obj)) end function circle:path(za, zb, nb) local P = path() for i = 0, nb do local t = i / nb local pt = self:point_arc(za, zb, t) P:add_point(pt) end return P end function circle:point_arc(za, zb, t) local anglezazb = get_angle_normalize_(self.center, za, zb) local phi = t * anglezazb return rotation_(self.center, phi, za) end ---=== PB CONTACT APOLLONIUS ===--- -- ----------------------------------------------------------- -- CPP CCP CLP CLL -- CCC(?) CCL(?) -- ================== CPP ========================== function circle:CPP(a, b) local o = self.center local th = self.through local u, v = mediator_(a, b) local pa_center = path:new() local pa_through = path:new() local n = 0 local function push(o, t) if o and t then pa_center:add_point(o) pa_through:add_point(t) n = n + 1 end end local ain = in_disk_(o, th, a) local bin = in_disk_(o, th, b) ------------------------------------------------------------ -- 0) One point inside and the other outside -> no solution ------------------------------------------------------------ if (ain and not bin) or (bin and not ain) then return pa_center, pa_through, 0 ------------------------------------------------------------ -- 1) Line (AB) tangent to the circle (self) ------------------------------------------------------------ elseif is_tangent_(a, b, o, th) and (not is_equidistant_(a, b, o)) then local w = circum_center_(a, b, o) local c, d = intersection_cc_(o, th, w, a) local i = intersection_ll_(c, d, a, b) local t, tp = tangent_from_(o, th, i) if is_parallel_(o, t, u, v) then local w1 = intersection_ll_(o, tp, u, v) push(w1, tp) elseif is_parallel_(o, tp, u, v) then local w1 = intersection_ll_(o, t, u, v) push(w1, t) else local w1 = intersection_ll_(u, v, o, t) local w2 = intersection_ll_(u, v, o, tp) push(w1, t) push(w2, tp) end return pa_center, pa_through, n ------------------------------------------------------------ -- 2) Line (AB) parallel to a tangent to the circle ------------------------------------------------------------ elseif is_equidistant_(a, b, o) then local t1, t2 = intersection_lc_(u, v, o, th) -- if (AB) is tangent to the circle, one of the points is on (AB): -- this degenerate solution is skipped for _, t in ipairs({ t1, t2 }) do if not tkz.aligned(a, b, t) then push(circum_center_(a, b, t), t) end end return pa_center, pa_through, n ------------------------------------------------------------ -- 3) Two EXTERNAL points: classical (general) case ------------------------------------------------------------ elseif (not ain) and (not bin) then local w = circum_center_(a, b, o) local c, d = intersection_cc_(o, th, w, a) local i = intersection_ll_(a, b, c, d) local t, tp = tangent_from_(o, th, i) local w1 = intersection_ll_(u, v, o, t) local w2 = intersection_ll_(u, v, o, tp) push(w1, t) push(w2, tp) return pa_center, pa_through, n ------------------------------------------------------------ -- 4) Two INTERNAL points -> inversion ------------------------------------------------------------ elseif ain and bin then local Cinv = circle:new(a, b) local Csol = Cinv:inversion(self) if not Csol then return pa_center, pa_through, 0 end local LTA, LTB = Csol:tangent_from(b) if LTA then local C1 = Cinv:inversion(LTA) if C1 then push(C1.center, a) end end if LTB then local C2 = Cinv:inversion(LTB) if C2 then push(C2.center, a) end end return pa_center, pa_through, n end end circle.c_c_pp = circle.CPP --====================== START of CCP ========--- ------------------------------------- -------------------------------------------------------- local function sp_c_cc_p(p, o1, t1, o2, t2, r1, r2, ex) local n = 0 local pa_center = path:new() local pa_through = path:new() local function push(o, t) if o and t then pa_center:add_point(o) pa_through:add_point(t) n = n + 1 end end -- NB: on suppose que report_ accepte (A,B,scale[,origin]) comme ailleurs. local x = report_(o1, o2, r1) -- point in the direction o1->o2, distance ~r1 local y = report_(o2, o1, r2) -- point in the direction o2->o1, distance ~r2 local o = circum_center_(p, x, y) -- center of the circle (p,x,y) -- intersections of the circles (o,p) and (o2,t2) local u, v = intersection_cc_(o, p, o2, t2) if not u or not v then return pa_center, pa_through, 0 end -- choose u such that u != y (swap otherwise) if u:identity(y) then u, v = v, u end -- ligne (y,u) croise (ex,p) en i local i = intersection_ll_(y, u, ex, p) if not i then return pa_center, pa_through, 0 end -- tangents from i to the circle (o2,t2) local tau1, tau2 = tangent_from_(o2, t2, i) -- renamed so as not to hide the parameters if not tau1 or not tau2 then return pa_center, pa_through, 0 end -- centers = intersections of the lines (o2,tau_k) and (p,o) local w1 = intersection_ll_(o2, tau1, p, o) local w2 = intersection_ll_(o2, tau2, p, o) if not w1 or not w2 then -- on pousse ce qui existe if w1 then push(w1, tau1) end if w2 then push(w2, tau2) end return pa_center, pa_through, n end -- order: w1 is the nearest to p if length_(w1, p) > length_(w2, p) then w1, w2, tau1, tau2 = w2, w1, tau2, tau1 end push(w1, tau1) push(w2, tau2) return pa_center, pa_through, n end -- =================================================== -- solve_intersect(Cself, C, p, mode) -- -> Circles tangent to Cself and C through p -- -> simple version, without inversion or external dependencies -- =================================================== -- Circle tangent to self and C, through p (general case) local function solve_intersect(Cself, C, p, mode) if Cself:on_circle(p) then local tA, tB = tangent_at_(Cself.center, p) local T = line:new(tA, tB) return C:CLP(T, p) elseif C:on_circle(p) then local tA, tB = tangent_at_(C.center, p) local T = line:new(tA, tB) return Cself:CLP(T, p) else local pa_center = path:new() local pa_through = path:new() local n = 0 local function push(o, t) if o and t then pa_center:add_point(o) pa_through:add_point(t) n = n + 1 end end local I = midpoint_(Cself.center, C.center) local K = circle:new(p, I) local A = K:inversion(Cself) local B = K:inversion(C) local T1, T2 = A:common_tangent_or_nil_(B, "external") local C1 = K:inversion(T1) local C2 = K:inversion(T2) local O1, T1 = C1:get() local O2, T2 = C2:get() push(O1, T1) push(O2, T2) return pa_center, pa_through, n end end -------------------------------------------------------- -- CCP: two internally tangent circles, solution circle through p -- (Uniformisation push/paths, logique intacte) -------------------------------------------------------- local function solve_inside_tangent(Cself, C, p, mode) if Cself:on_circle(p) then local tA, tB = tangent_at_(Cself.center, p) local T = line:new(tA, tB) return C:CLP(T, p) elseif C:on_circle(p) then local tA, tB = tangent_at_(C.center, p) local T = line:new(tA, tB) return Cself:CLP(T, p) else local T = intersection(Cself, C) -- point de tangence commun if p:identity(T) then return path:new(), path:new(), 0 -- P = T -> no solution end -- Inversion en p local K = circle:new(p, Cself.through or (p + point(1,0))) local A = K:inversion(Cself) -- usually a line (p on Cself?) local B = K:inversion(C) local L1, L2 = A:common_tangent_or_nil_(B, "all") local function valid(L) return L and (not L:on_line(p)) end local pa_center = path:new() local pa_through = path:new() local n = 0 local function push(o, t) if o and t then pa_center:add_point(o) pa_through:add_point(t) n = n + 1 end end if valid(L1) then local S1 = K:inversion(L1) if S1 and (S1.radius or 0) > tkz.epsilon then push(S1.center, S1.through) end end if valid(L2) then local S2 = K:inversion(L2) if S2 and (S2.radius or 0) > tkz.epsilon then push(S2.center, S2.through) end end return pa_center, pa_through, n end end -------------------------------------------------------- -- CCP: case "inside" -------------------------------------------------------- local function solve_inside(Cself, C, p, mode) local a = Cself:point_position(p) local b = C:point_position(p) -- Special cases (no solution) if (a == "OUT" and b == "ON") or (a == "IN" and b == "IN") then return path:new(), path:new(), 0 else -- Circle of inversion K local K = circle:new(C.through, Cself.through) local A = K:inversion(Cself) local B = K:inversion(C) -- Inverse of the point p local ip = K:inversion(p) -- Output paths local pa_center = path:new() local pa_through = path:new() local n = 0 local function push(o, t) if o and t then pa_center:add_point(o) pa_through:add_point(t) n = n + 1 end end local pc_inv, pt_inv, m = A:CLP(B, ip) if pc_inv and pt_inv and m and m > 0 then for i = 1, m do local o_inv = pc_inv:get(i) local t_inv = pt_inv:get(i) if o_inv and t_inv then -- Circle in the inverted plane local Cinv = circle:new(o_inv, t_inv) -- Back to the original plane local S = K:inversion(Cinv) push(S.center, S.through) end end end return pa_center, pa_through, n end end -------------------------------------------------------- -- CCP: case "outside" -- (same structure, with safeguards) -------------------------------------------------------- -- ==================== SPECIAL CASE r1 == r2 ==================== local function solve_outside_equal_radii(Cself, C, p, mode) mode = mode or "external" local pa_center = path:new() local pa_through = path:new() local n = 0 local function push(o, t) if o and t then pa_center:add_point(o) pa_through:add_point(t) n = n + 1 end end local M if C:on_circle(Cself.through) then M = antipode_(Cself.center, Cself.through) else M = Cself.through end local K = circle:new(M, C.center) local A = K:inversion(Cself) local B = K:inversion(C) local ip = K:inversion(p) local pc_inv, pt_inv, m = B:CLP(A, ip) if pc_inv and pt_inv and m then for i = 1, m do local oi = pc_inv:get(i) local ti = pt_inv:get(i) if oi and ti then local S = K:inversion(circle:new(oi, ti)) if S and S.center and S.through then push(S.center, S.through) end end end end return pa_center, pa_through, n end -- ================== FIN CAS r1 == r2 ================== local function solve_outside(Cself, C, p, mode) local o1, t1, r1 = Cself.center, Cself.through, Cself.radius local o2, t2, r2 = C.center, C.through, C.radius local pos1 = Cself:point_position(p) -- "IN", "ON", "OUT" local pos2 = C:point_position(p) if pos1 == "IN" or pos2 == "IN" then return path:new(), path:new(), 0 elseif pos1 == "ON" then local tA, tB = tangent_at_(o1, p) -- two points on the tangent to self at p local T = line:new(tA, tB) return C:CLP(T, p) elseif pos2 == "ON" then local tA, tB = tangent_at_(o2, p) -- two points on the tangent to self at p local T = line:new(tA, tB) return Cself:CLP(T, p) else if mode == "all" then local PCe, PTe, ne = Cself:CCP(C, p, "external") local PCi, PTi, ni = Cself:CCP(C, p, "internal") return (PCe + PCi), (PTe + PTi), (ne + ni) end -- ======= DELEGATION IF r1 == r2 ======= if tkz.approx(r1, r2) then return solve_outside_equal_radii(Cself, C, p, mode) end if mode == "internal" then local i = Cself:internal_similitude(C) local x = report_(o1, o2, -r1, o1) local y = report_(o2, o1, r2, o2) local o = circum_center_(x, y, p) local r, s = intersection_lc_(i, p, o, p) -- r != p if not r or not s then return path:new(), path:new(), 0 end if r:identity(p) then r, s = s, r end return Cself:CPP(p, r) else -- external local Cmid = Cself:midcircle(C) local ex = Cmid.center local q = Cmid:inversion(p) if p:identity(q) then return sp_c_cc_p(p, o1, t1, o2, t2, r1, r2, ex) else return Cself:CPP(p, q) end end end end -------------------------------------------------------- -- CCP dispatcher -------------------------------------------------------- function circle:CCP(C, p, mode) mode = mode or "all" local pos = self:circles_position(C) if pos == "DISJOINT_EXT" then return solve_outside(self, C, p, mode) elseif pos == "SECANT" then return solve_intersect(self, C, p, mode) elseif pos == "TANGENT_EXT" then return solve_outside(self, C, p, mode) elseif pos == "TANGENT_INT" then return solve_inside_tangent(self, C, p, mode) elseif pos == "DISJOINT_INT" then return solve_inside(self, C, p, mode) elseif pos == "CONCENTRIC" or pos == "IDENTICAL" then return false -- or specific handling else error("Unknown circle position in CCP: "..tostring(pos)) end end ----==== END of CCP ====== --===================== -- ======= START of CLP --===================== -- Circle tangent to one circle, on line and passing through a point -- Returns: PA.center, PA.through, n (paths of centers / through-points, and count) --===================== -- CLP : circle tangent to a line and passing through a point -- Retour : pa_center, pa_through, n --===================== function circle:CLP(l, p, which) which = which or "all" local pos_line = self:line_position(l) local os, ts, rs = self.center, self.through, self.radius local ip = projection_(l.pa, l.pb, p) -- or l:projection(p) local i = projection_(l.pa, l.pb, os) -- or l:projection(os) local no = report_(os, i, -rs) local so = report_(os, i, rs) local pos = intersection_ll_(l.pa, l.pb, os, p) local sO, sP = l:side_line(os), l:side_line(p) -- chemins + compteur local pa_center, pa_through = path:new(), path:new() local n = 0 local function push(o, t) if o and t then pa_center:add_point(o) pa_through:add_point(t) n = n + 1 end end if pos_line == "DISJOINT" then -- (1) P = i if p:identity(i) then local o1 = midpoint_(p, no) local o2 = midpoint_(p, so) push(o1, p) push(o2, p) return pa_center, pa_through, n -- (2) P on l (P != i) elseif l:on_line(p) then if not l.length or l.length == 0 then return pa_center, pa_through, n end local k = rs / l.length local u = orthogonal_at_(l.pa, l.pb, p, k) local v = orthogonal_at_(l.pa, l.pb, p, -k) local ux, uy = mediator_(os, u) local vx, vy = mediator_(os, v) local o1 = intersection_ll_(u, v, ux, uy) local o2 = intersection_ll_(u, v, vx, vy) push(o1, p) push(o2, p) return pa_center, pa_through, n -- (3) P on the wrong side elseif sO ~= sP then return pa_center, pa_through, n -- (4) P strictly inside the disk elseif self:in_disk_strict(p) then return pa_center, pa_through, n -- (5) P on C elseif self:on_circle(p) then if p:identity(no) then push(midpoint_(no, i), p) return pa_center, pa_through, n elseif p:identity(so) then push(midpoint_(so, i), p) return pa_center, pa_through, n else local i1 = intersection_ll_(l.pa, l.pb, no, p) local i2 = intersection_ll_(l.pa, l.pb, so, p) if i1 then local j1 = orthogonal_at_(l.pa, l.pb, i1) local o1 = intersection_ll_(i1, j1, os, p) push(o1, p) end if i2 then local j2 = orthogonal_at_(l.pa, l.pb, i2) local o2 = intersection_ll_(i2, j2, os, p) push(o2, p) end return pa_center, pa_through, n end -- (6) General case: inversion at P else local Cinv = circle:new(p, ip) local C1 = Cinv:inversion(self) local C2 = Cinv:inversion(l) if which == "external" or which == "all" then local L1, L2 = C1:common_tangent_or_nil_(C2, "external") if L1 then local S = Cinv:inversion(L1) if S and S.center and S.through then push(S.center, S.through) end end if L2 then local S = Cinv:inversion(L2) if S and S.center and S.through then push(S.center, S.through) end end end if which == "internal" or which == "all" then local L1, L2 = C1:common_tangent_or_nil_(C2, "internal") if L1 then local S = Cinv:inversion(L1) if S and S.center and S.through then push(S.center, S.through) end end if L2 then local S = Cinv:inversion(L2) if S and S.center and S.through then push(S.center, S.through) end end end return pa_center, pa_through, n end elseif pos_line == "TANGENT" then -- (1) P = i if p:identity(i) then push(os, i) return pa_center, pa_through, n -- (2) P on l (tangent configuration) elseif l:on_line(p) then local Cinv = circle:new(i, p) local Linv = Cinv:inversion(self) -- image of C: a line local pac, pat, m = l:LLP(Linv, p) for k = 1, (m or 0) do local o = pac:get(k) local t = pat:get(k) if o and t then local S1 = Cinv:inversion(circle:new(o, t)) if S1 then local oc, tc = S1:get() push(oc, tc) end end end return pa_center, pa_through, n -- (3) P on the wrong side elseif sO ~= sP then return pa_center, pa_through, n -- (4) P strictly inside the disk elseif self:in_disk_strict(p) then return pa_center, pa_through, n -- (5) P on C, at the singular point (no/so) elseif self:on_circle(p) and (p:identity(no) or p:identity(so)) then push(os, ts) return pa_center, pa_through, n -- (6) P on C, general case elseif self:on_circle(p) then local Cinv = circle:new(i, p) local Lself = Cinv:inversion(self) local mpt = midpoint_(p, ip) local S1 = Cinv:inversion(circle:new(mpt, p)) if S1 then local oc, tc = S1:get() push(oc, tc) end return pa_center, pa_through, n -- (7) General tangent case else local Cinv = circle:new(i, p) local Lself = Cinv:inversion(self) local pac, pat, m = l:LLP(Lself, p) for k = 1, (m or 0) do local o = pac:get(k) local t = pat:get(k) if o and t then local S1 = Cinv:inversion(circle:new(o, t)) if S1 then local oc, tc = S1:get() push(oc, tc) end end end return pa_center, pa_through, n end else -- "secant" if l:on_line(p) and self:on_circle(p) then return pa_center, pa_through, 0 elseif self:on_circle(p) then local Cinv = circle:new(p, ip) local Lself = Cinv:inversion(self) -- line local C2 = Cinv:inversion(l) -- circle if not (Lself and C2 and C2.center and C2.radius) then return pa_center, pa_through, n end local d = (Lself.pb - Lself.pa) local dn = point.abs(d) if dn <= tkz.epsilon then return pa_center, pa_through, n end -- normale unitaire : n = perp(d)/|d| local nux = (d.im / dn) local nuy = -(d.re / dn) local nvec = point(nux, nuy) -- Two points at distance r2 from the center of C2, along the normal to Lself local O2 = C2.center local r2 = C2.radius local Q1 = O2 + r2 * nvec local Q2 = O2 - r2 * nvec -- Two lines parallel to Lself ("tangent" to Lself in the required sense) and tangent to C2 local L1 = line:new(Q1, Q1 + d) local L2 = line:new(Q2, Q2 + d) -- Inverse again -> two solution circles (through P, tangent to C and l) if L1 then local S = Cinv:inversion(L1); if S and S.center and S.through then push(S.center, S.through) end end if L2 then local S = Cinv:inversion(L2); if S and S.center and S.through then push(S.center, S.through) end end return pa_center, pa_through, n elseif l:on_line(p) then local _, tp = through(p, 1) local Cinv = circle:new(p, tp) local C1 = Cinv:inversion(self) if not C1 or not C1.center or not C1.radius then return pa_center, pa_through, n end local o1, r1 = C1.center, C1.radius if not l.length or l.length == 0 then return pa_center, pa_through, n end local k = r1 / l.length local q1 = orthogonal_at_(l.pa, l.pb, o1, k) local q2 = orthogonal_at_(l.pa, l.pb, o1, -k) local dir = (l.pb - l.pa) local L1 = line:new(q1, q1 + dir) local L2 = line:new(q2, q2 + dir) if L1 then local S = Cinv:inversion(L1) if S and S.center and S.through then push(S.center, S.through) end end if L2 then local S = Cinv:inversion(L2) if S and S.center and S.through then push(S.center, S.through) end end return pa_center, pa_through, n else local Cinv = circle:new(p, ip) local C1 = Cinv:inversion(self) -- circle local C2 = Cinv:inversion(l) -- circle (p is not on l) local Le1, Le2 = C1:common_tangent_or_nil_(C2, "external") if Le1 then local S = Cinv:inversion(Le1) if S and S.center and S.through then push(S.center, S.through) end end if Le2 then local S = Cinv:inversion(Le2) if S and S.center and S.through then push(S.center, S.through) end end return pa_center, pa_through, n end end end -- Aliases circle.c_lc_p = circle.CLP circle.c_cl_p = circle.CLP --===================== -- ======= END of CLP --===================== --======================================== -- CLL: circle tangent to the two lines L1, L2 --======================================== --============================== -- Helper: degenerate circle? --============================== local function is_degenerate_circle(C, I) if not C then return true end if math.abs(C.radius or 0) < tkz.epsilon then return true end if C.through and C.through:identity(I) then return true end return false end local function circle_outside_strip_(o, r, a, b, c, d) local d1 = distance_(a, b, o) -- distance to L1 local d2 = distance_(c, d, o) -- distance to L2 local D = distance_(a, b, c) -- distance between L1 and L2 local EPS = 0.00001 if math.abs((d1 + d2) - D) <= EPS then return false end if d1 > r + EPS and d2 > r + EPS then return true else return false end end function circle:CLL_parallel(L1, L2) local pa_center = path:new() local pa_through = path:new() local n = 0 local function push(o, t) if o and t then pa_center:add_point(o) pa_through:add_point(t) n = n + 1 end end local o = self.center local t = self.through local r = self.radius local a, b = L1:get() local c, d = L2:get() -- missing: circle externally tangent outside the strip if circle_outside_strip_(o, r, a, b, c, d) then return path:new(), path:new(), 0 else local u = projection_(c, d, a) local v = projection_(c, d, b) local x = midpoint_(a, u) local y = midpoint_(b, v) local d = length_(a, u) / 2 local ta = report_(o, t, self.radius + d, o) local w1, w2 = intersection_lc_(x, y, o, ta) local t1 = projection_(a, b, w1) local t2 = projection_(a, b, w2) push(w1, t1) push(w2, t2) return pa_center, pa_through, n end end -- first sector u,i,v ; then v,i,up ; up,i,vp ; vp,i,u function circle:CLL(L1, L2, choice, inside) choice = choice or "all" local a, b = L1:get() local c, d = L2:get() if is_parallel_(a, b, c, d) then return self:CLL_parallel(L1, L2) else if choice == "all" then return self:CLL_all(L1, L2) else local a, b = L1:get() local c, d = L2:get() local o = self.center local r = self.radius local i = intersection_ll_(a, b, c, d) -- Relative orientation of the lines local s = (((b - a) ^ (d - c)) >= 0) and 1 or -1 local u = report_(a, b, 1, i) local up = report_(a, b, -1, i) local v = report_(c, d, s, i) local vp = report_(c, d, -s, i) -- Choice of the sector (1..4) local x, y if choice == 1 then x, y = u, v elseif choice == 2 then x, y = v, up elseif choice == 3 then x, y = up, vp elseif choice == 4 then x, y = vp, u else error("tkz-elements: circle:CLL: the choice must be 1, 2, 3, 4 or \"all\", got " .. tostring(choice), 2) end -- inside: the sign of the radius is changed inside = (inside == "inside") if inside then r = -r end -- Two lines at distance +-r from I towards x/y local LL1 = line:new(collinear_at_distance_(i, x, -r)) local LL2 = line:new(collinear_at_distance_(i, y, r)) -- LLP (paths version): centers on LL2, tangent to LL1, through o local pc, pt, m = LL2:c_ll_p(LL1, o) -- Solution circles (projection on the line (i,x)) local S1, S2 = nil, nil if (m or 0) >= 1 then local w1 = pc:get(1) local t1 = line(i, x):projection(w1) S1 = circle:new(w1, t1) end if (m or 0) >= 2 then local w2 = pc:get(2) local t2 = line(i, x):projection(w2) S2 = circle:new(w2, t2) end -- Remove degenerate circles local d1 = S1 and is_degenerate_circle(S1, i) local d2 = S2 and is_degenerate_circle(S2, i) if d1 and (not d2) and S2 then S1, S2 = S2, S1 d1, d2 = d2, d1 elseif d1 and d2 then S1, S2 = nil, nil end -- Store in paths with push (as in CLP) local pa_center = path:new() local pa_through = path:new() local n = 0 local function push(o, t) if o and t then pa_center:add_point(o) pa_through:add_point(t) n = n + 1 end end if S1 then push(S1.center, S1.through) end if S2 then push(S2.center, S2.through) end return pa_center, pa_through, n end end end circle.c_cll = circle.CLL --===================================================== -- CLL_all: all the solutions, sector by sector --===================================================== function circle:CLL_all(L1, L2) local sectors = self:lines_position(L1, L2) local pa_center = path:new() local pa_through = path:new() local n = 0 -- push (as in CLP) local function push(o_, t_) if o_ and t_ then pa_center:add_point(o_) pa_through:add_point(t_) n = n + 1 end end -- (pc, pt, m) -> local paths via push local function push_from_paths(pc, pt, m) if not pc or not pt or not m then return end for i = 1, m do push(pc:get(i), pt:get(i)) end end if not sectors or #sectors == 0 then return pa_center, pa_through, 0 end local function add_choice(choice, inside) local pc, pt, m = self:CLL(L1, L2, choice, inside) push_from_paths(pc, pt, m) end if #sectors == 1 then local m = sectors[1] add_choice(m, nil) -- 2 add_choice(m, "inside") -- 2 -> total 4 elseif #sectors == 2 then local m, n = sectors[1], sectors[2] -- adjacents add_choice(m, nil) -- 2 add_choice(n, nil) -- 2 -> total 4 elseif #sectors == 3 then local l, m, n = sectors[1], sectors[2], sectors[3] add_choice(l, nil) -- 2 add_choice(n, nil) -- 2 add_choice(m, nil) -- 2 add_choice(m, "inside") -- 2 -> total 8 else -- {1,2,3,4} -> 8 solutions (2 per sector) add_choice(1, nil) add_choice(2, nil) add_choice(3, nil) add_choice(4, nil) end return pa_center, pa_through, n end -- ========================================================= -- ========================================================= function circle:CCL_DDD(C2, D) local pa_center, pa_through = path:new(), path:new() local n = 0 local function push(o, t) if o and t then pa_center:add_point(o) pa_through:add_point(t) n = n + 1 end end local function circle_with_new_radius_keep_center(C, newR) local Tp = report_(C.center, C.through, newR) return circle:new(C.center, Tp) end local C1, O1, R1 = self, self.center, self.radius local O2, R2 = C2.center, C2.radius local EPS = 0.0001 local function signed_side_(L, P) local a, b = L:get() local z = (b - a) ^ (P - a) if math.abs(z) <= EPS then return 0 end return (z > 0) and 1 or -1 end local function line_separates_disks_(L, O1, R1, O2, R2) local s1 = signed_side_(L, O1) local s2 = signed_side_(L, O2) if s1 == 0 or s2 == 0 or s1 * s2 >= 0 then return false -- centers on the same side or on the line -> not "between" end local d1 = D:distance(O1) local d2 = D:distance(O2) return (d1 > R1 + EPS) and (d2 > R2 + EPS) -- strictly outside the two disks end if line_separates_disks_(D, O1, R1, O2, R2) then return pa_center, pa_through, 0 end ---------------------------------------------------------- -- ====== FAMILY R1 - R2: line translated by -R2 ====== ---------------------------------------------------------- do local newR = R1 - R2 local pc, pt, m local Lplus = D:collinear_at_distance(-R2) if math.abs(newR) < EPS then pc, pt, m = Lplus:LPP(O1, O2) else local Cminus = circle_with_new_radius_keep_center(C1, newR) pc, pt, m = Cminus:CLP(Lplus, O2) end for i = 1, m do local wi = pc:get(i) local ti = pt:get(i) local r0 = length_(wi, ti) local r = r0 - R2 local Ti_new = report_(wi, ti, r) local di = D:distance(wi) local li = length_(wi, Ti_new) if act_(O1, C1.through, wi, Ti_new) and act_(O2, C2.through, wi, Ti_new) and (math.abs(li - di) < EPS) then push(wi, Ti_new) end end end ---------------------------------------------------------- -- ====== FAMILY R1 - R2: line translated by +R2 ====== ---------------------------------------------------------- do local newR = R1 - R2 local pc, pt, m local Lplus = D:collinear_at_distance(R2) if math.abs(newR) < EPS then pc, pt, m = Lplus:LPP(O1, O2) else local Cminus = circle_with_new_radius_keep_center(C1, newR) pc, pt, m = Cminus:CLP(Lplus, O2) end for i = 1, m do local wi = pc:get(i) local ti = pt:get(i) local r0 = length_(wi, ti) local r = r0 + R2 local Ti_new = report_(wi, ti, r) local di = D:distance(wi) local li = length_(wi, Ti_new) if act_(O1, C1.through, wi, Ti_new) and act_(O2, C2.through, wi, Ti_new) and (math.abs(li - di) < EPS) then push(wi, Ti_new) end end end ---------------------------------------------------------- -- ====== FAMILY R1 + R2: line translated by -R2 ====== ---------------------------------------------------------- do local newR = R1 + R2 local Cplus = circle_with_new_radius_keep_center(C1, newR) local Lplus = D:collinear_at_distance(-R2) local pc, pt, m = Cplus:CLP(Lplus, O2) for i = 1, m do local wi = pc:get(i) local ti = pt:get(i) local r0 = length_(wi, ti) local r = r0 - R2 -- same principle as R1-R2 / -R2 local Ti_new = report_(wi, ti, r) local di = D:distance(wi) local li = length_(wi, Ti_new) if act_(O1, C1.through, wi, Ti_new) and act_(O2, C2.through, wi, Ti_new) and (math.abs(li - di) < EPS) then push(wi, Ti_new) end end end ---------------------------------------------------------- -- ====== FAMILY R1 + R2: line translated by +R2 ====== ---------------------------------------------------------- do local newR = R1 + R2 local Cplus = circle_with_new_radius_keep_center(C1, newR) local Lminus = D:collinear_at_distance( R2) local pc, pt, m = Cplus:CLP(Lminus, O2) for i = 1, m do local wi = pc:get(i) local ti = pt:get(i) local r0 = length_(wi, ti) local r = r0 + R2 -- same principle as R1-R2 / +R2 local Ti_new = report_(wi, ti, r) local di = D:distance(wi) local li = length_(wi, Ti_new) if act_(O1, C1.through, wi, Ti_new) and act_(O2, C2.through, wi, Ti_new) and (math.abs(li - di) < EPS) then push(wi, Ti_new) end end end return pa_center, pa_through, n end --============================================ ---======== CCL circles tangent inside --============================================ function circle:CCL_ITT(C2, D) local pa_center = path:new() local pa_through = path:new() local count = 0 local function push(o, t) if o and t then pa_center:add_point(o) pa_through:add_point(t) count = count + 1 end end local C1 = self local c1, t1 = C1:get() local c2, t2 = C2:get() local r1 = point.abs(c1 - t1) local r2 = point.abs(c2 - t2) -- circle of inversion: center O (intersection of the two circles), radius OC2 local O = tangent_point_two_circles_(c1, r1, c2, r2) local Cinv = circle:new(O, c2) -- inversion of the two circles and of the line D local L1 = Cinv:inversion(C1) -- ligne image de C1 local L2 = Cinv:inversion(C2) -- ligne image de C2 local CD = Cinv:inversion(D) -- circle, image of the line D -- circles tangent to CD and to the two lines L1, L2 local pc, pt, nsol = CD:CLL(L1, L2) -- each solution of CLL is inverted back for i = 1, nsol do local W = pc:get(i) local T = pt:get(i) if W and T then local CW = circle:new(W, T) -- solution circle in the inverted plane local w, t = Cinv:inversion(CW):get() -- back to the original plane push(w, t) end end return pa_center, pa_through, count end function circle:CCL_ISD(C2, D) local C1 = self local pa_center, pa_through = path:new(), path:new() local n = 0 local function push(o, t) if o and t then pa_center:add_point(o) pa_through:add_point(t) n = n + 1 end end local c1, t1 = C1:get() local a, b = D:get() local I1, I2 = intersection_lc_(a, b, c1, t1) local I = I1 or I2 local Cinv = circle:new(I, c1) local L1 = Cinv:inversion(C1) local C2i = Cinv:inversion(C2) -- circle C2' local pc, pt, m = C2i:CLL(L1, D) for i = 1, m do local W = pc:get(i) local T = pt:get(i) local Ctmp = circle:new(W, T) local S = Cinv:inversion(Ctmp) local O, Q = S:get() push(O, Q) end return pa_center, pa_through, n end function circle:CCL(C2, D, EPS) EPS = EPS or tkz.epsilon local pa_center, pa_through = path:new(), path:new() local n = 0 local function push(o, t) if o and t then pa_center:add_point(o) pa_through:add_point(t) n = n + 1 end end local C1 = self local c1 = C1.center local r1 = C1.radius local c2 = C2.center local r2 = C2.radius local posCC = circles_position_(c1, r1, c2, r2, EPS) local posC1L = line_position_(c1, r1, D.pa, D.pb, EPS) local posC2L = line_position_(c2, r2, D.pa, D.pb, EPS) -------------------------------------------------- if (posCC == "DISJOINT_INT" or posCC == "TANGENT_INT") and posC1L == "DISJOINT" then return path:new(), path:new(), 0 elseif (posCC == "DISJOINT_INT") and (posC1L == "SECANT") and ((posC2L == "DISJOINT") or (posC2L == "SECANT")) then return self:CCL_ISD(C2, D) elseif posCC == "TANGENT_INT" and posC2L == "TANGENT" then return self:CCL_ITT(C2, D) else return self:CCL_DDD(C2, D) end end ---===== START CCC by Viete =====------- -- === Test de tangence robuste (absolu ET relatif) ============================ -- === CCC: all the solutions (0..8) ====================================== function circle:CCC_viete_core(C2, C3, opts) opts = opts or {} local C1, R3, W3 = self, C3.radius, C3.center local pc, pt, n = path:new(), path:new(), 0 local function push(o,t) pc:add_point(o); pt:add_point(t); n = n + 1 end -- remove duplicates: compare center and radius with a tolerance local function exists_(o,t) local r = length_(o,t) for i=1,n do local oi, ti = pc:get(i), pt:get(i) if oi and ti then local same_center = length_(o, oi) <= (opts.abs or 1e-4) local same_radius = math.abs(r - length_(oi,ti)) <= (opts.abs or 1e-4) if same_center and same_radius then return true end end end return false end local function push_unique(o,t) if o and t and not exists_(o,t) then push(o,t) end end -- variant: keep the center, add dR to the radius local function grow_keep_center(C, dR) return circle:new(C.center, report_(C.center, C.through, dR, C.through)) end -- robust test "radius ~ 0" (absolute OR relative) local function radius_is_zero_(r, ref) local eps_abs = (opts.abs ~= nil) and opts.abs or 1e-4 local eps_rel = (opts.rel ~= nil) and opts.rel or 1e-6 return (math.abs(r) <= eps_abs) or (math.abs(r)/math.max(ref,1.0) <= eps_rel) end -- helper: from a list of centers (path pcv, nv), -- build the two lifts +/- R3 local function lift_from_centers_(pcv, nv) for i=1,nv do local O = pcv:get(i) if O then local Tplus = report_(O, W3, R3, W3) -- r = |OW3| + R3 (internally tangent to C3) local Tminus = report_(O, W3, -R3, W3) -- r = |OW3| - R3 (externally tangent to C3) if Tplus then local S = circle:new(O, Tplus) if act_(S.center, S.through, C1.center, C1.through, opts) and act_(S.center, S.through, C2.center, C2.through, opts) and act_(S.center, S.through, C3.center, C3.through, opts) then push_unique(O, Tplus) end end if Tminus then local S = circle:new(O, Tminus) if act_(S.center, S.through, C1.center, C1.through, opts) and act_(S.center, S.through, C2.center, C2.through, opts) and act_(S.center, S.through, C3.center, C3.through, opts) then push_unique(O, Tminus) end end end end end -- The 4 branches of Viete (+-R3 on C1 and C2) local mods = { {-R3,-R3}, {-R3, R3}, { R3,-R3}, { R3, R3} } for _, m in ipairs(mods) do local C1T = grow_keep_center(C1, m[1]) -- (O1, R1 + s1*R3) local C2T = grow_keep_center(C2, m[2]) -- (O2, R2 + s2*R3) local r1t = C1T.radius local r2t = C2T.radius local ref = math.max(C1.radius, C2.radius, R3, 1.0) local z1 = radius_is_zero_(r1t, ref) local z2 = radius_is_zero_(r2t, ref) if z1 and z2 then -- === PPP: both images are points => circle through O1, O2, O3 local pcv, ptv, nv = C1.center:PPP(C2.center, W3) lift_from_centers_(pcv, nv) elseif z1 then -- === CPP: only one image is a point => tangent to C2T and through O1 and O3 local pcv, ptv, nv = C2T:CPP(C1.center, W3) -- only the centers are used, to keep the same lifting logic lift_from_centers_(pcv, nv) elseif z2 then -- === CPP: only one image is a point => tangent to C1T and through O2 and O3 local pcv, ptv, nv = C1T:CPP(C2.center, W3) lift_from_centers_(pcv, nv) else -- === classical CCP: circle (preimage) through W3, tangent to C1T and C2T local pcv, _, nv = C1T:CCP(C2T, W3, "all") lift_from_centers_(pcv, nv) end end return pc, pt, n end ------------------------------------------------------------ -- CCC dispatcher: classification of the 3 circles -> resolution -- - "non-secant" cases -> CCC_viete_ok -- - "secant" cases -> dedicated treatment (CCC_secant) ------------------------------------------------------------ -- === Tolerances (abs + rel) ==================================== local function approx_(a, b, opts) opts = opts or {} local eps_abs = (opts.abs ~= nil) and opts.abs or 1e-4 local eps_rel = (opts.rel ~= nil) and opts.rel or 1e-6 local d = math.abs(a - b) return (d <= eps_abs) or (d <= eps_rel * math.max(math.abs(a), math.abs(b), 1.0)) end -- === Signature of the three pairs ================================= -- Relation between two circles Ci, Cj from circles_position_ local function circle_relation_from_positions(Ci, Cj, EPS) EPS = EPS or tkz.epsilon local c1, r1 = Ci.center, Ci.radius local c2, r2 = Cj.center, Cj.radius -- first: "same center" test, to separate concentric / equal circles local d = point.mod(c1 - c2) local same_center = tkz.approx(d, 0, EPS) if same_center then if tkz.approx(r1, r2, EPS) then return "equal" else return "concentric" end end -- otherwise, circles_position_ is used return circles_position_(c1, r1, c2, r2, EPS) end -- === Signature of the three pairs ================================= local function three_circles_signature(C1, C2, C3, opts) opts = opts or {} local EPS = opts.EPS or opts.eps or tkz.epsilon -- relation for each pair local r12 = circle_relation_from_positions(C1, C2, EPS) local r23 = circle_relation_from_positions(C2, C3, EPS) local r31 = circle_relation_from_positions(C3, C1, EPS) -- comptage global comme avant local counts = { secant = 0, tangent = 0, disjoint = 0, inside = 0, concentric = 0, equal = 0, } local function acc(tag) if tag == "SECANT" then counts.secant = counts.secant + 1 elseif tag == "TANGENT_EXT" or tag == "TANGENT_INT" then counts.tangent = counts.tangent + 1 elseif tag == "DISJOINT_EXT" then counts.disjoint = counts.disjoint + 1 elseif tag == "DISJOINT_INT" then counts.inside = counts.inside + 1 elseif tag == "CONCENTRIC" then counts.concentric = counts.concentric + 1 elseif tag == "IDENTICAL" then counts.equal = counts.equal + 1 else -- optional, but useful to avoid silent failures -- error("three_circles_signature: unexpected tag "..tostring(tag)) end end acc(r12) acc(r23) acc(r31) return { r12 = r12, r23 = r23, r31 = r31, counts = counts } end local function CCC_viete_ok(C1, C2, C3, opts) return C1:CCC_viete_core(C2, C3, opts) end -- === 2) Secant cases ========= -- CCC case: the three circles are pairwise secant -- Strategy: inversion at X = C1 n C2 -> C1, C2 become lines -- solve (S3:CLL(L1,L2)) or (L1:LLL(L2,L3)) -> inverse back local function CCC_secant(C1, C2, C3, opts) opts = opts or {} -- --- outputs (always returned) local pa_center = path:new() local pa_through = path:new() local n = 0 -- --- single push local function push(o, t) if o and t then pa_center:add_point(o) pa_through:add_point(t) n = n + 1 end end -- --- 1) intersections of C1 and C2 (they exist in this case) local X, Y = intersection(C1, C2) if (not X) or (not Y) then -- safeguard: not secant -> not handled here return pa_center, pa_through, 0 end -- --- 2) circle of inversion: center X, arbitrary radius (XY) local Cinv = circle:new(X, Y) -- --- 3) inverted objects local L1 = Cinv:inversion(C1) -- attendu : ligne local L2 = Cinv:inversion(C2) -- attendu : ligne local S3 = Cinv:inversion(C3) -- circle OR line (if C3 passes through X) -- --- helper: inverse back a solution given by (o_i, t_i) in the inverted plane local function push_back_from_inverted(o_i, t_i) if (not o_i) or (not t_i) then return end local S_i = circle:new(o_i, t_i) -- solution circle in the inverted plane local S_o = Cinv:inversion(S_i) -- inverse back -> solution circle local t_out = Cinv:inversion(t_i) -- inverse back the tangency (through) point if S_o and S_o.center and t_out then push(S_o.center, t_out) end end -- --- 4) Resolution in the inverted plane, then inverse back if S3 and S3.center then -- 4.a) S3 is a circle: CLL problem (circle tangent to the lines L1, L2 and to S3) -- All the solutions (sectors); see opts.which / opts.inside. local pc, pt, m = S3:CLL(L1, L2, "all") if m and m > 0 then for i = 1, m do push_back_from_inverted(pc:get(i), pt:get(i)) -- pt[i] is the point of tangency on S3 end end elseif S3 then -- 4.b) S3 is a line (C3 passes through X): LLL problem (three lines) local pc, pt, m = L1:LLL(L2, S3, "all") if m and m > 0 then for i = 1, m do push_back_from_inverted(pc:get(i), pt:get(i)) -- pt[i]: foot on L1 end end else -- inversion of C3 failed (very unlikely) -> no solution return pa_center, pa_through, 0 end return pa_center, pa_through, n end -- === 3) Special case: "two tangent + one disjoint" ============== local function CCC_two_tangent_one_disjoint(C1, C2, C3, sig, opts) -- Ex: r12 tangent, r23 tangent, r31 disjoint (or any permutation) return CCC_viete_ok(C1, C2, C3, opts) end local function three_distinct_radii(C1, C2, C3, opts) opts = opts or {} local eps = opts.EPS or opts.eps or tkz.epsilon local R1, R2, R3 = C1.radius, C2.radius, C3.radius return (not tkz.approx(R1, R2, eps)) and (not tkz.approx(R2, R3, eps)) and (not tkz.approx(R3, R1, eps)) end -- === 4) Dispatcher principal ===================================== function circle:CCC(C2, C3, opts) local C1 = self local sig = three_circles_signature(C1, C2, C3, opts) local k = sig.counts -- (A) If AT LEAST one pair is secant -> separate treatment if k.secant > 0 then return CCC_secant(C1, C2, C3, opts) -- (B) Special case "two tangent + one disjoint" elseif (k.tangent == 2) and (k.disjoint == 1) and (k.inside == 0) and (k.concentric == 0) then return CCC_two_tangent_one_disjoint(C1, C2, C3, sig, opts) elseif (k.tangent == 3) and ( (sig.r12 == "inside tangent" and sig.r23 == "outside tangent" and sig.r31 == "outside tangent") or (sig.r23 == "inside tangent" and sig.r31 == "outside tangent" and sig.r12 == "outside tangent") or (sig.r31 == "inside tangent" and sig.r12 == "outside tangent" and sig.r23 == "outside tangent") ) then return C1:CCC_gergonne(C2, C3, opts) else return CCC_viete_ok(C1, C2, C3, opts) end end ----================================= ----===========