; docformat = 'rst'
;
; NAME: 
;    interLineCircle
; PURPOSE:
;    Evaluate the coordinates of the intersection(s) between a circle
;    and a segment or a line
;
;+
; :Description:
;    Evaluate the coordinates of the intersection(s) between a circle
;    and a segment or a line
;
; :Categories:
;    Maths, Intersection
;    
; :Params:
;    A: in, required, type=fltarr(2)
;        Cartesian coordinates of the first point of segment [A-B]
;    B: in, required, type=fltarr(2)
;        Cartesian coordinates of the second point of segment [A-B]
;    PARAM: in, required, type=float
;        The parameters of the circle::
;          P[0] = X center
;          P[1] = Y center
;          P[2] = Radius
;
; :Returns: The cartesian coordinates of the intersection
;
; :Keywords:
;    FLAG: out, optional, type=integer
;        Provides a summary of the intersection result::
;          0: No intersection
;          1: One Intersection
;          2: Two Intersections
;          3: Null Segment or Degenerate Circle
;          4: Bug
;    LINE: in, optional, type=boolean, default=0
;        Treat segment A-B as a line
;    DIST: out, optional, type=float
;        Distance of the intersection from A
;
; :Author:
;    B.Carry (IMCCE)
;
; :History:
;     Change History::
;        Original Version written in June 2013, B. Carry (IMCCE)
;        2013 Sep. - B.Carry (IMCCE) - Header cleaned
;-
function interLineCircle, A, B, param, flag=flag, line=line, dist=dist

 ;--------------------------------------------------------------------------------
 ;--I-- Initialization And Input Verification  -----------------------------------
 ;--I.1-- Check Presence of Arguments
  if N_params() LT 3 then begin
    print,'Syntax -  Y =  interLineCircle( A, B, P )'
    return, -1
 endif

 ;--I.2-- Flag Initialization
  flag=4

 ;--I.3-- Null Segment or Degenerate Circle
  if (A(0) eq B(0) and A(1) eq B(1)) or param(2) eq 0 then begin
    flag = 3
    return, [0,0]
  endif


 ;--------------------------------------------------------------------------------
 ;--II-- Shift and Rotate the System  --------------------------------------------
 ;--II.1-- Translate the System so A is at the Origin
  inB = B - float(A)
  C = param(0:1) - float(A)

 ;--II.2-- Length of Segment A-B
  distAB = sqrt( inB(0)*inB(0) + inB(1)*inB(1))

 ;--II.3-- Rotate the System so B is on the positive X axis.
  theCos = inB(0)/distAB
  theSin = inB(1)/distAB

  newX =  C(0)*theCos + C(1)*theSin
  C(1) = -C(0)*theSin + C(1)*theCos
  C(0) = newX


 ;--------------------------------------------------------------------------------
 ;--III-- Search for Intersections  -----------------------------------------------
  delta = ( param(2)*param(2)-C(1)*C(1) )  ;-Discriminant

 ;------------------------------------------------------------
 ;--III.1-- No intersection
  if delta lt 0 then begin
    flag=0
    return, [0,0]

  endif else begin

   ;------------------------------------------------------------
   ;--III.2-- Single intersection
    if delta eq 0 then begin

     ;--III.2.1-- Position in the A-B coordinate system
      crossPos = C(0)
      dist=crossPos

      ;--III.2.2-- Position in the Original Coordinate System
      xOut = A(0) + crossPos*theCos
      yOut = A(1) + crossPos*theSin

     ;--III.2.3-- Intersection in the Segment-Circle or Line-Circle Cases
      if ( crossPos ge 0 and crossPos le distAB ) OR keyword_set(line) then begin
        flag=1
        return, [xOut,yOut]

     ;--III.2.4-- Intersection Outside the Segment AND Segment-Circle Case
      endif else begin
        flag=0
        return, [0,0]
      endelse


   ;------------------------------------------------------------
   ;--III.3-- Two intersections
    endif else begin

     ;--III.3.1-- Positions in the A-B Coordinate System
      crossPos1 = C(0) + sqrt( delta )
      crossPos2 = C(0) - sqrt( delta )
      dist=[crossPos1, crossPos2]

     ;--III.3.2-- Position in the Original Coordinate System
      xOut1 = A(0) + crossPos1*theCos
      yOut1 = A(1) + crossPos1*theSin
      xOut2 = A(0) + crossPos2*theCos
      yOut2 = A(1) + crossPos2*theSin

     ;--III.3.3-- Double Intersection in the Segment-Circle or Line-Circle Cases
      if keyword_set(line) or $
        ( ( crossPos1 ge 0 and crossPos1 le distAB ) AND $
          ( crossPos2 ge 0 and crossPos2 le distAB ) ) then begin
       flag=2
       return, [ [xOut1,yOut1], [xOut2,yOut2] ]

      endif else begin

       ;--III.3.4-- Both Intersections Outside the Segment AND Segment-Circle Case
        if ( ( crossPos1 lt 0 and crossPos1 gt distAB ) AND $
             ( crossPos2 lt 0 and crossPos2 gt distAB ) ) then begin
          flag=0
          return, [0,0]
    
       ;--III.3.5-- One Intersection inside the Segment AND Segment-Circle Case
        endif else begin
    
         ;--III.3.5.1-- Intersection Corresponds to the First Point
          if ( crossPos1 ge 0 and crossPos1 le distAB ) then begin
            flag=1
            return, [xOut1,yOut1]

         ;--III.3.5.2-- Intersection Corresponds to the Second Point
          endif else begin
            if ( crossPos2 ge 0 and crossPos2 le distAB ) then begin
              flag=1
              return, [xOut2,yOut2]

           ;--III.3.5.3-- Deadend: No Intersection in the Segment
            endif else begin
              flag=0
              return, [0,0]

            endelse

          endelse

        endelse

      endelse

    endelse

  endelse




end