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

-- -------------------------------------------------------------------------
--                           Lines
-- -------------------------------------------------------------------------
line = {}
function line:new(za, zb)
    local type             = 'line'
    local mid              = (za+zb)/2
    local north_pa         = rotation_ (za,math.pi/2,zb)
    local south_pa         = rotation_ (za,-math.pi/2,zb)
    local north_pb         = rotation_ (zb,-math.pi/2,za)
    local south_pb         = rotation_ (zb,math.pi/2,za)
    local west             = rotation_ (za,math.pi/2,north_pa)
    local east             = rotation_ (zb,math.pi/2,south_pb)
    local slope            = angle_normalize_(point.arg(zb-za))
    local length           = point.mod(zb-za)
    local vec              = vector : new (za,zb)
    local o = {pa          = za, 
               pb          = zb,
               north_pa    = north_pa, 
               south_pa    = south_pa,
               west        = west, 
               east        = east,
               north_pb    = north_pb, 
               south_pb    = south_pb,
               slope       = slope,
               mid         = mid, 
               type        = type, 
               vec         = vec,
               length      = length}
    setmetatable(o, self)
    self.__index = self
    return o
end

-------------------
-- Result -> real
-------------------
function line:distance(pt)   
    return point.mod(projection(self,pt)-pt)
end

function line:slope()
   return slope_(self.pa,self.pb)
end
-------------------
-- Result -> boolean
-------------------
function line:in_out(pt)
    return math.abs((pt - self.pa) ^ (pt - self.pb)) <= tkz_epsilon
end

function line:in_out_segment(pt)
    return point.mod(pt - self.pa) + point.mod(pt - self.pb) - point.mod(self.pb - self.pa) <= tkz_epsilon
end

function line:is_parallel(L)
    return math.abs(self.slope - L.slope) < tkz_epsilon
end

function line:is_orthogonal(L)
    return math.abs(self.slope * L.slope + 1) < tkz_epsilon
end

function line:is_equidistant(p)
    return math.abs( (point.mod(self.pa-p)-(point.mod(self.pb-p)))) < tkz_epsilon
end
-------------------
-- Result -> point
-------------------
function line:barycenter(ka,kb)
    return barycenter_({self.pa,ka},{self.pb,kb})
end

function line:point(t) --  t=o A  t=1 B  t = AM / AB 
    return barycenter_({self.pa,1-t},{self.pb,(t)})
end

function line:midpoint()
    return (self.pa+self.pb)/2
end

function line:harmonic_int(pt)
    return div_harmonic_int_(self.pa,self.pb,pt)
end

function line:harmonic_ext(pt)
    return div_harmonic_ext_(self.pa,self.pb,pt)
end

function line:harmonic_both(k)
    return div_harmonic_both_(self.pa,self.pb,k)
end

function line:gold_ratio()
   return self.pa + (self.pb-self.pa)*tkzinvphi
end

function line:normalize()
   return  self.pa+(self.pb-self.pa)/point.mod(self.pb-self.pa)
end

function line:normalize_inv()
   return normalize_ (self.pb,self.pa)
end

function line:_east(d)
   local d = d or 1
   return self.pb+ d/self.length * (self.pb-self.pa)
end

function line:_west(d)
   local d = d or 1
   return self.pa+ d/self.length * (self.pa-self.pb)
end

function line:_north_pa(d)
   local d = d or 1
   return d/self.length * ( self.north_pa - self.pa ) + self.pa
end

function line:_south_pa(d)
   local d = d or 1
   return d/self.length *( self.south_pa - self.pa ) + self.pa
end

function line:_south_pb(d)
   local d = d or 1
   return d/self.length *( self.south_pb - self.pb ) + self.pb
end

function line:_north_pb(d)
   local d = d or 1
   return d/self.length *( self.north_pb - self.pb ) + self.pb
end

function line:report(d, pt)
    if not self.length or self.length == 0 then
        error("self.length must be non-zero")
    end
    local t = d / self.length
    local result = barycenter_({self.pa, 1 - t}, {self.pb, t})
    if pt then
        return result + pt - self.pa
    else
        return result
    end
end


function line:colinear_at (pt,k)
   if k == nil  
      then  
         return colinear_at_ (self.pa,self.pb,pt,1)
      else 
         return colinear_at_ (self.pa,self.pb,pt,k)
      end
