;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
; NAME: linearequation
;
; PURPOSE:calcule une equation de droite du type ax+by+c=0 a partir
; des coordonnees de 2 points. Rq: on peut avoir un tableau de couple
; de points.
;
; CATEGORY:petit truc qui peut etre utile (sans boucles, ca va de soit!)
;
; CALLING SEQUENCE:abc=linearequation(point1, point2)
; 
; INPUTS: point1 et point2 dont deux point de(s) la droite(s) dont on
; veut calculer l(es) equations(s).
; deux possibilites sont possibles:
;      1) point est un complexe ou un tableau de complexes, ou chaque
;      element du complexe est les coordonnees du point
;      2) points est un tableau de reels de dimensions 2
;      ,nbre_de_droite. ou pour chaque ligne du tableau on a les
;      coordonnees du point.
;
; KEYWORD PARAMETERS:
;
; OUTPUTS:abc c''est un tableau de dimensions 3, nbre_de_droite, ou
; pour chaque ligne du tableau on obtient les 3 parametres a, b, c de
; l'equation de la droite ax+by+c=0
;
; COMMON BLOCKS:
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;
; EXAMPLE:
;    IDL> abc=linearequation(complex(1,2),[3,4])
;    IDL> print, abc[0]*1+abc[1]*2+abc[2]
;          0.00000
;
; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr)
;                       10 juin 2000
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
FUNCTION linearequation, point1, point2

   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