-- File: tkz_elements_triangle.lua
-- Copyright (c) 2026 Alain Matthes
-- SPDX-License-Identifier: LPPL-1.3c
-- Maintainer: Alain Matthes

triangle = {}
triangle.__index = triangle


function triangle:new(za, zb, zc)
	local eps = tkz.epsilon
	-- ===== Guardrails against degeneration =====
	-- sall vertices combined
	if point.abs(za - zb) <= eps or point.abs(zb - zc) <= eps or point.abs(zc - za) <= eps then
		tex.error("Degenerate triangle: two vertices coincide.")
		return
	end
	-- aligned vertices
	if tkz.is_linear(za, zb, zc) then
		tex.error("Degenerate triangle: vertices aligned.")
		return
	end

	-- ===== Side lengths =====
	local c = point.abs(zb - za)
	local a = point.abs(zc - zb)
	local b = point.abs(za - zc)

	-- ===== Line objects =====
	local ab = line:new(za, zb)
	local bc = line:new(zb, zc)
	local ca = line:new(zc, za)

	-- ===== Angles =====
	local alpha = math.abs(get_angle_(za, zb, zc))
	local beta = math.abs(get_angle_(zb, zc, za))
	local gamma = math.abs(get_angle_(zc, za, zb))

	local alpha_ = angle:new(za, zb, zc)
	local beta_ = angle:new(zb, zc, za)
	local gamma_ = angle:new(zc, za, zb)

	local type = "triangle"
	local circumcenter = circum_center_(za, zb, zc)
	local centroid = barycenter_({ za, 1 }, { zb, 1 }, { zc, 1 })
	local incenter = in_center_(za, zb, zc)
	local orthocenter = ortho_center_(za, zb, zc)
	local eulercenter = euler_center_(za, zb, zc)
	local spiekercenter = spieker_center_(za, zb, zc)
	local orientation = find_orientation_(za, zb, zc)

	local semiperimeter = (a + b + c) / 2
	local cross = (zb - za) ^ (zc - za)
	local area = 0.5 * math.abs((zb - za) ^ (zc - za))
	local inradius = area / semiperimeter
	local circumradius = (a * b * c) / (4 * area)
  local ma = midpoint_(zb, zc)
	local mb = midpoint_(zc, za)
	local mc = midpoint_(za, zb)
  local ha = projection_(zb, zc, za)
	local hb = projection_(zc, za, zb)
	local hc = projection_(za, zb, zc)
	local ia = find_foot_bisector_(za, zb, zc)
	local ib = find_foot_bisector_(zb, zc, za)
	local ic = find_foot_bisector_(zc, za, zb)
	local tr = {
		pa = za,
		pb = zb,
		pc = zc,
		type = type,
		circumcenter = circumcenter,
		centroid = centroid,
		incenter = incenter,
		eulercenter = eulercenter,
		orthocenter = orthocenter,
		spiekercenter = spiekercenter,
		a = a,
		b = b,
		c = c,
		ab = ab,
		ca = ca,
		bc = bc,
		ma = ma,
		mb = mb,
		mc = mc,
		ha = ha,
		hb = hb,
		hc = hc,
		ia = ia,
		ib = ib,
		ic = ic,
		alpha = alpha,
		beta = beta,
		gamma = gamma,
		alpha_ = alpha_,
		beta_ = beta_,
		gamma_ = gamma_,
		cross = cross,
		semiperimeter = semiperimeter,
		area = area,
		inradius = inradius,
		circumradius = circumradius,
		orientation = orientation,
	}
	setmetatable(tr, self)
	return tr
end

setmetatable(triangle, {
	__call = function(cls, ...)
		return cls:new(...)
	end,
})

function triangle:get(i)
	if i == nil then
		return self.pa, self.pb, self.pc
	elseif i == 1 then
		return self.pa
	elseif i == 2 then
		return self.pb
	elseif i == 3 then
		return self.pc
	else
		return nil
	end
end

-----------------------
-- Result -> boolean
-----------------------
-- added strict
function triangle:position(p, EPS)
	EPS = EPS or tkz.epsilon
	local x, y, z = barycentric_coordinates_(self.pa, self.pb, self.pc, p)

	local in_  = (x >  EPS and y >  EPS and z >  EPS)
	local out_ = (x < -EPS or  y < -EPS or  z < -EPS)

	if in_ then
		return "IN"
	elseif out_ then
		return "OUT"
	else
		return "ON" -- zone de bord (tolérance)
	end
end


function triangle:in_out(pt, strict)
	-- ne change pas → tous les anciens documents restent valides
	return in_out_(self.pa, self.pb, self.pc, pt, strict)
end

function triangle:on_triangle(pt, EPS)
	EPS = EPS or tkz.epsilon
	local a, b, c = self.pa, self.pb, self.pc
	local x, y, z = barycentric_coordinates_(a, b, c, pt)
	-- boundary cases
	local onx = math.abs(x) <= EPS
	local ony = math.abs(y) <= EPS
	local onz = math.abs(z) <= EPS

	if onx or ony or onz then
		return true
	end

	return false
end

function triangle:check_equilateral(EPS)
	return check_equilateral_(self.pa, self.pb, self.pc, EPS)
end
triangle.is_equilateral = triangle.check_equilateral

function triangle:check_acutangle(EPS)
return check_acutangle_(self.pa, self.pb, self.pc, EPS)
end
triangle.is_acute = triangle.check_acutangle

function triangle:check_isosceles(EPS)
return check_isosceles_(self.pa, self.pb, self.pc, EPS)
end

triangle.is_isosceles = triangle.check_isosceles

-----------------------
-- Result -> real
-----------------------

function triangle:area() -- obsolete
	return area_(self.pa, self.pb, self.pc)
end


function triangle:barycentric_coordinates(pt)
	local xA, yA = self.pa:get()
	local xB, yB = self.pb:get()
	local xC, yC = self.pc:get()
	local xP, yP = pt:get()
	local orient2d_ABC = oriented_area2_(xA, yA, xB, yB, xC, yC)
	local l_A = oriented_area2_(xP, yP, xB, yB, xC, yC) / orient2d_ABC
	local l_B = oriented_area2_(xA, yA, xP, yP, xC, yC) / orient2d_ABC
	local l_C = oriented_area2_(xA, yA, xB, yB, xP, yP) / orient2d_ABC
	return l_A, l_B, l_C
end

function triangle:trilinear_coordinates(pt)
	local la, lb, lc = self:barycentric_coordinates(pt)
	local ta, tb, tc = la / self.a, lb / self.b, lc / self.c
	-- normalize
	local N = ta + tb + tc
	return ta / N, tb / N, tc / N
end

function triangle:trilinear_to_d(x, y, z)
	return trilinear_to_d_(x, y, z, self.a, self.b, self.c)
end

function triangle:get_angle(arg)
	local i = resolve_triangle_index(self, arg)  -- 0,1,2
	if i < 0 or i > 2 then
		tex.error("triangle:get_angle -> invalid argument (use 0/1/2, 'a'/'b'/'c', or a vertex point)")
		return nil
	end

	if i == 0 then return self.alpha end  -- angle en B dans ton alpha = angle(za,zb,zc)
	if i == 1 then return self.beta  end
	return self.gamma
end