end
-------------- transformations -------------
function line:translation_pt( pt )
    return translation_ ( self.pb-self.pa,pt )
end

function line:translation_C( obj )
   local pa,pb,x,y
   pa = obj.center
   pb = obj.through
   x,y = set_translation_( self.pb-self.pa,pa,pb )
   return circle : new  (x,y)
end

function line: translation_T( obj )
   local pa,pb,pc,x,y,z
   pa    = obj.pa
   pb    = obj.pb
   pc    = obj.pc
   x,y,z = set_translation_( self.pb-self.pa,pa,pb,pc )
   return triangle : new  (x,y,z)
end

function line: translation_L( obj )
   local pa,pb,x,y
   pa = obj.pa
   pb = obj.pb
   x,y = set_translation_ ( self.pb-self.pa,pa,pb )
   return line : new  (x,y)
end

function line:translation(...)
   local obj, nb, t
   local tp = table.pack(...)
   obj = tp[1]
   nb = tp.n

   -- If only one object is passed
   if nb == 1 then
      if obj.type == "point" then
         return translation_(self.pb - self.pa, obj)  -- Translate point
      elseif obj.type == "line" then
         return self:translation_L(obj)  -- Translate line
      elseif obj.type == "triangle" then
         return self:translation_T(obj)  -- Translate triangle
      elseif obj.type == "circle" then
         return self:translation_C(obj)  -- Translate circle
      else
         error("Unsupported object type for translation")
      end
   else
      -- If multiple objects are passed, translate each one
      t = {}
      for i = 1, nb do
         -- Translate each object using the translation vector
         table.insert(t, translation_(self.pb - self.pa, tp[i]))
      end
      return table.unpack(t)  -- Return the translated objects
   end
end


function line: set_translation ( ...)
    return set_translation_ ( self.pb-self.pa,... )
end

function line:projection(...)
   local obj,nb,t
   local tp = table.pack(...)
   obj = tp[1]
   nb = tp.n
    if nb == 1 then
          return projection_ ( self.pa, self.pb, obj )
    else
        t = {}
        for i=1,tp.n do
            table.insert( t , projection_ (self.pa, self.pb, tp[i])  ) 
         end
      return table.unpack ( t )      
    end
end

function line:set_projection(...)
	local tp = table.pack(...)
	local i
    local t = {}
	for i=1,tp.n do
        table.insert( t , projection_ (self.pa,self.pb , tp[i])  ) 
	end
  return table.unpack ( t )
end

function line:symmetry_axial_L( obj )
   local pa,pb,x,y
   pa = obj.pa
   pb = obj.pb
   x,y = self:set_reflection(pa, pb)
   return line : new  (x,y)
end
function line:symmetry_axial_T( obj )
   local pa,pb,pc,x,y,z
   pa    = obj.pa
   pb    = obj.pb
   pc    = obj.pc
   x,y,z = self:set_reflection (pa, pb, pc)
   return triangle : new  (x,y,z)
end

function line:symmetry_axial_C( obj )
   local pa,pb,x,y
   pa = obj.center
   pb = obj.through
   x,y = self:set_reflection( pa,pb )
   return circle : new  (x,y)
end

function line:reflection(...)
   local obj,nb,t
   local tp = table.pack(...)
   obj = tp[1]
   nb = tp.n
    if nb == 1 then
       if obj.type == "point" then
          return symmetry_axial_ ( self.pa,self.pb,obj )
       elseif  obj.type == "line" then
          return self: symmetry_axial_L (obj)
       elseif obj.type == "triangle" then
          return self: symmetry_axial_T (obj)
       else
         return self: symmetry_axial_C (obj)
       end
    else
        t = {}
        for i=1,tp.n do
            table.insert( t , symmetry_axial_ ( self.pa,self.pb , tp[i])  ) 
         end
      return table.unpack ( t )      
    end
end

function line:set_reflection (...)
    return set_symmetry_axial_ ( self.pb,self.pa,... )
end

-------------------
-- Result -> line
-------------------
function line:ll_from( pt )
	return line : new (pt,pt+self.pb-self.pa) 
end

function line:ortho_from( pt )
	return  line : new (pt+(self.pb-self.pa)*point(0,-1),pt+(self.pb-self.pa)*point(0,1))
