-- tkz_elements_intersections.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.

-------------------------------------------------------------------------
-- intersection of lines
-------------------------------------------------------------------------
function intersection_ll (la,lb)
    return intersection_ll_ (la.pa,la.pb,lb.pa,lb.pb)
end
---------------------------------------------------------------------------
-- intersection of a line and a circle
---------------------------------------------------------------------------
function intersection_lc (D,C )
    return intersection_lc_ ( D.pa,D.pb ,C.center,C.through )
end -- function
---------------------------------------------------------------------------
-- intersection of two circles
---------------------------------------------------------------------------
function intersection_cc (Ca , Cb )
   return intersection_cc_(Ca.center,Ca.through,Cb.center,Cb.through)
end -- function

--  line ellipse
function intersection_le (L,E)
   local a,b,c,d,t1,t2,z1,z2,A,B,Bx,By,Ax,Ay,Rx,Ry,sd
   A = (L.pa - E.center)*(point(math.cos(E.slope),-math.sin(E.slope)))
   B = (L.pb - E.center)*(point(math.cos(E.slope),-math.sin(E.slope)))
   Rx  = E.Rx
   Ry  = E.Ry
   Ax  = A.re
   Ay  = A.im
   Bx  = B.re
   By  = B.im
    a  = Rx^2  * (By-Ay)^2 +Ry^2 * (Bx-Ax)^2
    b  = 2 * Rx^2 * Ay * (By-Ay)+ 2  * Ry^2 * Ax * (Bx-Ax)
    c  = Rx^2 * Ay^2 + Ry^2 * Ax^2 -  Rx^2 * Ry^2
    d  =  b^2 - 4 * a * c

   if d > 0 then
      sd = math.sqrt(d)
      t1 = (-(b)+sd)/(2*a)
      t2 = (-(b)-sd)/(2*a)
      z1 = point ( Ax + (Bx-Ax)*t1 ,  Ay + (By-Ay)*t1 )
      z2 = point ( Ax + (Bx-Ax)*t2 ,  Ay + (By-Ay)*t2 )
      if angle_normalize (point.arg(z1)) < angle_normalize (point.arg(z2))
      then
      return z1*(point(math.cos(E.slope),math.sin(E.slope))) + E.center,
             z2*(point(math.cos(E.slope),math.sin(E.slope))) + E.center
          else
             return z2*(point(math.cos(E.slope),math.sin(E.slope))) + E.center,
                    z1*(point(math.cos(E.slope),math.sin(E.slope))) + E.center
            end -- if 
       elseif  math.abs(d) < tkz_epsilon 
       then
           t1 = (-(b))/(2*a)
           z1 = point ( Ax + (Bx-Ax)*t1 ,  Ay + (By-Ay)*t1 )           
          return z1*(point(math.cos(E.slope),math.sin(E.slope))) + E.center,
                 z1*(point(math.cos(E.slope),math.sin(E.slope))) + E.center
       else
         return false,false
      end
end
 
function intersection_ll_(a, b, c, d)
    local x1, y1, x2, y2, x3, y3, x4, y4
    local DN, NX, NY

    x1, y1 = a.re, a.im
    x2, y2 = b.re, b.im
    x3, y3 = c.re, c.im
    x4, y4 = d.re, d.im

    DN = (x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4)

    if math.abs(DN) < tkz_epsilon then
        return false
    end

    NX = (x1 * y2 - y1 * x2) * (x3 - x4) - (x1 - x2) * (x3 * y4 - y3 * x4)
    NY = (x1 * y2 - y1 * x2) * (y3 - y4) - (y1 - y2) * (x3 * y4 - y3 * x4)

    return point(NX / DN, NY / DN)
end