function triangle:parameter(p)
	local a, b, c = self.a, self.b, self.c -- a=BC, b=CA, c=AB
	local A, B, C = self.pa, self.pb, self.pc
	local AB, BC, CA = self.ab, self.bc, self.ca
	local ptotal = a + b + c

	if AB:on_segment(p) then
		return length_(A, p) / ptotal
	elseif BC:on_segment(p) then
		return (c + length_(B, p)) / ptotal
	elseif CA:on_segment(p) then
		return (c + a + length_(C, p)) / ptotal
	else
		tex.error("The point is not on the boundary of the triangle.")
	end
end

-----------------------
-- Result -> point
-----------------------
function triangle:barycenter(ka, kb, kc)
	return barycenter_({ self.pa, ka }, { self.pb, kb }, { self.pc, kc })
end
triangle.barycentric = triangle.barycenter

function triangle:base(u, v) -- (ab,ac) base coord u,v
	return barycenter_({ self.pa, (1 - u - v) }, { self.pb, u }, { self.pc, v })
end

function triangle:point(t)
	local p = self.a + self.b + self.c -- périmètre
	local t1 = self.c / p -- fin de AB
	local t2 = (self.c + self.a) / p -- fin de AB + BC

	if t <= t1 then
		return self.ab:point(t / t1)
	elseif t <= t2 then
		return self.bc:point((t - t1) / (t2 - t1))
	else
		return self.ca:point((t - t2) / (1 - t2))
	end
end

function triangle:random(inside)
	inside = (inside == "inside")
	if inside then
		local x = self.bc:random()
		local L = line:new(self.pa, x)
		return L:random()
	else
		return self:point(math.random())
	end
end

function triangle:trilinear(a, b, c)
	return barycenter_({ self.pa, a * self.a }, { self.pb, b * self.b }, { self.pc, c * self.c })
end

function triangle:kimberling(n)
	local cos, sin, tan, pi = math.cos, math.sin, math.tan, math.pi
	local A, B, C = self.alpha, self.beta, self.gamma
	local a, b, c = self.a, self.b, self.c
	local pi = math.pi
	if n == 1 then
		return self:trilinear(1, 1, 1) -- incenter
	elseif n == 2 then
		return self:barycentric(1, 1, 1) -- centroid
	elseif n == 3 then
		return self:trilinear(cos(A), cos(B), cos(C)) --circumcenter
	elseif n == 4 then
		return self:barycentric(tan(A), tan(B), tan(C)) -- orthocenter
	elseif n == 5 then
		return self:trilinear(cos(B - C), cos(C - A), cos(A - B)) --nine
	elseif n == 6 then
		return self:trilinear(a, b, c) -- lemoine
	elseif n == 7 then
		return self:barycentric(1 / (b + c - a), 1 / (a + c - b), 1 / (a + b - c)) -- gergonne
	elseif n == 8 then
		return self:barycentric(b + c - a, a + c - b, a + b - c) -- nagel
	elseif n == 9 then
		return self:barycentric(a * (b + c - a), b * (a + c - b), c * (a + b - c)) -- mittenpunkt
	elseif n == 10 then
		return self:barycentric(b + c, a + c, a + b) -- spieker
	elseif n == 11 then
		return self:trilinear(1 - cos(B - C), 1 - cos(C - A), 1 - cos(A - B)) -- feuerbach
	elseif n == 13 then
		return self:trilinear(1 / sin(A + pi / 3), 1 / sin(B + pi / 3), 1 / sin(C + pi / 3)) -- first fermat
	elseif n == 14 then
		return self:trilinear(1 / sin(A - pi / 3), 1 / sin(B - pi / 3), 1 / sin(C - pi / 3)) -- second fermat
	elseif n == 19 then
		return self:trilinear(tan(A), tan(B), tan(C)) -- clawson
	elseif n == 20 then
		return self:barycentric(tan(B) + tan(C) - tan(A), tan(A) + tan(C) - tan(B), tan(A) + tan(B) - tan(C)) --de Longchamps
	elseif n == 39 then
		return self:trilinear(a * (b ^ 2 + c ^ 2), b * (c ^ 2 + a ^ 2), c * (a ^ 2 + b ^ 2)) --brocard midpoint
	elseif n == 55 then
		return self:trilinear(a * (b + c - a), b * (c + a - b), c * (a + b - c))
	elseif n == 56 then
		return self:trilinear(a / (b + c - a), b / (c + a - b), c / (a + b - c))
	elseif n == 110 then
		return self:trilinear(a / (b ^ 2 - c ^ 2), b / (c ^ 2 - a ^ 2), c / (a ^ 2 - b ^ 2)) -- kiepert parabola
	elseif n == 115 then
		return self:barycentric((b ^ 2 - c ^ 2) ^ 2, (c ^ 2 - a ^ 2) ^ 2, (a ^ 2 - b ^ 2) ^ 2) -- kiepert hyperbola
	elseif n == 175 then
		local ra, rb, rc = radius_excircle_(self.pa, self.pb, self.pc)
		return self:barycentric(a - ra, b - rb, c - rc)
	elseif n == 176 then
		ra, rb, rc = radius_excircle_(self.pa, self.pb, self.pc)
		return self:barycentric(a + ra, b + rb, c + rc)
	elseif n == 181 then
	  return self:trilinear((a * (b + c) ^ 2) / (b + c - a), (b * (c + a) ^ 2) / (c + a - b), (c * (a + b) ^ 2) / (a + b - c))--apollonius point
  elseif n == 213 then
		return self:trilinear((b + c) * a ^ 2, (c + a) * b ^ 2, (a + b) * c ^ 2)
	elseif n == 264 then
		return self:barycentric(1 / sin(2 * A), 1 / sin(2 * B), 1 / sin(2 * C))
	elseif n == 371 then
		return self:trilinear(cos(A - pi / 4), cos(B - pi / 4), cos(C - pi / 4))
	elseif n == 598 then
		return self:trilinear(
			b * c / (a ^ 2 - 2 * b ^ 2 - 2 * c ^ 2),
			c * a / (b ^ 2 - 2 * c ^ 2 - 2 * a ^ 2),
			a * b / (c ^ 2 - 2 * a ^ 2 - 2 * b ^ 2)
		) -- to get lemoine triangle
	else
		tex.error({ " Your number is not yet in the list" })
	end
end

function triangle:isogonal(pt)
	local pa = self.bc:reflection(pt)
	local pb = self.ca:reflection(pt)
	local pc = self.ab:reflection(pt)
	return circum_center_(pa, pb, pc)
end

function triangle:excenter(pt)
	if pt == self.pa then
		return self:trilinear(-1, 1, 1)
	elseif pt == self.pb then
		return self:trilinear(1, -1, 1)
	elseif pt == self.pc then
		return self:trilinear(1, 1, -1)
	else
		tex.error("The argument must be a vertex of the triangle")
	end
end

function triangle:projections(p)
	local p1 = projection_(self.pb, self.pc, p)
	local p2 = projection_(self.pa, self.pc, p)
	local p3 = projection_(self.pa, self.pb, p)
	return p1, p2, p3
end


function triangle:reflections(p)
	local p1 = symmetry_axial_(self.pa, self.pb, p)
	local p2 = symmetry_axial_(self.pb, self.pc, p)
	local p3 = symmetry_axial_(self.pc, self.pa, p)
	return p2, p3, p1
end

function triangle:parallelogram()
	local x = self.pc + self.pa - self.pb
	return x
end

function triangle:bevan_point()
	return circum_center_(excentral_tr_(self.pa, self.pb, self.pc))
end