end

function line:mediator() 
   local m = midpoint_ (self.pa,self.pb)
  return line : new (rotation_ (m,-math.pi/2,self.pb),rotation_ (m,math.pi/2,self.pb)) 
end

function line:perpendicular_bisector ()
   local  m = midpoint_ (self.pa,self.pb)
  return line : new (rotation_ (m,-math.pi/2,self.pb),rotation_ (m,math.pi/2,self.pb))
end

-------------------
-- Result -> circle
-------------------
function line:circle(swap)
  swap = swap or false
  if swap then
    return circle:new(self.pb,self.pa)
  else
    return circle:new(self.pa,self.pb)
  end
end

function line:circle_swap()   
    return circle : new (self.pb,self.pa)
end

function line:diameter()
   local c = midpoint_(self.pa,self.pb)
  return circle:new (c,self.pb)
end

function line:apollonius(k)
   local z1,z2,c
    z1     = barycenter_({self.pa,1},{self.pb,k})
    z2     = barycenter_({self.pa,1},{self.pb,-k})
    c      = midpoint_ (z1,z2)
  return circle : new  (c,z2)
end

function line:test(x,y) 
  
end


-- Circle tangent to a line passing through two points
-- In general, there are two solutions
function line:c_l_pp(a, b) -- a and b on the same side
  -- Initialisation
  local lab = line:new(a, b) -- Line through a and b
  local Cab = circle:diameter(a, b) -- Circle with a and b diameters
  local i = intersection(lab, self) -- Intersection with current line
 
  -- One point on the line  (a)
   if self : in_out(a)  and not self:in_out(b) then
    local lmed  = lab : mediator()
     local laperp = self:ortho_from(a)
     local o = intersection(lmed,laperp)
     return circle:new(o,a),
            circle:new(o,a)
      end
      -- One point on the line  (b)
       if self:in_out(b) and not self:in_out(a) then
         local lmed  = lab:mediator()
         local laperp = self:ortho_from(b)
         local o = intersection(lmed,laperp)
         return circle:new(o,b),
                circle:new(o,b)
       end    
  -- Check: if the intersection exists and lies on the segment [a, b].
   if i  and lab:in_out_segment(i) then
    return nil, nil -- No circle is possible
   end

  -- If the current line is orthogonal to lab
  if self:is_orthogonal(lab) then
    local lmed = lab:mediator() 
    local m = midpoint(a, b) 
    local r = length(m, i) 
    local pt1 = lab:isosceles_s(r)
    local pt2 = lab:isosceles_s(r, true)
    return circle:new(pt1, a),
           circle:new(pt2, a)
  end

  -- If the two lines are parallel
  if lab:is_parallel(self) then
    local mid = midpoint(a, b) -- Midpoint of segment [a, b]
    local proj = self:projection(mid) -- Mid projection on the running line

    return circle:new(circum_center_(a, b, proj), proj),
           circle:new(circum_center_(a, b, proj), proj)
  end

  -- General case
  local t = Cab:tangent_from(i).pb
  local x, y = intersection(self, circle:new(i, t))
  return circle:new(intersection(self:ortho_from(x), lab:mediator()), x),
         circle:new(intersection(self:ortho_from(y), lab:mediator()), y)
end


-- Circle tangent to two straight lines passing through a given point
function line:c_ll_p(a, p)
  
  -- Compute the bisector of the triangle formed by self.pa, self.pb, and a
  local lbi = bisector(self.pa, self.pb, a)

  if lbi:in_out(p) then
    -- Orthogonal projection of p onto the bisector
    local lp = lbi:ortho_from(p)

    -- Intersection of line from p to its projection with self.pa and self.pb
    local i = intersection_ll_(p, lp.pb, self.pa, self.pb)

    -- Intersection points of the line with the circle defined by (i, p)
    local t1, t2 = intersection_lc_(self.pa, self.pb, i, p)

    -- Create the main line and find orthogonal projections from t1 and t2
    local lab = line:new(self.pa, self.pb)
    local x = lab:ortho_from(t1).pb
    local y = lab:ortho_from(t2).pb

    -- Return two circles based on the orthogonal projections and points t1, t2
    return circle:new(intersection_ll_(x, t1, self.pa, p), t1),
           circle:new(intersection_ll_(y, t2, self.pa, p), t2)
  else
    -- Reflection of p across the bisector
    local q = lbi : reflection (p)

    -- Compute circles from the Wallis construction
    local c1, c2 = self:c_l_pp(p, q)

    -- Return two circles with centers and points on their circumference
    return  c1,c2
  end
