;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
; NAME: lineintersection
;
; PURPOSE: Calcule les coordonnees de l''intersection de 2 droites ou
; d''une serie de 2 droites
;
; CATEGORY:petit truc qui peut etre utile (sans boucles, ca va de soit!)
;
; CALLING SEQUENCE: point=lineintersection(abc1, abc2)
; 
; INPUTS: abc est un tableau de dimensions 3,
; nbre_de_couple_de_droites, dont chaque ligne contient les 3
; parametres a, b, c de l''equation de droite du type ax+by+c=0
;
; KEYWORD PARAMETERS:
;       /FLOAT: pour retourner l''output sous forme de tableau de
;       reel plutot que de vecteur decomplexes (par defaut)
;
; OUTPUTS:2 possibilites:
;      1) par defaut: c''est une vecteur de complexe dont chaque
;      element est les coordonnees du point d''intersection d''un
;      couple de droites.
;      2) si /FLOAT est active, c'est un tableau de reels de
;      dimensiones 2, nbre_de_couple_de_droites dont chaque
;      ligne est les coordonnees du point d''intersection d''un
;      couple de droites
;
; COMMON BLOCKS:
;
; SIDE EFFECTS:si les deux droites dont paralleles, on retourne les
; coordonnes (!values.f_nan,!values.f_nan)
;
; RESTRICTIONS:attention a la precision de la machine qui fait que
; les coordonnees calcules de verifient peut-etre pas exactement les
; equations du couple de droites.
;
; EXAMPLE:
;     IDL> abc1=linearequation(complex(1,2),[3,4])
;     IDL> abc2=linearequation(complex(1,2),[8,15])
;     IDL> print, lineintersection(abc1, abc2)
;     (      1.00000,      2.00000)
;     IDL> print, lineintersection(abc1, abc2,/float)
;           1.00000      2.00000
;
; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr)
;                      10 juin 2000
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
FUNCTION lineintersection, abc1, abc2, FLOAT = float
;
   a1 = float(reform(abc1[0, *]))
   b1 = float(reform(abc1[1, *]))
   c1 = float(reform(abc1[2, *]))
   a2 = float(reform(abc2[0, *]))
   b2 = float(reform(abc2[1, *]))
   c2 = float(reform(abc2[2, *]))
;
   determinant = a1*b2-a2*b1
   nan = where(determinant EQ 0)
   if nan[0] NE -1 THEN determinant = !values.f_nan
;
   x = (b1*c2-c1*b2)/determinant
   y = (c1*a2-a1*c2)/determinant
;
   if keyword_set(float) then begin
      npts = n_elements(x) 
      res = [reform(x, 1, npts, /over), reform(y, 1, npts, /over)]
   ENDIF ELSE res = complex(x, y)
   return, res
end