-- tkz_elements_functions_circles.lua
-- date 2025/03/04
-- version 3.34c
-- Copyright 2025  Alain Matthes
-- This work may be distributed and/or modified under the
-- conditions of the LaTeX Project Public License, either version 1.3
-- of this license or (at your option) any later version.
-- The latest version of this license is in
-- http://www.latex-project.org/lppl.txt
-- and version 1.3 or later is part of all distributions of LaTeX
-- version 2005/12/01 or later.
-- This work has the LPPL maintenance status “maintained”.
-- The Current Maintainer of this work is Alain Matthes.

--  define a circle by the center and a radius
function through(c,r,an)
  if an then   
    return c, c + point(r * math.cos(an), r * math.sin(an))
  else
    return c, c + point(r, 0)
  end
end

function diameter(a, b, swap)
	if swap then
		return midpoint_(a, b), a
	else
		return midpoint_(a, b), b
	end
end

function midarc_(o, a, b) -- a -> b
	local phi = 0.5 * get_angle_(a, o, b)
	return rotation_(o, phi, b)
end

function tangent_from_(c, p, pt)
	local o
	o = midpoint_(c, pt)
	return intersection_cc_(o, c, c, p)
end

function tangent_at_(a, p)
	return rotation_(p, math.pi / 2, a), rotation_(p, -math.pi / 2, a)
end

function orthogonal_from_(a, b, p)
	return tangent_from_(a, b, p)
end

function orthogonal_through_(a, b, x, y)
	local d = 1 / point.mod(x - a)
	local z = a + (b - a) * d
	return circum_center_(x, y, z)
end

function inversion_(c, p, pt)
	local ry = point.abs(c - p)
	local d = point.abs(c - pt)
	local r = (ry * ry) / d
	return c + polar_(r, point.arg(pt - c))
end

function circles_position_(c1, r1, c2, r2)
	local d = point.mod(c1 - c2)
	local max = r1 + r2
	local min = math.abs(r1 - r2)

	if d > max then
		return "outside"
	elseif math.abs(d - max) < tkz_epsilon then
		return "outside tangent"
	elseif math.abs(d - min) < tkz_epsilon then
		return "inside tangent"
	elseif d < min then
		return "inside"
	else
		return "intersect"
	end
end

function radical_axis_(c1, p1, c2, p2)
	local r1 = point.abs(c1 - p1)
	local r2 = point.abs(c2 - p2)
	local d = point.abs(c1 - c2)
	local h = (r1 * r1 - r2 * r2 + d * d) / (2 * d)

	local ck = radical_center_(c1, p1, c2, p2)
	local cj = rotation_(ck, -math.pi / 2, c1)
	local ci = symmetry_(ck, cj)

	return cj, ci
end

function radical_center_(c1, p1, c2, p2)
	local d, r1, r2, h
	r1 = point.abs(c1 - p1)
	r2 = point.abs(c2 - p2)
	d = point.abs(c1 - c2)
	h = (r1 * r1 - r2 * r2 + d * d) / (2 * d)
	return h * (c2 - c1) / d + c1
end

function radical_center3(C1, C2, C3)
	local t1, t2, t3, t4
	t1, t2 = radical_axis_(C1.center, C1.through, C2.center, C2.through)
	if C3 == nil then
		return intersection_ll_(t1, t2, C1.center, C2.center)
	else
		t3, t4 = radical_axis_(C3.center, C3.through, C2.center, C2.through)
		return intersection_ll_(t1, t2, t3, t4)
	end
end

function south_pole_(c, p)
	return c - point(0, point.abs(c - p))
end

function north_pole_(c, p)
	return c + point(0, point.abs(c - p))
end

function antipode_(c, pt)
	return 2 * c - pt
end

function internal_similitude_(c1, r1, c2, r2)
	return barycenter_({ c2, r1 }, { c1, r2 })
end

function external_similitude_(c1, r1, c2, r2)
	return barycenter_({ c2, r1 }, { c1, -r2 })
end

function circlepoint_(c, t, k)
	local phi = 2 * k * math.pi
	return rotation_(c, phi, t)
end

function midcircle_(C1, C2)
	local state, r, s, t1, t2, T1, T2, p, a, b, c, d, Cx, Cy, i, j
	state = circles_position_(C1.center, C1.radius, C2.center, C2.radius)
	i = barycenter_({ C2.center, C1.radius }, { C1.center, -C2.radius })
	j = barycenter_({ C2.center, C1.radius }, { C1.center, C2.radius })
	t1, t2 = tangent_from_(C1.center, C1.through, i)
	T1, T2 = tangent_from_(C2.center, C2.through, i)

	if (state == "outside") or (state == "outside tangent") then
		p = math.sqrt(point.mod(i - t1) * point.mod(i - T1))
		return circle:radius(i, p)
	elseif state == "intersect" then
		r, s = intersection(C1, C2)
		return circle:radius(i, point.mod(r - i)), circle:radius(j, point.mod(r - j))
	elseif (state == "inside") or (state == "inside tangent") then
		a, b = intersection_lc_(C1.center, C2.center, C1.center, C1.through)
		c, d = intersection_lc_(C1.center, C2.center, C2.center, C2.through)

		-- Ensure the smaller radius circle is used first
		if C1.radius < C2.radius then
			z.u, z.v, z.r, z.s = a, b, c, d
		else
			z.u, z.v, z.r, z.s = c, d, a, b
		end

		-- Determine circle orientation and return orthogonal from j
		if in_segment_(z.s, z.v, z.u) then
			Cx = circle:diameter(z.r, z.v)
			Cy = circle:diameter(z.u, z.s)
		else
			Cx = circle:diameter(z.s, z.v)
			Cy = circle:diameter(z.u, z.r)
		end

		-- Return the circle with the smaller radius orthogonal from j
		if Cx.radius < Cy.radius then
			return Cy:orthogonal_from(j)
		else
			return Cx:orthogonal_from(j)
		end
	end
end