-- 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
		error("tkz-elements: triangle: degenerate triangle (two vertices are identical)", 2)
	end
	-- aligned vertices
	if tkz.is_linear(za, zb, zc) then
		error("tkz-elements: triangle: degenerate triangle (the vertices are collinear)", 2)
	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" -- border zone (tolerance)
	end
end


-- strict = true: the boundary is excluded
function triangle:in_out(pt, strict)
	if strict == true then
		return in_out_strict_(self.pa, self.pb, self.pc, pt)
	end
	return in_out_(self.pa, self.pb, self.pc, pt)
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

-- get_angle(arg): measure (radians) of the angle at a vertex;
-- arg = 0, 1, 2 (vertex A, B, C) or a vertex given as a point.
function triangle:get_angle(arg)
	local i = resolve_triangle_index_(self, arg)  -- 0, 1, 2
	if i == 0 then return self.alpha end  -- angle at A
	if i == 1 then return self.beta  end  -- angle at B
	return self.gamma                     -- angle at C
end

-- angle(i): the angle object (class angle) at a vertex;
-- i = 1, 2, 3 (as for get(i)) or a vertex given as a point.
function triangle:angle(i)
	local angles = { self.alpha_, self.beta_, self.gamma_ }
	if type(i) == "number" then
		if i ~= 1 and i ~= 2 and i ~= 3 then
			error("tkz-elements: triangle:angle: the index must be 1, 2 or 3, got " .. tostring(i), 2)
		end
		return angles[i]
	elseif type(i) == "table" and i.type == "point" then
		for k, v in ipairs({ self.pa, self.pb, self.pc }) do
			if point.abs(i - v) <= tkz.epsilon then
				return angles[k]
			end
		end
		error("tkz-elements: triangle:angle: the point is not a vertex of the triangle", 2)
	end
	error("tkz-elements: triangle:angle: expected 1, 2, 3 or a vertex, got " .. tostring(i), 2)
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
		error("tkz-elements: triangle:parameter: the point is not on the boundary of the triangle", 2)
	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 -- perimeter
	local t1 = self.c / p -- end of AB
	local t2 = (self.c + self.a) / p -- end of 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
		error("tkz-elements: triangle:kimberling: X(" .. tostring(n) .. ") is not available", 2)
	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
		error("tkz-elements: triangle:excenter: the argument must be a vertex of the triangle", 2)
	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 p, err = poncelet_point_(self.pa, self.pb, self.pc, pt)
	if not p then
		error("tkz-elements: triangle:poncelet_point: " .. err, 2)
	end
	return p
end

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

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

	local p, q = l:get()
	if point.abs(p - q) <= EPS then
		error("tkz-elements: triangle:orthopole: degenerate line", 2)
	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)

	-- feet on the sides CA and BC
	local a1, a2 = self.bc:get()   -- line (BC)
	local c1, c2 = self.ca:get()   -- line (CA)

	if point.abs(a2 - a1) <= EPS or point.abs(c2 - c1) <= EPS then
		error("tkz-elements: triangle:orthopole: degenerate side of the triangle", 2)
	end

	local app = projection_(a1, a2, ap)  -- foot of ap on (BC)
	local bpp = projection_(c1, c2, bp)  -- foot of bp on (CA)

	-- orthopole = intersection of the lines (ap, app) and (bp, bpp)
	return intersection(line:new(ap, app), line:new(bp, bpp))
end



function triangle:isodynamic_points()
	-- Special case: equilateral triangle
	if self:is_equilateral() then
		-- The centroid G = (A + B + C)/3 is used
		local A, B, C = self.pa, self.pb, self.pc
		local G = (A + B + C) / 3
		-- The same point is returned twice, to keep the same signature
		return G, G
	end

	-- General case
	local a, b, c = self.a, self.b, self.c
	local ab, bc = self.ab, self.bc, self.ca  -- self.ca is not used here

	-- Apollonius circles for the ratios b/a on AB and c/b on BC
	local Cab = ab:apollonius(b / a)
	local Cbc = bc:apollonius(c / b)

	-- Their intersection gives the two isodynamic points
	return intersection(Cab, Cbc)
end

triangle.neuberg_points = triangle.isodynamic_points

-- Points of intersection of the Apollonius circle (b/a, side AB)
-- with the line (AB), the point on the segment [AB] first.

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
			error("tkz-elements: triangle:apollonius_points: invalid side " .. tostring(side)
				.. " (use \"ab\", \"bc\" or \"ca\")", 2)
	end

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

     return L:harmonic(k)