function intersection_lc_(pa, pb, c, p)
    local zh, dh, arg_ab, phi, c1, c2, r, test

    r = point.mod(c - p)
    zh = projection_(pa, pb, c)
    dh = point.abs(c - zh)
    arg_ab = point.arg(pa - pb)

    if dh < tkz_epsilon then
        -- Le centre du cercle est sur la droite
        return c + polar_(r, math.pi + arg_ab), c + polar_(r, arg_ab)
    elseif math.abs(r - dh) < tkz_epsilon then
        -- La droite est tangente au cercle
        return zh, zh
    elseif dh > r then
        -- Aucune intersection
        return false, false
    else
        -- Il y a une intersection, calcul de l'angle
        phi = math.asin(dh / r)
        test = (pa - pb) * point.conj(c - zh)
        if test.im < 0 then
            phi = math.pi + phi
        end

        c1 = angle_normalize(arg_ab + phi)
        c2 = angle_normalize(math.pi + arg_ab - phi)

        -- Retourner les deux points d'intersection
        if c2 < c1 then
            return c + polar_(r, c2), c + polar_(r, c1)
        else
            return c + polar_(r, c1), c + polar_(r, c2)
        end
    end
end

function intersection_cc_(ca, pa, cb, pb)
    local d, cosphi, phi, ra, rb, c1, c2, epsilon

    -- Précision pour arrondir les résultats
    epsilon = 12
    -- Distance entre les centres des cercles
    d = point.abs(ca - cb)
    -- Rayons des cercles
    ra = point.abs(ca - pa)
    rb = point.abs(cb - pb)

    -- Calcul du cosinus de l'angle phi entre les centres et les points sur les cercles
    cosphi = tkzround((ra * ra + d * d - rb * rb) / (2 * ra * d), epsilon)
    
    -- Calcul de l'angle phi
    phi = tkzround(math.acos(cosphi), epsilon)

    -- Si phi est invalide (par exemple, cosphi > 1 ou < -1), aucune intersection
    if not phi then
        return false, false
    elseif phi == 0 then
        -- Les cercles sont tangents l'un à l'autre, retourne le même point pour les deux intersections
        return ca + polar_(ra, point.arg(cb - ca)), ca + polar_(ra, point.arg(cb - ca))
    else
        -- Calcul des angles des points d'intersection
        c1 = angle_normalize(phi + point.arg(cb - ca))
        c2 = angle_normalize(-phi + point.arg(cb - ca))

        -- Retourner les points d'intersection dans l'ordre croissant des angles
        if c1 < c2 then
            return ca + polar_(ra, c1), ca + polar_(ra, c2)
        else
            return ca + polar_(ra, c2), ca + polar_(ra, c1)
        end
    end
end

function intersection(X, Y)
   local t = {}  -- Tableau pour stocker les points d'intersection

   -- Cas où X est un cercle
   if X.type == 'circle' then
      if Y.type == 'circle' then
         -- Intersection entre deux cercles
         local z1, z2 = intersection_cc(X, Y)
         table.insert(t, z1)
         table.insert(t, z2)
      else -- Y est une droite
         local z1, z2 = intersection_lc(Y, X)  -- Intersection entre un cercle et une droite
         table.insert(t, z1)
         table.insert(t, z2)
      end

   -- Cas où X est une droite
   elseif X.type == 'line' then
      if Y.type == 'circle' then
         -- Intersection entre une droite et un cercle
         local z1, z2 = intersection_lc(X, Y)
         table.insert(t, z1)
         table.insert(t, z2)
      elseif Y.type == 'line' then
         -- Intersection entre deux droites
         local z1 = intersection_ll(X, Y)
         table.insert(t, z1)
      else -- Y est une ellipse
         local z1, z2 = intersection_le(X, Y)  -- Intersection entre une droite et une ellipse
         table.insert(t, z1)
         table.insert(t, z2)
      end

   -- Cas où X est une ellipse
   elseif X.type == 'ellipse' then
      -- Intersection entre une ellipse et l'autre objet (cercle, droite ou ellipse)
      local z1, z2 = intersection_le(Y, X)
      table.insert(t, z1)
      table.insert(t, z2)
   end

   -- Retourner les résultats sous forme de valeurs distinctes
   return table.unpack(t)
end