-- tkz_elements_functions_circles.lua
-- date 2025/01/06
-- version 3.10
-- Copyright 2024  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 circle_cr ( c, r )
--     return c +  point(r,0)
-- 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