function triangle:mittenpunkt_point()
	return lemoine_point_(excentral_tr_(self.pa, self.pb, self.pc))
end

function triangle:gergonne_point()
	return gergonne_point_(self.pa, self.pb, self.pc)
end

function triangle:nagel_point()
	return nagel_point_(self.pa, self.pb, self.pc)
end

function triangle:feuerbach_point()
	return feuerbach_point_(self.pa, self.pb, self.pc)
end

function triangle:symmedian_point()
	return lemoine_point_(self.pa, self.pb, self.pc)
end

triangle.lemoine_point = triangle.symmedian_point
triangle.grebe_point = triangle.symmedian_point

function triangle:spieker_center()
	return spieker_center_(self.pa, self.pb, self.pc)
end

function triangle:euler_points()
	local H = self.orthocenter
	return midpoint_(H, self.pa), midpoint_(H, self.pb), midpoint_(H, self.pc)
end

function triangle:orthic_points()
	return  orthic_tr_(self.pa, self.pb, self.pc)
end

function triangle:medial_points()
	return medial_tr_(self.pa, self.pb, self.pc)
end

function triangle:tangential_points()
	return tangential_tr_(self.pa, self.pb, self.pc)
end

function triangle:nine_points()
	local ma, mb, mc, ha, hb, hc
	-- Calculate the medial triangle
	ma, mb, mc = medial_tr_(self.pa, self.pb, self.pc)
	-- Calculate the orthic triangle
	ha, hb, hc = orthic_tr_(self.pa, self.pb, self.pc)
	-- Calculate the orthocenter
	local H = self.orthocenter
	-- Return the points of the nine-point circle
	return ma, mb, mc, ha, hb, hc, midpoint_(H, self.pa), midpoint_(H, self.pb), midpoint_(H, self.pc)
end

function triangle:first_lemoine_points()
	local lc = self:lemoine_point()
	local p1 = intersection_ll_(lc, ll_from_(lc, self.pa, self.pb), self.pa, self.pc)
	local p2 = intersection_ll_(lc, ll_from_(lc, self.pb, self.pc), self.pb, self.pa)
	local p3 = intersection_ll_(lc, ll_from_(lc, self.pc, self.pa), self.pc, self.pb)
	local p5 = intersection_ll_(lc, ll_from_(lc, self.pa, self.pb), self.pb, self.pc)
	local p6 = intersection_ll_(lc, ll_from_(lc, self.pb, self.pc), self.pa, self.pc)
	local p4 = intersection_ll_(lc, ll_from_(lc, self.pc, self.pa), self.pa, self.pb)
	return p1, p2, p3, p4, p5, p6
end

function triangle:second_lemoine_points()
	local lc = self:lemoine_point()
	local ha, hb, hc = orthic_tr_(self.pa, self.pb, self.pc)
	local p1 = intersection_ll_(lc, ll_from_(lc, ha, hb), self.pa, self.pc)
	local p2 = intersection_ll_(lc, ll_from_(lc, hb, hc), self.pb, self.pa)
	local p3 = intersection_ll_(lc, ll_from_(lc, hc, ha), self.pc, self.pb)
	local p5 = intersection_ll_(lc, ll_from_(lc, ha, hb), self.pb, self.pc)
	local p6 = intersection_ll_(lc, ll_from_(lc, hb, hc), self.pa, self.pc)
	local p4 = intersection_ll_(lc, ll_from_(lc, hc, ha), self.pa, self.pb)
	return p1, p2, p3, p4, p5, p6
end

function triangle:soddy_points()
	return soddy_center_(self.pa, self.pb, self.pc)
end

-- soddy centers
function triangle:soddy_center(outer)
	outer = (outer == "outer")
	local index = outer and 175 or 176
	return self:kimberling(index)
end

function triangle:conway_points()
	local a1 = report_(self.pb, self.pa, length_(self.pb, self.pc), self.pa)
	local a2 = report_(self.pc, self.pa, length_(self.pb, self.pc), self.pa)
	local b1 = report_(self.pa, self.pb, length_(self.pa, self.pc), self.pb)
	local b2 = report_(self.pc, self.pb, length_(self.pa, self.pc), self.pb)
	local c1 = report_(self.pb, self.pc, length_(self.pb, self.pa), self.pc)
	local c2 = report_(self.pa, self.pc, length_(self.pb, self.pa), self.pc)
	return a1, a2, b1, b2, c1, c2
end

function triangle:conway_center()
	return self.incenter
end

function triangle:first_fermat_point()
	local x = equilateral_tr_(self.pb, self.pa)
	local y = equilateral_tr_(self.pa, self.pc)
	return intersection_ll_(x, self.pc, y, self.pb)
end

function triangle:second_fermat_point()
	local x = equilateral_tr_(self.pa, self.pb)
	local y = equilateral_tr_(self.pc, self.pa)
	return intersection_ll_(x, self.pc, y, self.pb)
end

function triangle:orthic_axis_points()
	return orthic_axis_(self.pa, self.pb, self.pc)
end

function triangle:kenmotu_point()
	return self:kimberling(371)
end

function triangle:kenmotu_points()
	local ck = self:kenmotu_circle()
	local p1, p2 = intersection_lc_(self.pa, self.pb, ck.center, ck.through)
	local p3, p4 = intersection_lc_(self.pb, self.pc, ck.center, ck.through)
	local p5, p6 = intersection_lc_(self.pc, self.pa, ck.center, ck.through)
	return p1, p2, p3, p4, p5, p6
end

function triangle:taylor_points()
	local a, b, c = orthic_tr_(self.pa, self.pb, self.pc)
	return projection_(self.pa, self.pc, a),
		projection_(self.pa, self.pb, a),
		projection_(self.pb, self.pa, b),
		projection_(self.pb, self.pc, b),
		projection_(self.pc, self.pb, c),
		projection_(self.pc, self.pa, c)
end

function triangle:adams_points()
	local i = self:kimberling(1)
	local g = self:kimberling(7)
	local l = projection_(self.pa, self.pc, i)
	local n = projection_(self.pa, self.pb, i)
	local m = projection_(self.pb, self.pc, i)
	local x = ll_from_(g, m, n)
	local y = ll_from_(g, n, l)
	local z = ll_from_(g, m, l)
	local t1 = intersection_ll_(x, g, self.pa, self.pb)
	local t4 = intersection_ll_(x, g, self.pb, self.pc)
	local t2 = intersection_ll_(y, g, self.pa, self.pb)
	local t5 = intersection_ll_(y, g, self.pa, self.pc)
	local t3 = intersection_ll_(z, g, self.pb, self.pc)
	local t6 = intersection_ll_(z, g, self.pa, self.pc)
	return t1, t2, t3, t4, t5, t6
end

function triangle:lamoen_points()
	local ma, mb, mc = medial_tr_(self.pa, self.pb, self.pc)
	local g = centroid_(self.pa, self.pb, self.pc)
	local p1 = circum_center_(self.pa, mc, g)
	local p2 = circum_center_(self.pb, ma, g)
	local p3 = circum_center_(self.pb, mc, g)
	local p4 = circum_center_(self.pc, mb, g)
	local p5 = circum_center_(self.pc, ma, g)
	local p6 = circum_center_(self.pa, mb, g)
	return p1, p2, p3, p4, p5, p6
end