end


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] -- starting vertex
	local B = S[(i + 1) % 3 + 1] -- premier sommet adjacent
	local C = S[(i + 2) % 3 + 1] -- second sommet adjacent

	-- u, v are the reflections of B and C about the bisector at 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
		error("tkz-elements: triangle:euler_line: no Euler line for an equilateral triangle", 2)
	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() -- to be checked
	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
				error("tkz-elements: triangle:apollonius_circle: invalid side " .. tostring(side)
					.. " (use \"ab\", \"bc\" or \"ca\")", 2)
		end

		-- case k = 1 (approximately): handled on purpose
		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:projections(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)          -- excenter opposite A
	local T = projection_(B, C, O)         -- point of tangency on (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    -- reordered vertices
	local side       -- side containing p (line object)

	--------------------------------------------------
	-- 1. Side containing p, and (a,b,c)
	--    a = vertex opposite the side containing p
	--------------------------------------------------
	if self.bc:on_segment(p) then
		-- p on BC -> opposite vertex A
		a, b, c = pa, pb, pc
		side    = self.bc

	elseif self.ca:on_segment(p) then
		-- p on CA -> opposite vertex B
		a, b, c = pb, pc, pa
		side    = self.ca

	elseif self.ab:on_segment(p) then
		-- p on AB -> opposite vertex C
		a, b, c = pc, pa, pb
		side    = self.ab

	else
		error("tkz-elements: triangle:thebault_circles: the point is not on a side of the triangle", 2)
	end

	--------------------------------------------------
	-- 2. Construction of the centers
	--------------------------------------------------
	local i  = self.incenter        -- incenter of the triangle ABC
	local ic = in_center_(p, a, c)  -- incenter of the triangle (p,a,c)
	local ib = in_center_(p, a, b)  -- incenter of the triangle (p,a,b)

	-- projections of p on (ic i) and (ib i)
	local pc_ = projection_(p, ic, i)
	local pb_ = projection_(p, ib, i)

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

	-- perpendiculars to (b c) through tc and tb
	local xc = ortho_from_(tc, b, c)
	local xb = ortho_from_(tb, b, c)

	-- centers of the two circles
	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
		error("tkz-elements: triangle:thebault_circles: the centers of the circles cannot be constructed", 2)
	end

	--------------------------------------------------
	-- 3. Points of tangency on the side, and circles
	--------------------------------------------------
	-- side:projection(pt1, pt2) returns the 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. Feuerbach points
	 local Ea, Eb, Ec = self:feuerbach():get()
	-- -- 2. Excenter Ja
	local C_JaEa = self:ex_circle()
	-- 4. Spieker center
	local S = self.spiekercenter
	-- 5. Orthogonal circle
	local Cortho = C_JaEa:orthogonal_from(S)
	-- 6. Euler circle
	local Ceuler = self:euler_circle()
	-- 7. Inversion of the circle
	local InvEuler = Cortho:inversion(Ceuler)
	local O = InvEuler.center
	-- 8. Inverted points
	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
		error("tkz-elements: triangle:kiepert_parabola: undefined for an isosceles triangle (no X(110))", 2)
	else
		local L = self:euler_line()
		local F = self:kimberling(110)
		return conic:new(F, L, 1)
	end
end

--lemoine inellipse
-- bifocal needs the two foci and a point on the ellipse, or 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
		error("tkz-elements: triangle:macbeath_inellipse: the triangle must not have an obtuse angle", 2)
	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
		error("tkz-elements: triangle:orthic_inellipse: the triangle must not have an obtuse angle", 2)
	end
end

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

	-- List of the sides and vertices
	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 },
	}

	-- Indices of the segments containing p1 and p2
	local function locate(p)
		for i, s in ipairs(sides) do
			if s.seg:on_segment(p) then
				return i
			end
		end
		error("tkz-elements: triangle:path: the point is not on the boundary of the triangle", 3)
	end

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

	-- Simple case: same 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) .. ")")

	-- Intermediate vertices, counterclockwise
	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.
-- Optional: which in { nil/"all", "in", "pa", "pb", "pc", "ex_a", "ex_b", "ex_c" }
function triangle:LLL(which)
	which = which or "all"
	-- --- selector
	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

	error("tkz-elements: triangle:LLL: unknown selector " .. tostring(which)
		.. " (use \"all\", \"in\", \"pa\", \"pb\" or \"pc\")", 2)

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
		error("tkz-elements: triangle:cevian_point: the points must lie on two distinct sides of the triangle", 2)
	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
		error("tkz-elements: triangle:cevian_point: the cevian point cannot be constructed", 2)
	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
		error("tkz-elements: triangle:ellipse_inscribed: the Brianchon point is required", 2)
	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
		error("tkz-elements: triangle:ellipse_inscribed: a specification table is expected", 2)
	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

	error("tkz-elements: triangle:ellipse_inscribed: the table must contain brianchon, center or focus", 2)
end


return triangle