source: trunk/SRC/Utilities/linearequation.pro @ 508

;+
;
; @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 of complex, 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 line 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 = real_part(point1)
    y1 = imaginary(point1)
  ENDIF ELSE BEGIN
    if size(point1, /type) EQ 5 then begin
      x1 = double(reform(point1[0, *]))
      y1 = double(reform(point1[1, *]))
    ENDIF ELSE BEGIN
      x1 = float(reform(point1[0, *]))
      y1 = float(reform(point1[1, *]))
    ENDELSE
  ENDELSE
  
  if size(point2, /type) EQ 6 OR size(point2, /type) EQ 9 then begin
    x2 = real_part(point2)
    y2 = imaginary(point2)
  ENDIF ELSE BEGIN
    if size(point2, /type) EQ 5 then begin
      x2 = double(reform(point2[0, *]))
      y2 = double(reform(point2[1, *]))
    ENDIF ELSE BEGIN
      x2 = float(reform(point2[0, *]))
      y2 = float(reform(point2[1, *]))
    ENDELSE
  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