end



----------------------
-- Result -> triangle
----------------------
function line:equilateral(swap)
    swap = swap or false
    if swap then
        return triangle:new(self.pa, self.pb, rotation_(self.pa, -math.pi / 3, self.pb))
    else
        return triangle:new(self.pa, self.pb, rotation_(self.pa, math.pi / 3, self.pb))
    end
end

function line:isosceles(phi,swap)
  local pta,ptb
   swap = swap or false
  if swap  then 
     pta = rotation_(self.pa,-phi,self.pb)
     ptb = rotation_(self.pb,phi,self.pa)
   return triangle : new (self.pa,self.pb, intersection_ll_(self.pa,pta,self.pb,ptb ))
  else
    pta = rotation_(self.pa,phi,self.pb)
    ptb = rotation_(self.pb,-phi,self.pa)
  return triangle : new (self.pa,self.pb, intersection_ll_ (self.pa,pta,self.pb,ptb ))
end
end
line.isosceles_a = line.isosceles

function line:isosceles_s(a,swap)
     local c1,c2,pta,ptb,pt1,pt2
  c1 = circle : radius (self.pa,a)
  c2 = circle : radius (self.pb,a)
  pta,ptb = intersection_cc (c1,c2)
   if get_angle(self.pa,self.pb,pta) < get_angle(self.pa,self.pb,ptb) then
     pt1=pta pt2=ptb
      else  pt1=ptb pt2=pta end  
   swap = swap or false
  if swap  then
   return triangle : new (self.pa,self.pb,pt2 )
  else
  return triangle : new (self.pa,self.pb,pt1)
end
end

function line:two_angles(alpha,beta,swap)
   local pta,ptb,pt
   swap = swap or false
   if swap  then 
     pta = rotation_(self.pa,-alpha,self.pb)
     ptb = rotation_(self.pb,beta,self.pa)
   else
   pta = rotation_(self.pa,alpha,self.pb)
   ptb = rotation_(self.pb,-beta,self.pa)
 end
   pt = intersection_ll_(self.pa,pta,self.pb,ptb)
   return triangle:new(self.pa,self.pb,pt)
end

function line:school(swap)
   local pta,ptb,pt
   swap = swap or false
   if swap  then 
     pta = rotation_(self.pa,-math.pi/6,self.pb)
     ptb = rotation_(self.pb,math.pi/3,self.pa)
   else
   pta = rotation_(self.pa,math.pi/6,self.pb)
   ptb = rotation_(self.pb,-math.pi/3,self.pa)
   end
   pt = intersection_ll_(self.pa,pta,self.pb,ptb)
   return triangle:new(self.pa,self.pb,pt)
end

function line:half(swap)
   local x,pt
   x  = midpoint_(self.pa,self.pb)
   swap = swap or false
   if swap  then 
    pt = rotation_(self.pb,math.pi/2,x)  
   else
   pt = rotation_(self.pb,-math.pi/2,x)
    end
   return triangle:new(self.pa,self.pb,pt)
end

function line:sss(a,b,swap)
   local pta,ptb,i,j
   swap = swap or false
    pta = self.pa + point ( a,  0 )
    ptb = self.pb + point ( -b , 0)
    i,j = intersection_cc_ (self.pa,pta,self.pb,ptb)
       if swap  then 
   return triangle : new (self.pa,self.pb,j)
 else 
   return triangle : new (self.pa,self.pb,i)
 end
end

function line:ssa(a, phi,swap)
    local x, y, i, j
     swap = swap or false
    x = rotation_(self.pb, -phi, self.pa)
    y = self.pa + polar_(a, self.slope)
    i, j = intersection_lc_(self.pb, x, self.pa, y)
    if swap  then  
     return  triangle:new(self.pa, self.pb, j)
    else
      return triangle:new(self.pa, self.pb, i)
    end
end