function triangle:yiu_centers()
	local ap = symmetry_axial_(self.pb, self.pc, self.pa)
	local bp = symmetry_axial_(self.pa, self.pc, self.pb)
	local cp = symmetry_axial_(self.pa, self.pb, self.pc)
	local oa = circum_center_(self.pa, bp, cp)
	local ob = circum_center_(self.pb, ap, cp)
	local oc = circum_center_(self.pc, bp, ap)
	return oa, ob, oc
end

function triangle:first_brocard_point()
	return self:trilinear(self.c / self.b, self.a / self.c, self.b / self.a)
end

function triangle:second_brocard_point()
	return self:trilinear(self.b / self.c, self.c / self.a, self.a / self.b)
end

function triangle:brocard_midpoint()
	return self:kimberling(39)
end

function triangle:macbeath_point()
	return self:kimberling(264)
end

function triangle:poncelet_point(pt)
	local A, B, C, D = self.pa, self.pb, self.pc, pt
	local e1 = euler_center_(A, B, C)
	local e2 = euler_center_(A, C, D)
	local m = midpoint_(A, C)
	local x, y = intersection_cc_(e1, m, e2, m)

	if not x or not y then
		tex.error("Poncelet construction failed: no intersection.")
	end

	if x == m then
		return y
	else
		return x
	end
end

