;+
;
; @file_comments
; Calculate coordinates of the intersection between 2 straight lines
; or of a succession of 2 straight lines.
;
; @categories
; Utilities
;
; @param ABC1 {in}{required}{type=3d array}
; is the first array of dimension 3, number_of_pairs_of_straight_lines,
; whose each line contain the 3 parameters a, b and c of the first linear
; equation of the type ax+by+c=0
;
; @param ABC2 {in}{required}{type=3d array}
; is second array of dimension 3, number_of_pairs_of_straight_lines,
; whose each line contain the 3 parameters a, b and c of the second linear
; equation of the type ax+by+c=0
;
; @keyword FLOAT
; To return the output as a array of real numbers instead of vectors of
; complex (by default)
;
; @returns
; 2 possibilities:
;  1) by default: it is a vector of complex whose each element is the 
;     coordinates of the intersection point of a pair of straight lines.
;  2) if FLOAT is activated, it is a array of reals of dimension 2,
;     number_of_pairs_of_straight_lines whose each row is the coordinates
;     of the intersection point of a pair of straight line.
;
; @restrictions
; If the 2 straight lines are parallel, we return coordinates
; (!values.f_nan,!values.f_nan)
;
; Beware of the precision of the machine which make
; that calculated coordinates may not exactly verify
; equations of the pair of straight lines.
;
; @examples
;
;   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
;
; @history
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
;          10 juin 2000
;
; @version
; $Id$
;
;-
FUNCTION lineintersection, abc1, abc2, FLOAT=float
;
  compile_opt idl2, strictarrsubs
;
   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