function line:sas(a, phi,swap)
    local x, pt
    swap = swap or false
    x = self.pa + polar_(a, self.slope)
     if swap  then 
        pt = rotation_(self.pa, -phi, x) 
     else
     pt = rotation_(self.pa, phi, x)
  end
   return triangle:new(self.pa, self.pb, pt)
end

function line:asa(alpha,beta,swap)
   local pta,ptb,pt
   swap = swap or false
    if swap  then 
      pta = rotation_ (self.pa,-alpha,self.pb)
      ptb = rotation_ (self.pb,beta,self.pa)
    else
   pta = rotation_ (self.pa,alpha,self.pb)
   ptb = rotation_ (self.pb,-beta,self.pa)
 end
 pt = intersection_ll_ (self.pa,pta,self.pb,ptb)
   return triangle : new (self.pa,self.pb,pt)
end
---- sacred triangles ----

function line:gold(swap)
    local pt
     swap = swap or false
     if swap  then  
        pt = rotation_ (self.pa,-math.pi/2,self.pb)
      return triangle : new (self.pa,self.pb, self.pa + (pt-self.pa) * tkzinvphi)
    else
         pt = rotation_ (self.pa,math.pi/2,self.pb)
      return triangle : new (self.pa,self.pb, self.pa + (pt-self.pa) * tkzinvphi)
    end
end

function line:sublime(swap)
   local pta,ptb,pt
   swap = swap or false
   if swap  then  
     pta = rotation_(self.pa,-2*math.pi/5,self.pb)
     ptb = rotation_(self.pb,2*math.pi/5,self.pa)
     pt = intersection_ll_(self.pa,pta,self.pb,ptb)
     return triangle:new(self.pa,self.pb,pt)
   else
  pta = rotation_(self.pa,2*math.pi/5,self.pb)
  ptb = rotation_(self.pb,-2*math.pi/5,self.pa)
  pt = intersection_ll_(self.pa,pta,self.pb,ptb)
  return triangle:new(self.pa,self.pb,pt)
end
end

line.euclid = line.sublime

function line:euclide (swap)
  swap = swap or false
  if swap  then  
     return triangle : new (self.pa,self.pb, rotation_(self.pa,-math.pi/5,self.pb))
  else
     return triangle : new (self.pa,self.pb, rotation_(self.pa,math.pi/5,self.pb))
     end
  end

  function line:divine(swap)
      local pta, ptb, pt
      swap = swap or false
      if swap  then  
        pta = rotation_(self.pa, -math.pi / 5, self.pb)
        ptb = rotation_(self.pb, math.pi / 5, self.pa)
        pt = intersection_ll_(self.pa, pta, self.pb, ptb)
        return triangle:new(self.pa, self.pb, pt)
      else
      pta = rotation_(self.pa, math.pi / 5, self.pb)
      ptb = rotation_(self.pb, -math.pi / 5, self.pa)
      pt = intersection_ll_(self.pa, pta, self.pb, ptb)
      return triangle:new(self.pa, self.pb, pt)
    end
  end

function line:cheops(swap)
    local m, n, pt
    m = midpoint_(self.pa, self.pb)
    swap = swap or false
    if swap  then 
      n = rotation_(m, math.pi / 2, self.pa)
     pt = m + (n - m) * tkzsqrtphi 
else
      n = rotation_(m, -math.pi / 2, self.pa)
     pt = m + (n - m) * tkzsqrtphi
  end
    return triangle:new(self.pa, self.pb, pt)
end


function line:egyptian(swap)
    local n, pt
    swap = swap or false
    if swap  then 
      n = rotation_(self.pb, math.pi / 2, self.pa)
      pt = self.pb + (n - self.pb) / point.mod(n - self.pb) * self.length * 0.75
   else  
    n = rotation_(self.pb, -math.pi / 2, self.pa)
    pt = self.pb + (n - self.pb) / point.mod(n - self.pb) * self.length * 0.75
  end
    return triangle:new(self.pa, self.pb, pt)
end

line.pythagoras = line.egyptian
line.isis = line.egyptian
line.golden = line.sublime
line.golden_gnomon = line.divine

------------------------------
-- Result -> square
------------------------------
function line:square (swap)
   swap = swap or false
   if swap  
   then
      return square : side (self.pa,self.pb,indirect)
   else
        return square : side (self.pa,self.pb)
   
   end
end


return line