source: trunk/SRC/Utilities/lineintersection.pro @ 192

;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
;
; @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 reels 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 line are parallel, we return coordinates 
; (!values.f_nan,!values.f_nan)
;
; @restrictions 
; 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