;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
;
; @file_comments 
; Calculate a linear equation of the type ax+by+c=0
; thanks to coordinates of 2 points.
; comment: we can have a table with pairs of points.
;
; @categories utilities
; 
; @param point1 {in}{required} This is the first point of(the) straight 
; line(s) whose we want to calculate equation(s)
;
; @param point2 {in}{required} This is the second point of(the) straight 
; line(s) whose we want to calculate equation(s)
;
;    There is 2 possibilities:
;      1) point is a complex or a table ofcomplex, where each element is the coordinates of the point.
;      2) point is a table of real of dimension 2,number_of_straight_line.
;         For each row of the table, we have coordinates of the point.
;
; @returns abc is a table of dimension 3, number_of_straight_line, 
;          where for each lign of the table we obtain the 3 parameters
;          a, b and c of the linear equation ax+by+c=0
;
; @examples IDL> abc=linearequation(complex(1,2),[3,4])
;           IDL> print, abc[0]*1+abc[1]*2+abc[2]
;           0.00000
;
; @history Sebastien Masson (smasson@lodyc.jussieu.fr)
;          10 juin 2000
;
; @version $Id$
;
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
FUNCTION linearequation, point1, point2
;
  compile_opt idl2, strictarrsubs
;

   if size(point1, /type) EQ 6 OR size(point1, /type) EQ 9 then begin
      x1 = float(point1)
      y1 = imaginary(point1)
   ENDIF ELSE BEGIN
      x1 = float(reform(point1[0, *]))
      y1 = float(reform(point1[1, *]))
   ENDELSE

   if size(point2, /type) EQ 6 OR size(point2, /type) EQ 9 then begin
      x2 = float(point2)
      y2 = imaginary(point2)
   ENDIF ELSE BEGIN
      x2 = float(reform(point2[0, *]))
      y2 = float(reform(point2[1, *]))
   ENDELSE

   vertical = where(x1 EQ x2)
   novertical = where(x1 NE x2)
   abc = fltarr(3, n_elements(x1)) 

   IF novertical[0] NE -1 then BEGIN
; y=mx+p
      nele = n_elements(novertical) 
      m = (y2[novertical]-y1[novertical])/(x2[novertical]-x1[novertical])
      p = (x2[novertical]*y1[novertical]-y2[novertical]*x1[novertical])/(x2[novertical]-x1[novertical])
      abc[*, novertical] = [reform(-m, 1, nele), replicate(1, 1, nele), reform(-p, 1, nele)]
   ENDIF
   IF vertical[0] NE -1 then BEGIN
; x=ny+p
      nele = n_elements(vertical) 
      n = (x2[vertical]-x1[vertical])/(y2[vertical]-y1[vertical])
      p = (y2[vertical]*x1[vertical]-x2[vertical]*y1[vertical])/(y2[vertical]-y1[vertical])
      abc[*, vertical] = [replicate(1, 1, nele), reform(-n, 1, nele), reform(-p, 1, nele)]
   ENDIF

   return, abc
end