; docformat = 'rst'
;
; NAME: 
;    INTERLINE
; PURPOSE:
;    Compute the position of the intersection between two lines or
;    segments of lines.
;
;+
; :Description:
;    Compute the position of the intersection between two lines or
;    segments of lines.
;
;    From the cartesian coordinates of two pairs of points, the code
;    compute where the two lines intersect and return the cartesian
;    coordinates of the intersection. The intersection can be searched
;    only inside one, or both, pair of point (line-segment or
;    segment-segment intersection) 
;
; :Categories:
;    Maths, Intersection
;    
; :Params:
;    A: in, required, type=fltarr(2)
;        Cartesian coordinates of the first point of line #1 (A-B)
;    B: in, required, type=fltarr(2)
;        Cartesian coordinates of the second point of line #1 (A-B)
;    C: in, required, type=fltarr(2)
;        Cartesian coordinates of the first point of line #2 (C-D)
;    D: in, required, type=fltarr(2)
;        Cartesian coordinates of the second point of line #2 (C-D)
;
; :Returns: The cartesian coordinates of the intersection
;
; :Keywords:
;    flag: out, optional, type=integer
;        Provides a summary of the intersection result::
;          -1: Invalid input
;          0: No intersection
;          1: Intersection
;          2: Share a point
;          3: Zero length segment
;          4: Bug
;    distance: out, optional, type=float
;        Distance of the intersection from A (positive toward B)
;    line1: in, optional, type=boolean, default=0
;        Treat segment A-B as a line
;    line2: in, optional, type=boolean, default=0
;        Treat segment C-D as a line
;
; :Examples:
;    Search for the intersection between line (A-B) and (C-D), and its
;    distance to A::
;      IDL> A=[0,0]
;      IDL> B=[2,2]
;      IDL> C=[1,0]
;      IDL> D=[1,2]
;      IDL> cross = interLine(A,B,C,D,distance=fromA)
;      IDL> print, cross, fromA
;
; :Author:
;    B.Carry (OCA)
;
; :History:
;     Change History::
;        Original Version written in June 2012, B. Carry (ESA)
;        2013 Nov. - B.Carry (IMCCE) - Changed I/O keywords
;        2016 Aug. - B.Carry (IMCCE) - idl2 added - input check added
;-
function interLine, A, B, C, D, flag=flag, distance=distance, line1=line1, line2=line2
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---    I -- Initialization And Input Verification              -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
  COMPILE_OPT hidden, idl2
  flag=4
  distance=0

 ;--I.1-- Valid Input Variable  ----------------------------------------------------------------
  if n_params() lt 4 then begin
    message, /ioError, 'Syntax: xy = interLine( A, B, C, D [, flag=, distance=, /line1, /line2)'
    flag = -1
    return, [0,0]
  endif

 ;--I.2-- Check against zero-length segments  --------------------------------------------------
  if (A[0] eq B[0] and A[1] eq B[1]) or $
     (C[0] eq D[0] and C[1] eq D[1]) then begin
    flag = 3
    return, [0,0]
  endif

 ;--I.3-- Check if Segments share a point  -----------------------------------------------------
 if (A[0] eq C[0] and A[1] eq C[1]) or $
    (A[0] eq D[0] and A[1] eq D[1]) then begin
    flag = 2
    return, A
  endif
 if (B[0] eq C[0] and B[1] eq C[1]) or $
    (B[0] eq D[0] and B[1] eq D[1]) then begin
    flag = 2
    return, B
  endif

;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---   II -- Rotate and Translate the System                    -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
 ;--II.1-- Translate system so A is at origin  -------------------------------------------------
  B -= A
  C -= A
  D -= A

 ;--II.2-- Length of [AB]  ---------------------------------------------------------------------
  distAB = sqrt( B[0]*B[0] + B[1]*B[1]);

 ;--II.3-- Rotate system so B is on positive X axis  -------------------------------------------
  theCos = B[0]/distAB;
  theSin = B[1]/distAB;

  newX = C[0]*theCos+C[1]*theSin
  C[1] = C[1]*theCos-C[0]*theSin
  C[0] = newX

  newX = D[0]*theCos+D[1]*theSin
  D[1] = D[1]*theCos-D[0]*theSin
  D[0] = newX


 ;--II.4-- Find intersection along line (AB)  --------------------------------------------------
  crossPos= D[0]+(C[0]-D[0])*D[1]/(D[1]-C[1])
  distance = crossPos

 ;--II.5-- Rotate back the system  -------------------------------------------------------------
  xOut = A[0] + crossPos*theCos;
  yOut = A[1] + crossPos*theSin;



;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---  III -- Case of [CD] is a segment                          -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
  if not keyword_set(line2) then begin

   ;--III.1-- No Intersection  -----------------------------------------------------------------
    if (C[1] lt 0 and D[1] lt 0) or $
       (C[1] gt 0 and D[1] gt 0) then begin
      flag=0
      return, [0,0]
    endif

   ;--III.2-- Intersection Outside [AB]  -------------------------------------------------------
    if ( crossPos lt 0 or crossPos gt distAB ) then begin

     ;--III.2.1-- (AB) is a Line  ----------------------------------------------------------------
      if keyword_set(line1) then begin
        flag = 1
        return, [xOut, yOut]

     ;--III.2.2-- [AB] is a Segment  -------------------------------------------------------------
      endif else begin
        flag=0
        return, [0,0]
      endelse

   ;--III.3-- Intersection Inside [AB]  --------------------------------------------------------
    endif else begin
      flag = 1
      return, [xOut, yOut]
    endelse


;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---   IV -- Case of (CD) is a Line                             -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
  endif else begin

   ;--IV.1-- Intersection Outside [AB]  --------------------------------------------------------
    if ( crossPos lt 0 or crossPos gt distAB ) then begin

     ;--IV.1.1-- (AB) is a Line  ---------------------------------------------------------------
      if keyword_set(line1) then begin
        flag = 1
        return, [xOut, yOut]

     ;--IV.1.2-- [AB] is a Segment  ------------------------------------------------------------
      endif else begin
        flag=0
        return, [0,0]
      endelse


   ;--IV.2-- Intersection Inside [AB]  ---------------------------------------------------------
    endif else begin
       flag=1
       return, [xOut, yOut]
    endelse


  endelse
end