function triangle:orthopole(l, EPS)
	EPS = EPS or tkz.epsilon

	-- garde-fous
	if tkz.is_linear(self.pa, self.pb, self.pc) then
		tex.error("orthopole: degenerate triangle")
		return
	end

	local p, q = l:get()
	if point.abs(p - q) <= EPS then
		tex.error("orthopole: degenerate line")
		return
	end

	-- projections de A,B sur l (on n'a pas besoin de C ici)
	local ap = projection_(p, q, self.pa)
	local bp = projection_(p, q, self.pb)

	-- pieds sur les côtés CA et BC
	local a1, a2 = self.bc:get()   -- droite (BC)
	local c1, c2 = self.ca:get()   -- droite (CA)

	if point.abs(a2 - a1) <= EPS or point.abs(c2 - c1) <= EPS then
		tex.error("orthopole: degenerate triangle side")
		return
	end

	local app = projection_(a1, a2, ap)  -- pied de ap sur (BC)
	local bpp = projection_(c1, c2, bp)  -- pied de bp sur (CA)

	-- orthopole = intersection des droites (ap, app) et (bp, bpp)
	return intersection(line:new(ap, app), line:new(bp, bpp))
end



function triangle:isodynamic_points()
	-- Cas particulier : triangle équilatéral
	if self:is_equilateral() then
		-- On prend le centre de gravité G = (A + B + C)/3
		local A, B, C = self.pa, self.pb, self.pc
		local G = (A + B + C) / 3
		-- On renvoie deux fois le même point pour garder la même signature
		return G, G
	end

	-- Cas général
	local a, b, c = self.a, self.b, self.c
	local ab, bc = self.ab, self.bc, self.ca  -- self.ca est inutile ici mais cohérent

	-- Cercles d’Apollonius associés aux rapports b/a sur AB et c/b sur BC
	local Cab = ab:apollonius(b / a)
	local Cbc = bc:apollonius(c / b)

	-- L’intersection donne les deux points isodynamiques
	return intersection(Cab, Cbc)
end

triangle.neuberg_points = triangle.isodynamic_points

-- Points d’intersection du cercle d’Apollonius (b/a, côté AB)
-- avec la droite (AB), en mettant en premier le point sur le segment [AB].

function triangle:apollonius_points(side, EPS)
	EPS = EPS or tkz.epsilon

	local a, b, c = self.a, self.b, self.c  -- a = BC, b = CA
	local A, B, C = self.pa, self.pb, self.pc
  local ab, bc, ca = self.ab, self.bc, self.ca
   local x, y, X, Y, L, k
		if side == "ab"  then
			 L = self.ab
			 k = b / a
			x, y = a, b
			X, Y = A, B     -- AC / BC
	elseif side == "bc"  then
			L = self.bc
			k = c / b
			x, y = b, c
			X, Y = B, C         -- AB / CA
	elseif side == "ca"  then
			L = self.ca
			k = a / c
			x, y = c, a
			X, Y = C, A           -- BC / AB
	else
			tex.error("triangle:apollonius_circle -> invalid side: use 'ab', 'bc', or 'ca'")
			return
	end

		if math.abs(x - y) <= EPS then
		local M = (X + Y) / 2
		return M, M
	end

     return L:harmonic(k)
end

-- Ancien nom si tu veux conserver la compatibilité :
-- triangle.bisector_points = triangle.apollonius_points

function triangle:the_apollonius_center()
	local lbrocard = self:brocard_axis()
	local n = self.eulercenter
	local s = self.spiekercenter
	local lns = line(n,s)
	return intersection(lbrocard, lns)
end

function triangle:apollonius_point()
	return kimberling(181)
end

function triangle:isotomic(p)
	local pta, ptb, ptc = cevian_(self.pa, self.pb, self.pc, p)

	local pta_s = 2*self.ma - pta
	local ptb_s = 2*self.mb - ptb

	return intersection_ll_(self.pa,pta_s,self.pb,ptb_s)
end

function triangle:anticomplement(p)
return self.centroid:homothety(-2,p)
end

function triangle:complement(p)
return self.centroid:homothety(-1/2,p)
end
-------------------
-- Result -> line
-------------------

function triangle:symmedian_line(arg)
	local pts = { self.pa, self.pb, self.pc }
	local i = resolve_triangle_index(self, arg)
	if i < 0 or i > 2 then
		return nil
	end

	local A = pts[i + 1]
	local B = pts[(i + 1) % 3 + 1]
	local C = pts[(i + 2) % 3 + 1]
	local L = self:lemoine_point()

	local P = intersection_ll_(A, L, B, C)
	return line:new(A, P)
end


function triangle:altitude(arg)
	local pts = { self.pa, self.pb, self.pc }
	local i = resolve_triangle_index(self, arg)
	if i < 0 or i > 2 then
		return nil
	end

	local A = pts[i + 1]
	local B = pts[(i + 1) % 3 + 1]
	local C = pts[(i + 2) % 3 + 1]

	return line:new(A, projection_(B, C, A))
end


function triangle:bisector(arg)
	local i = resolve_triangle_index(self, arg)
	if i < 0 or i > 2 then
		return nil
	end

	local pts = { self.pa, self.pb, self.pc }
	local I = self.incenter

	local A = pts[i + 1]
	local B = pts[(i + 1) % 3 + 1]
	local C = pts[(i + 2) % 3 + 1]

	local P = intersection_ll_(A, I, B, C)
	return line:new(A, P)
end


function triangle:bisector_ext(arg)
	local i = resolve_triangle_index(self, arg)
	if i < 0 or i > 2 then
		return nil
	end

	local pts = { self.pa, self.pb, self.pc }
	local A = pts[i + 1]
	local B = pts[(i + 1) % 3 + 1]
	local C = pts[(i + 2) % 3 + 1]

	return line:new(A, bisector_ext_(A, B, C))
end


function triangle:mediator(arg)
	local i = resolve_triangle_index(self, arg)
	if i < 0 or i > 2 then
		return nil
	end

	local pts = { self.pa, self.pb, self.pc }
	local B = pts[(i + 1) % 3 + 1]
	local C = pts[(i + 2) % 3 + 1]

	return line:new(B, C):mediator()
end


function triangle:trisector(arg)
	local pts = { self.pa, self.pb, self.pc }
	local i = resolve_triangle_index(self, arg)
	if i < 0 or i > 2 then
		return nil
	end

	local A = pts[i + 1]
	local B = pts[(i + 1) % 3 + 1]
	local C = pts[(i + 2) % 3 + 1]

	local an = get_angle_(A, B, C)
	local pti = report_(A, B, 1)
	local ptj = rotation_(A, an / 3, pti)
	local ptk = rotation_(A, an / 3, ptj)

	return line(A, ptj), line(A, ptk)
end

function triangle:antiparallel(pt, arg)
	local a, b, c = self.pa, self.pb, self.pc
	local I = self.incenter
	local i = resolve_triangle_index(self, arg)

	local S = { a, b, c }
	local A = S[i + 1] -- sommet de départ
	local B = S[(i + 1) % 3 + 1] -- premier sommet adjacent
	local C = S[(i + 2) % 3 + 1] -- second sommet adjacent

	-- u, v sont les symétriques de B et C par rapport à la bissectrice en I
	local u = symmetry_axial_(A, I, B)
	local v = symmetry_axial_(A, I, C)
	local w = ll_from_(pt, u, v)

	return line:new(intersection_ll_(pt, w, B, A), intersection_ll_(pt, w, C, A))
end

function triangle:orthic_axis()
	local x, y, z = orthic_axis_(self.pa, self.pb, self.pc)
	return line:new(x, z)
end

-- N,H,G,O -- check_equilateral_ (a,b,c)
function triangle:euler_line()
	if check_equilateral_(self.pa, self.pb, self.pc) then
		tex.error("An error has occurred", { "No euler line with equilateral triangle" })
	else
		local a, b, c = orthic_axis_(self.pa, self.pb, self.pc)
		local x = projection_(a, c, self.eulercenter)
		return line:new(self.eulercenter, x)
	end
end

-- with pt on the circumcircle
function triangle:steiner_line(pt)
	local u = symmetry_axial_(self.pa, self.pb, pt)
	local v = symmetry_axial_(self.pa, self.pc, pt)
	return line:new(u, v)
end

function triangle:lemoine_axis() -- revoir car problème
	local C = self:circum_circle()
	local pl = self:lemoine_point()
	return C:polar(pl)
end

function triangle:fermat_axis()
	local x = self:first_fermat_point()
	local y = self:second_fermat_point()
	return line:new(x, y)
end

function triangle:brocard_axis()
	return line:new(self.circumcenter, self:lemoine_point())
end

function triangle:simson_line(pt) -- pt on circumcircle
	local x = self.ab:projection(pt)
	local y = self.bc:projection(pt)
	return line:new(x, y)
end

function triangle:soddy_line()
	return line:new(self:soddy_center(), self:soddy_center("outer"))
end

-----------------------
--- Result -> circles --
-----------------------


function triangle:apollonius_circle(side, EPS)
		 EPS = EPS or tkz.epsilon
		-- lengths
		local a, b, c = self.a, self.b, self.c
		local L, k

		if side == "ab"  then
				L = self.ab
				k = b / a        -- AC / BC
		elseif side == "bc"  then
				L = self.bc
				k = c / b        -- AB / CA
		elseif side == "ca"  then
				L = self.ca
				k = a / c        -- BC / AB
		else
				tex.error("triangle:apollonius_circle -> invalid side: use 'ab', 'bc', or 'ca'")
				return
		end

		-- cas k ≈ 1 : volontairement traité
		if math.abs(k - 1) <= EPS then
    	return L:bisector()
		else
		  return L:apollonius(k)
		end
end



function triangle:three_apollonius_circles()
	local a, b, c = self.a, self.b, self.c
	local ab, bc, ca = self.ab, self.bc, self.ca
	local Cab = ab:apollonius(b / a)
	local Cbc = bc:apollonius(c / b)
	local Cca = ca:apollonius(a / c)
	return Cab, Cbc, Cca
end


function triangle:three_tangent_circles()
	local i = self.incenter
	local a, b, c = self:projection(i)
	return circle:new(self.pa, c), circle:new(self.pb, a), circle:new(self.pc, b)
end


function triangle:euler_circle()
	return circle:new(euler_center_(self.pa, self.pb, self.pc), midpoint_(self.pb, self.pc))
end

--aliase
triangle.feuerbach = triangle.euler_circle

function triangle:circum_circle()
	return circle:new(circum_center_(self.pa, self.pb, self.pc),self.pa)
end

function triangle:in_circle()
	return circle:new(self.incenter, projection_(self.pb, self.pc, self.incenter))
end
-- see excenter ??

function triangle:ex_circle(arg)
	local pts = { self.pa, self.pb, self.pc }
	local i = resolve_triangle_index(self, arg)
	if i < 0 or i > 2 then return nil end

	local A = pts[i + 1]
	local B = pts[(i + 1) % 3 + 1]
	local C = pts[(i + 2) % 3 + 1]

	local O = ex_center_(A, B, C)          -- excentre opposé à A
	local T = projection_(B, C, O)         -- point de tangence sur (BC)

	return circle:new(O, T)
end

function triangle:radical_circle()
	local Ca, Cb, Cc
	Ca = self:ex_circle(self.pa)
	Cb = self:ex_circle(self.pb)
	Cc = self:ex_circle(self.pc)
	return Ca:radical_circle(Cb, Cc)
end

function triangle:spieker_circle()
	local ma, mb, mc, sp
	ma, mb, mc = medial_tr_(self.pa, self.pb, self.pc)
	sp = in_center_(ma, mb, mc)
	return circle:new(sp, projection_(ma, mb, sp))
end

function triangle:cevian_circle(p)
	local pta, ptb, ptc = cevian_(self.pa, self.pb, self.pc, p)
	return circle:new(circum_center_(pta, ptb, ptc), pta)
end

function triangle:symmedial_circle()
	local pta, ptb, ptc, p
	p = lemoine_point_(self.pa, self.pb, self.pc)
	pta, ptb, ptc = cevian_(self.pa, self.pb, self.pc, p)
	return circle:new(circum_center_(pta, ptb, ptc), pta)
end

function triangle:conway_circle()
	local t
	t = report_(self.pb, self.pa, length_(self.pb, self.pc), self.pa)
	return circle:new(self.incenter, t)
end

function triangle:pedal_circle(pt)
	local x, y, z, c
	x = projection_(self.pb, self.pc, pt)
	y = projection_(self.pa, self.pc, pt)
	z = projection_(self.pa, self.pb, pt)
	c = circum_center_(x, y, z)
	return circle:new(c, x)
end

function triangle:first_lemoine_circle()
	local lc, oc
	lc = self:lemoine_point()
	oc = self.circumcenter
	return circle:new(midpoint_(lc, oc), intersection_ll_(lc, ll_from_(lc, self.pa, self.pb), self.pa, self.pc))
end

function triangle:second_lemoine_circle()
	local lc, a, b, c, r, th
	lc = self:lemoine_point()
	a = point.abs(self.pc - self.pb)
	b = point.abs(self.pa - self.pc)
	c = point.abs(self.pb - self.pa)
	r = (a * b * c) / (a * a + b * b + c * c)
	th = lc + point(r, 0)
	return circle:new(lc, th)
end

function triangle:bevan_circle()
	local o, r, s, t
	o = circum_center_(excentral_tr_(self.pa, self.pb, self.pc))
	r, s, t = excentral_tr_(self.pa, self.pb, self.pc)
	return circle:new(o, r)
end

function triangle:taylor_circle()
	local a, b, c = orthic_tr_(self.pa, self.pb, self.pc)
	local d = projection_(self.pa, self.pb, a)
	local e = projection_(self.pb, self.pc, b)
	local f = projection_(self.pc, self.pa, c)
	return circle:new(circum_center_(d, e, f), d)
end

function triangle:adams_circle()
	local i = self:kimberling(1)
	local g = self:kimberling(7)
	local m = projection_(self.pb, self.pc, i)
	local n = projection_(self.pa, self.pb, i)
	local x = ll_from_(g, m, n)
	local t = intersection_ll_(x, g, self.pa, self.pb)
	return circle:new(i, t)
end

function triangle:lamoen_circle()
	local ma, mb, mc = medial_tr_(self.pa, self.pb, self.pc)
	local g = centroid_(self.pa, self.pb, self.pc)
	local p1 = circum_center_(self.pa, mb, g)
	local p2 = circum_center_(self.pb, ma, g)
	local p3 = circum_center_(self.pc, ma, g)
	local c = circum_center_(p1, p2, p3)
	return circle:new(c, p1)
end

-- soddy circles
function triangle:soddy_circle(outer)
	outer = (outer == "outer")
	local index = outer and 175 or 176
	local c = self:kimberling(index)
	local ra = (self.b + self.c - self.a) / 2
	local ta = report_(self.pa, self.pb, ra, self.pa)
	local x, y = intersection_lc_(self.pa, c, self.pa, ta)
	local t
	if outer then
		t = (length_(x, c) > length_(y, c)) and x or y
	else
		t = (length_(x, c) < length_(y, c)) and x or y
	end
	return circle:new(c, t)
end

function triangle:yiu_circles()
	local x, y, z = self:yiu_centers()
	return circle:new(x, self.pa), circle:new(y, self.pb), circle:new(z, self.pc)
end

function triangle:kenmotu_circle()
	local a, b, c = self.a, self.b, self.c
	local area = self.area
	local rho = (math.sqrt(2) * a * b * c) / (4 * area + (a ^ 2 + b ^ 2 + c ^ 2))
	local center = self:kimberling(371)
	return circle:new(through(center, rho))
end


-- Circle tangent to the sides originating from p and to the circle passing through the other two vertices.
-- argument the vertex of the chosen triangle, the center of the circle passing through the other two vertices.

function triangle:c_c(p, C)
	local a, b, c
	local i = resolve_triangle_index(self, p)
	if i == 0 then
		a, b, c = self.pa, self.pb, self.pc
	elseif i == 1 then
		a, b, c = self.pb, self.pc, self.pa
	elseif i == 2 then
		a, b, c = self.pc, self.pa, self.pb
	end
	local Lm = self:mediator(p)
	local n = intersection(Lm, C, {near = p})
	local i = self.incenter
	local q = line:new(c, i):reflection(n)
	local L = line:new(c, q):parallel_from(i)
	local Lac = line:new(a, c)
	local e = intersection(L, Lac)
	local Lia = line:new(i, a)
	local Le = Lac:orthogonal_from(e)
	return circle(intersection(Lia, Le), e)
end
triangle.thebault = triangle.c_c

function triangle:thebault_circles(p)
	local pa, pb, pc = self.pa, self.pb, self.pc
	local a, b, c    -- sommets réordonnés
	local side       -- côté sur lequel se trouve p (objet line)

	--------------------------------------------------
	-- 1. Déterminer le côté contenant p et (a,b,c)
	--    a = sommet opposé au côté portant p
	--------------------------------------------------
	if self.bc:on_segment(p) then
		-- p sur BC -> sommet opposé A
		a, b, c = pa, pb, pc
		side    = self.bc

	elseif self.ca:on_segment(p) then
		-- p sur CA -> sommet opposé B
		a, b, c = pb, pc, pa
		side    = self.ca

	elseif self.ab:on_segment(p) then
		-- p sur AB -> sommet opposé C
		a, b, c = pc, pa, pb
		side    = self.ab

	else
		tex.error("thebault_circles: point is not on any side of the triangle.")
		return
	end

	--------------------------------------------------
	-- 2. Construction des centres (comme dans ta version)
	--------------------------------------------------
	local i  = self.incenter        -- incenter du triangle ABC
	local ic = in_center_(p, a, c)  -- incenter du triangle (p,a,c)
	local ib = in_center_(p, a, b)  -- incenter du triangle (p,a,b)

	-- projections de p sur (ic i) et (ib i)
	local pc_ = projection_(p, ic, i)
	local pb_ = projection_(p, ib, i)

	-- intersections avec (b c)
	local tc = intersection_ll_(i, pc_, b, c)
	local tb = intersection_ll_(i, pb_, b, c)

	-- perpendiculaires à (b c) par tc et tb
	local xc = ortho_from_(tc, b, c)
	local xb = ortho_from_(tb, b, c)

	-- centres des deux cercles
	local center_c = intersection_ll_(p, ic, tc, xc)
	local center_b = intersection_ll_(p, ib, tb, xb)

	if not center_b or not center_c then
		tex.error("thebault_circles: failed to construct circle centers.")
		return
	end

	--------------------------------------------------
	-- 3. Points de tangence sur le côté et cercles
	--------------------------------------------------
	-- side:projection(pt1, pt2) renvoie les projections
	local tb_tan, tc_tan = side:projection(center_b, center_c)

	local Cb = circle:new(center_b, tb_tan)
	local Cc = circle:new(center_c, tc_tan)

	return Cb, Cc
end



function triangle:mixtilinear_incircle(arg)
	local i = resolve_triangle_index(self, arg)

	local incenter = self.incenter
	local pts = { self.pa, self.pb, self.pc }
	local segs = { self.ab, self.bc, self.ca }

	local A = pts[i + 1]
	local B = pts[(i + 1) % 3 + 1]
	local C = pts[(i + 2) % 3 + 1]

	local AB = segs[i + 1]
	local BC = segs[(i + 1) % 3 + 1]
	local CA = segs[(i + 2) % 3 + 1]

	local l = line:new(A, incenter)
	local lp = l:orthogonal_from(incenter)
	local P1 = intersection(lp, CA)
	local P2 = intersection(lp, AB)
	local l1 = CA:orthogonal_from(P1)
	local l2 = AB:orthogonal_from(P2)
	local W = intersection(l1, l2)

	return circle(W, P1)
end

function triangle:feuerbach_apollonius(EPS)
	EPS = EPS or tkz.epsilon
	-- 1. Points de Feuerbach
	 local Ea, Eb, Ec = self:feuerbach():get()
	-- -- 2. Excentre Ja
	local C_JaEa = self:ex_circle()
	-- 4. Spieker center
	local S = self.spiekercenter
	-- 5. Cercle orthogonal
	local Cortho = C_JaEa:orthogonal_from(S)
	-- 6. Cercle d'Euler
	local Ceuler = self:euler_circle()
	-- 7. Inversion du cercle
	local InvEuler = Cortho:inversion(Ceuler)
	local O = InvEuler.center
	-- 8. Points inversés
	local xa, xb, xc = Cortho:inversion(Ea, Eb, Ec)
  return circle:new(O, xa)
end

function triangle:feuerbach_apollonius_k181(EPS)
	EPS = EPS or tkz.epsilon
	local Ea, Eb, Ec = self:feuerbach():get() --Feuerbach points
	local A = self:kimberling(181)-- Apollonius point
	local S = self.spiekercenter
  local xa = intersection_ll_(self.pa,A,S,Ea)
	local xb = intersection_ll_(self.pb,A,S,Eb)
	local xc = intersection_ll_(self.pc,A,S,Ec)
	return circle:new(circum_circle_(xa, xb, xc))
end

function triangle:the_apollonius_circle()
	local ce = self:the_apollonius_center()
	local apo = self:kimberling(181)
	local ll = line(apo, self.pa)
	local cxa = self:ex_circle()
	local u,v = intersection(ll, cxa,{near = self.pa})
	return circle:new(ce, v)
end

-------------------
-- Result -> triangle
-------------------
function triangle:orthic()
	return triangle:new(orthic_tr_(self.pa, self.pb, self.pc))
end

function triangle:medial()
	return triangle:new(medial_tr_(self.pa, self.pb, self.pc))
end

function triangle:incentral()
	return triangle:new(incentral_tr_(self.pa, self.pb, self.pc))
end

function triangle:excentral()
	return triangle:new(excentral_tr_(self.pa, self.pb, self.pc))
end

function triangle:intouch()
	return triangle:new(intouch_tr_(self.pa, self.pb, self.pc))
end

function triangle:contact()
	return triangle:new(intouch_tr_(self.pa, self.pb, self.pc))
end
triangle.contact = triangle.intouch
triangle.gergonne = triangle.intouch

function triangle:extouch()
	return triangle:new(extouch_tr_(self.pa, self.pb, self.pc))
end
triangle.nagel = triangle.extouch

function triangle:feuerbach()
	return triangle:new(feuerbach_tr_(self.pa, self.pb, self.pc))
end

function triangle:anti()
	return triangle:new(anti_tr_(self.pa, self.pb, self.pc))
end

triangle.anticomplementary = triangle.anti
triangle.similar = triangle.anti

function triangle:tangential()
	return triangle:new(tangential_tr_(self.pa, self.pb, self.pc))
end

function triangle:cevian(p)
	return triangle:new(cevian_(self.pa, self.pb, self.pc, p))
end

function triangle:symmedian()
	local p = lemoine_point_(self.pa, self.pb, self.pc)
	return triangle:new(cevian_(self.pa, self.pb, self.pc, p))
end
triangle.symmedial = triangle.symmedian

function triangle:lemoine()
	local X = self:kimberling(598)
	return triangle:new(cevian_(self.pa, self.pb, self.pc, X))
end

function triangle:macbeath()
	local X = self:kimberling(264)
	return triangle:new(cevian_(self.pa, self.pb, self.pc, X))
end

function triangle:euler()
	return triangle:new(euler_points_(self.pa, self.pb, self.pc))
end

function triangle:pedal(pt)
	return triangle:new(self:projections(pt))
end

function triangle:yiu()
	return triangle:new(self:yiu_centers())
end


function triangle:reflection()
	return triangle:new(
		symmetry_axial_(self.pb, self.pc, self.pa),
		symmetry_axial_(self.pa, self.pc, self.pb),
		symmetry_axial_(self.pa, self.pb, self.pc)
	)
end

function triangle:projection()
	return triangle:new(
		projection_(self.pb, self.pc, self.pa),
		projection_(self.pa, self.pc, self.pb),
		projection_(self.pa, self.pb, self.pc)
	)
end

function triangle:circumcevian(pt)
	local circum = self:circum_circle()
	local x = intersection_lc_(self.pa, pt, circum.center, circum.through)
	local y = intersection_lc_(self.pb, pt, circum.center, circum.through)
	local z = intersection_lc_(self.pc, pt, circum.center, circum.through)
	return triangle:new(x, y, z)
end


function triangle:morley()
	local la1, la2 = self:trisector(self.pa)
	local lb1, lb2 = self:trisector(self.pb)
	local lc1, lc2 = self:trisector(self.pc)
	local D = intersection(la1, lb2)
	local E = intersection(lb1, lc2)
	local F = intersection(lc1, la2)
	return triangle:new(D, E, F)
end


function triangle:napoleon(swap)
	swap = (swap == "swap")
	if swap then
		local a = self.ab:equilateral().centroid
		local b = self.bc:equilateral().centroid
		local c = self.ca:equilateral().centroid
		return triangle:new(a, b, c)
	else
		local a = self.ab:equilateral("swap").centroid
		local b = self.bc:equilateral("swap").centroid
		local c = self.ca:equilateral("swap").centroid
		return triangle:new(a, b, c)
	end
end


function triangle:soddy(outer)
	local ci = self:soddy_circle(outer)
	local a, b, c = self.pa, self.pb, self.pc
	outer = (outer == "outer")

	local ce = ci.center
	local ti = ci.through
	local sa = intersection_lc_(ce, a, ce, ti)
	local sb = intersection_lc_(ce, b, ce, ti)
	local sc = intersection_lc_(ce, c, ce, ti)
	return triangle:new(sa, sb, sc)
end

-------------------
-- Result -> square
-------------------

-- Method for inscribing a square in a triangle with vertex rotation
function triangle:square_inscribed(permutation)
	-- Set a default value for rotation if not supplied
	permutation = permutation or 0
	permutation = permutation + 1
	local vertices = { self.pa, self.pb, self.pc }
	local i = (permutation - 1) % 3 + 1
	local j = i % 3 + 1
	local k = j % 3 + 1
	local a = vertices[i]
	local b = vertices[j]
	local c = vertices[k]
	return square:new(square_inscribed_(a, b, c))
end
-------------------
-- Result -> conic
-------------------
function triangle:steiner_inellipse()
	return steiner_(self.pa, self.pb, self.pc)
end

function triangle:steiner_circumellipse()
	local e = steiner_(self.pa, self.pb, self.pc)
	local v = report_(e.center, e.vertex, e.a, e.vertex)
	local cv = report_(e.center, e.covertex, e.b, e.covertex)
	return conic:new(EL_points(e.center, v, cv))
end

function triangle:euler_ellipse()
	if self:check_acutangle() then
		local x, y = intersection_lc_(self.orthocenter, self.circumcenter, self.eulercenter, self.ab.mid)
		local a = 0.5 * length_(x, y)
		return conic:new(EL_bifocal(self.orthocenter, self.circumcenter, a))
	end
end

--  kiepert hyperbola
function triangle:kiepert_hyperbola()
	-- center
	local center = self:kimberling(115)
	local O = self.circumcenter
	local G = self.centroid
	--  Brocard axis
	local L = self:brocard_axis()
	local M, N = intersection_lc_(L.pa, L.pb, O, self.pa)
	-- simson lines
	local Lsa = self:simson_line(M)
	local Lsb = self:simson_line(N)
	local axis = tkz.bisector(center, Lsa.pa, Lsb.pa)
	-- new frame system
	local minor = axis:ortho_from(center)
	local sys = occs:new(minor, center)
	local x, y = sys:coordinates(self.pa)
	local a = math.sqrt(x ^ 2 - y ^ 2)
	local c = math.sqrt(2) * a
	local Fa = report_(center, axis.pb, -c, center)
	local Fb = report_(center, axis.pb, c, center)
	local K = report_(center, axis.pb, -a ^ 2 / c, center)
	local directrix = axis:ortho_from(K)
	return conic:new(Fa, directrix, math.sqrt(2))
end

--  kiepert parabola
function triangle:kiepert_parabola()
	if (self.a == self.b) or (self.a == self.c) or (self.c == self.b) then
		tex.error("An error has occurred", { "No X(110)" })
	else
		L = self:euler_line()
		F = self:kimberling(110)
		return conic:new(F, L, 1)
	end
end

--lemoine inellipse
-- bifocal a besoin des deux foyers et d'un point sur l'ellipse ou bien de a
function triangle:lemoine_inellipse()
	local G = self.centroid
	local L = self:lemoine_point()
	local Xa = self:lemoine().pa
	return conic:new(EL_bifocal(G, L, Xa))
end

--brocard inellipse
function triangle:brocard_inellipse()
	local B1 = self:first_brocard_point()
	local B2 = self:second_brocard_point()
	local Ka = self:symmedian().pa
	return conic:new(EL_bifocal(B2, B1, Ka))
end

-- orthic inellipse
function triangle:macbeath_inellipse()
	if self:check_acutangle() then
		local O = self.circumcenter
		local H = self.orthocenter
		local Xa = self:macbeath().pa
		return conic:new(EL_bifocal(O, H, Xa))
	else
		tex.error("This function is only valid if none of the  angles is obtuse.")
	end
end

function triangle:mandart_inellipse()
	local T = self
	local M = T:mittenpunkt_point()
	local Xa, Yb, Zc = T:extouch():get()
	local X, Y = M:symmetry(Xa, Yb)
	local coefficients = search_ellipse_(X, Y, Xa, Yb, Zc)
	local center, ra, rb, angle = ellipse_axes_angle(coefficients)
	return conic:new(EL_radii(center, ra, rb, angle))
end

function triangle:orthic_inellipse()
	if self:check_acutangle() then
		local T = self
		local K = T:symmedian_point()
		local Ha, Hb, Hc = T:orthic():get()
		local Sa, Sb = K:symmetry(Ha, Hb)
		local coefficients = search_ellipse_(Ha, Hb, Hc, Sa, Sb)
		local center, ra, rb, angle = ellipse_axes_angle(coefficients)
		return conic:new(EL_radii(center, ra, rb, angle))
	else
		tex.error("This function is only valid if none of the angles is obtuse.")
	end
end

function triangle:path(p1, p2)
	local list = {}

	-- Liste des côtés et sommets
	local sides = {
		{ seg = self.ab, A = self.pa, B = self.pb },
		{ seg = self.bc, A = self.pb, B = self.pc },
		{ seg = self.ca, A = self.pc, B = self.pa },
	}

	-- Trouver les indices des segments contenant p1 et p2
	local function locate(p)
		for i, s in ipairs(sides) do
			if s.seg:on_segment(p) then
				return i
			end
		end
		tex.error("Point not on triangle boundary.")
	end

	local i1 = locate(p1)
	local i2 = locate(p2)

	-- Cas simple : même segment
	if i1 == i2 then
		table.insert(list, "(" .. checknumber_(p1.re) .. "," .. checknumber_(p1.im) .. ")")
		table.insert(list, "(" .. checknumber_(p2.re) .. "," .. checknumber_(p2.im) .. ")")
		return path:new(list)
	end

	-- Sinon, commencer par p1
	table.insert(list, "(" .. checknumber_(p1.re) .. "," .. checknumber_(p1.im) .. ")")

	-- Parcourir les sommets intermédiaires dans l’ordre trigonométrique
	local i = i1
	repeat
		local s = sides[i]
		table.insert(list, "(" .. checknumber_(s.B.re) .. "," .. checknumber_(s.B.im) .. ")")
		i = i % 3 + 1
	until i == i2

	-- Terminer par p2
	table.insert(list, "(" .. checknumber_(p2.re) .. "," .. checknumber_(p2.im) .. ")")

	return path:new(list)
end

--- CONTACT PROBLEM APOLLONIUS

-- Circle tangent to the three lines of the triangle:
-- returns (in this order): incircle, excircle at pa, excircle at pb, excircle at pc.
-- Optionnel: which ∈ { nil/"all", "in", "pa", "pb", "pc", "ex_a", "ex_b", "ex_c" }
function triangle:LLL(which)
	which = which or "all"
	-- --- sélecteur
	local C_in = self:in_circle()
	local C_ea = self:ex_circle(0) -- pa
	local C_eb = self:ex_circle(1) -- pb
	local C_ec = self:ex_circle(2) -- pc
	if which == "all" or which == nil then
		return C_in, C_ea, C_eb, C_ec
	end

	-- alias pratiques
	if which == "in" then return C_in end
	if which == "pa" or which == "ex_a" then return C_ea end
	if which == "pb" or which == "ex_b" then return C_eb end
	if which == "pc" or which == "ex_c" then return C_ec end

	tex.error("Unknown selector for triangle:c_lll: "..tostring(which))

end
triangle.c_lll = triangle.LLL



-- code for an inellipse in a triangle
local function side_of_point_(pa, pb, pc, pt, eps)
	eps = eps or tkz.epsilon

	if line_on_(pa, pb, pt, eps) then return "ab" end
	if line_on_(pb, pc, pt, eps) then return "bc" end
	if line_on_(pc, pa, pt, eps) then return "ca" end

	return nil
end

function triangle:cevian_point(p, q, eps)
	eps = eps or tkz.epsilon

	local sp = side_of_point_(self.pa, self.pb, self.pc, p, eps)
	local sq = side_of_point_(self.pa, self.pb, self.pc, q, eps)

	if not sp or not sq or sp == sq then
		tex.error("cevian_point: points must lie on two distinct sides of the triangle.")
		return nil
	end

	local br
	if (sp == "ab" and sq == "ca") or (sp == "ca" and sq == "ab") then
		br = intersection(line(self.pc, p), line(self.pb, q))
	elseif (sp == "ab" and sq == "bc") or (sp == "bc" and sq == "ab") then
		br = intersection(line(self.pc, p), line(self.pa, q))
	else
		br = intersection(line(self.pb, p), line(self.pa, q))
	end

	if br == nil then
		tex.error("cevian_point: unable to construct the cevian point.")
		return nil
	end

	local x, y, z = cevian_(self.pa, self.pb, self.pc, br)
	return br, x, y, z
end


local function ellipse_from_brianchon_(T, br, eps)
	eps = eps or tkz.epsilon

	if br == nil then
		tex.error("ellipse_inscribed: Brianchon's point is required.")
		return nil
	end

	local a, b, c = cevian_(T.pa, T.pb, T.pc, br)
	local k = T:isotomic(br)
	local o = T:complement(k)
	local u, v = o:symmetry(a, b)

	return conic:new(EL_five_points(a, b, c, u, v))
end


function triangle:ellipse_inscribed(spec, eps)
	eps = eps or tkz.epsilon

	if type(spec) ~= "table" then
		tex.error("ellipse_inscribed: a specification table is expected.")
		return nil
	end

	if spec.brianchon ~= nil then
		return ellipse_from_brianchon_(self, spec.brianchon, eps)
	end

	if spec.center ~= nil then
		local o = spec.center
		local k = self:anticomplement(o)
		local br = self:isotomic(k)
		return ellipse_from_brianchon_(self, br, eps)
	end

	if spec.focus ~= nil then
		local f1 = spec.focus
		local f2 = self:isogonal(f1)
		local o = midpoint_(f1, f2)
		local k = self:anticomplement(o)
		local br = self:isotomic(k)
		return ellipse_from_brianchon_(self, br, eps)
	end

	tex.error("ellipse_inscribed: expected brianchon, center, or focus.")
	return nil
end


return triangle