Module CartQuadXsectArea

  implicit none

  real(8), parameter :: ZERO = 0._8
  real(8), parameter :: ONEHALF = .5_8
  real(8), parameter :: ONE = 1._8

  type FPOINT
    real(8) :: x,y
  end type FPOINT

  type LINESEGMENT
    type(FPOINT) :: a,b
  end type LINESEGMENT

  CONTAINS

  function pnt_in_triangle(Tri_x, Tri_y) result(InTriangle)

  real(8), intent(in)  :: Tri_x(0:3), Tri_y(0:3)
  logical              :: InTriangle

  real(8) :: a, b, c

  a = (Tri_x(0) - Tri_x(3))*(Tri_y(1) - Tri_y(3)) - (Tri_x(1) - Tri_x(3))*(Tri_y(0) - Tri_y(3)) 
  b = (Tri_x(1) - Tri_x(3))*(Tri_y(2) - Tri_y(3)) - (Tri_x(2) - Tri_x(3))*(Tri_y(1) - Tri_y(3)) 
  c = (Tri_x(2) - Tri_x(3))*(Tri_y(0) - Tri_y(3)) - (Tri_x(0) - Tri_x(3))*(Tri_y(2) - Tri_y(3)) 

  InTriangle = a*b >= ZERO .and. b*c >= ZERO

  end function pnt_in_triangle

  function pnt_in_quad(Pnt, Quad_x, Quad_y) result(InQuad)

  real(8), intent(in)  :: Pnt(0:1)
  real(8), intent(in)  :: Quad_x(0:3)
  real(8), intent(in)  :: Quad_y(0:3)
  logical              :: InQuad

  real(8) :: Tri_x(0:3), Tri_y(0:3)
  logical :: InTriangle

  Tri_x(0:2) = Quad_x(0:2)
  Tri_y(0:2) = Quad_y(0:2)
  Tri_x(3) = pnt(0)
  Tri_y(3) = pnt(1)
  InTriangle = pnt_in_triangle(Tri_x,Tri_y)
  if( .not. InTriangle ) then
    Tri_x(1:2) = Quad_x(2:3)
    Tri_y(1:2) = Quad_y(2:3)
    InTriangle = pnt_in_triangle(Tri_x,Tri_y)
  endif

  InQuad = InTriangle

  end function pnt_in_quad

  function poly_area(x, y) result(Area)

    real(8), intent(in)  :: x(0:)
    real(8), intent(in)  :: y(0:)
    real(8)              :: Area

    integer :: nv, i,im1,ip1

    nv = size( x )
    Area = ZERO
    do i = 0,nv-1
      ip1 = mod((i+1),nv)
      im1 = i - 1
      if( im1 < 0 ) im1 = nv - 1
      Area = Area + (x(ip1) - x(im1))*y(i)
    enddo

    Area = -ONEHALF*Area

  end function poly_area

  function quadintersectarea(Quad1_x, Quad1_y, Quad2_x, Quad2_y) result(XsectArea)
! -----------------------------------------------------------
!  find area of intersection of two quadrilaterals in x,y plane
! -----------------------------------------------------------
 
  real(8), intent(in)  :: Quad1_x(0:3), Quad1_y(0:3)
  real(8), intent(in)  :: Quad2_x(0:3), Quad2_y(0:3)
  real(8)              :: XsectArea

   type(FPOINT)      :: matchPnt
   type(FPOINT)      :: intersect(0:3)
   type(FPOINT)      :: points(0:15)
   type(LINESEGMENT), target  :: Quad(0:3,0:1)
   type(LINESEGMENT), pointer :: TestLine
   type(LINESEGMENT), pointer :: MatchLine
   type(LINESEGMENT) :: lines(0:15)

   integer :: lineIndex
   integer :: i, n, p, r
   integer :: cnt, inum
   integer :: quadIndex
   integer :: index
   integer :: vtxcnt
   integer :: pointNumber
   real(8) :: A
   real(8) :: Area
   real(8) :: sAn, sAnp1
   real(8) :: vertex(0:1)
   real(8) :: aArray(0:3), sortedA(0:3)
   real(8) :: poly_x(0:15), poly_y(0:15)
   logical :: InQuad
   logical :: alreadyDone
   logical :: mask(0:15)

!------------------------------------------------------------
!  first a few diagnostics
!------------------------------------------------------------
   area = poly_area(Quad1_x, Quad1_y)
   area = poly_area(Quad2_x, Quad2_y)

! -----------------------------------------------------------
!  check for complete inclusion
! -----------------------------------------------------------
   p = 0
   r = 0
!------------------------------------------------------------
!  is quad1 "inside" quad2?
!------------------------------------------------------------
   do i = 0,3
     vertex(:) = (/ Quad1_x(i),Quad1_y(i) /)
     InQuad = pnt_in_quad(vertex,Quad2_x,Quad2_y)
     if( .not. InQuad ) then
       exit
     else
       p = p + 1
     endif
   enddo

!------------------------------------------------------------
!  is quad2 "inside" quad1?
!------------------------------------------------------------
   if( p /= 4 ) then
     p = 0
     do i = 0,3
       vertex(:) = (/ Quad2_x(i),Quad2_y(i) /)
       InQuad = pnt_in_quad(vertex,Quad1_x,Quad1_y)
       if( .not. InQuad ) then
         exit
       else
         p = p + 1
       endif
     enddo
     if( p == 4 ) r = 1
   endif
!  if( p == 4 ) then
   if( p == 40 ) then
     if( r == 0 ) then
       XsectArea = poly_area(Quad1_x, Quad1_y)
     else
       XsectArea = poly_area(Quad2_x, Quad2_y)
     endif
     return
   endif

! -----------------------------------------------------------
!  Transfer external quad to internal quad representation
! -----------------------------------------------------------
   do i = 0,3
     p = mod(i+1,4)
!     Quad one
     quad(i,0)%a%x = Quad1_x(i)            ! Point A
     quad(i,0)%a%y = Quad1_y(i)
     quad(i,0)%b%x = Quad1_x(p)            ! Point B
     quad(i,0)%b%y = Quad1_y(p)
!     Quad two
     quad(i,1)%a%x = Quad2_x(i)            ! Point A
     quad(i,1)%a%y = Quad2_y(i)
     quad(i,1)%b%x = Quad2_x(p)            ! Point B
     quad(i,1)%b%y = Quad2_y(p)
   enddo

   lineIndex = 0
! -----------------------------------------------------------
!  Get the intersection line segments for the quads...
! -----------------------------------------------------------
Match_quad_loop: &
   do quadIndex = 0,1
Test_quad_loop: &
     do i = 0,3
       MatchLine => quad(i,quadIndex)
! -----------------------------------------------------------
!  Check intersctions for every line in other quad...
! -----------------------------------------------------------
       index = 0
       do n = 0,3
         TestLine => Quad(n,1-quadIndex)
         if( .not. Parallel(TestLine,MatchLine) ) then
	   A = LineTest(TestLine, MatchLine)
	   if(A >= ZERO .and. A <= ONE) then
             A = LineTest(MatchLine, TestLine)
             intersect(index)%x = GetX(MatchLine, A)
             intersect(index)%y = GetY(MatchLine, A)
             aArray(index) = A
             index = index + 1
	   endif
         endif
       enddo
       if (index == 0) then
	 cycle Test_quad_loop
       endif
! -----------------------------------------------------------
!  Remove xsection point doubles...
! -----------------------------------------------------------
       mask(:) = .false.
       do p = 0,index-2
         if( .not. mask(p) ) then
           do r = p+1,index-1
             mask(r) = PointEquals(intersect(r),intersect(p))
           enddo
         endif
       enddo

       cnt = 1

       do p = 1,index-1
         if( .not. mask(p) ) then
           intersect(cnt) = intersect(p)
           aArray(cnt)    = aArray(p)
           cnt = cnt + 1
         endif
       enddo

       index = cnt
       if( index == 0 ) then
	 cycle Test_quad_loop
       endif
        
! -----------------------------------------------------------
!  Sort the lists smallest a value to largest a value...
! -----------------------------------------------------------
       do n = 0,index-1
         do inum = 0,n-1
	   if (aArray(n) < sortedA(inum)) then
! -----------------------------------------------------------
!  shift list to fit segment in correct place...
! -----------------------------------------------------------
	     do p = n,inum+1,-1
               sortedA(p) = sortedA(p - 1)
	     enddo
	     sortedA(inum) = aArray(n)
	     exit
	   endif
         enddo
	 if (inum == n) then
	   sortedA(inum) = aArray(n)
	 endif
       enddo
! -----------------------------------------------------------
!  create line segments and clip line segments
!  to fit current quad and add to a list...
! -----------------------------------------------------------
LineSeg_loop: &
       do n = 0,index-1,2
         sAn   = sortedA(n)
         sAnp1 = sortedA(n+1)
	 if (sAn < ZERO) then
	   if (sAnp1 > ONE) then
	     lines(lineIndex) = MatchLine
	   elseif (sAnp1 >= ZERO) then
	     lines(lineIndex)%a%x = GetX(MatchLine, sAnp1)
             lines(lineIndex)%a%y = GetY(MatchLine, sAnp1)
             lines(lineIndex)%b = MatchLine%b
	   endif
	 elseif (sAn > ONE) then
	   if (sAnp1 <= ONE .and. sAnp1 >= ZERO) then
	     lines(lineIndex)%b%x = GetX(MatchLine, sAnp1)
             lines(lineIndex)%b%y = GetY(MatchLine, sAnp1)
             lines(lineIndex)%a = MatchLine%a
	   endif
	 else
	   if (sAnp1 > ONE) then
             lines(lineIndex)%b%x = GetX(MatchLine, sAn)
             lines(lineIndex)%b%y = GetY(MatchLine, sAn)
             lines(lineIndex)%a = MatchLine%a
	   elseif (sAnp1 < ZERO) then
             lines(lineIndex)%a%x = GetX(MatchLine, sAn)
             lines(lineIndex)%a%y = GetY(MatchLine, sAn)
             lines(lineIndex)%b = MatchLine%b
	   else
             lines(lineIndex)%a%x = GetX(MatchLine, sAn)
             lines(lineIndex)%a%y = GetY(MatchLine, sAn)
             lines(lineIndex)%b%x = GetX(MatchLine, sAnp1)
             lines(lineIndex)%b%y = GetY(MatchLine, sAnp1)
	   endif
	 endif
         lineIndex = lineIndex + 1
       enddo LineSeg_loop
     enddo Test_quad_loop
   enddo Match_quad_loop

! ------------------------------------------------------
!  if there are no intersections area is zero...
! ------------------------------------------------------
   if (lineIndex == 0) then
     XsectArea = ZERO
     return 
   endif
! ------------------------------------------------------
!  Remove lines with zero length
! ------------------------------------------------------
   do i = 0,lineIndex-1
     mask(i) = PointEquals(lines(i)%a, lines(i)%b)
   enddo

   r = 0
   do p = 0,lineIndex-1
     if( .not. mask(p) ) then
       lines(r) = lines(p)
       r = r + 1
     endif
   enddo
   lineIndex = r
   if (lineIndex == 0) then
     XsectArea = ZERO
     return 
   endif
! ------------------------------------------------------
!  Remove line doubles
! ------------------------------------------------------
   mask(0:lineIndex-1) = .false.
   do i = 0,lineIndex-1
     if( .not. mask(i) ) then
       do n = i+1,lineIndex-1
         if( .not. mask(n) ) then
           mask(n) = LineEquals(lines(i), lines(n))
         endif
       enddo
     endif
   enddo
   r = 0
   do p = 0,lineIndex-1
     if( .not. mask(p) ) then
       lines(r) = lines(p)
       r = r + 1
     endif
   enddo
   lineIndex = r
! ------------------------------------------------------
!  if there are no intersections area is zero...
! ------------------------------------------------------
   if (lineIndex == 0) then
     XsectArea = ZERO
     return 
   endif
! ------------------------------------------------------
!  find adjacent lines to points...
! ------------------------------------------------------
   vtxcnt = 0
   pointNumber = 0
   do i = 0,lineIndex-1
     alreadyDone = .false.
     do n = 0,pointNumber-1
       if (PointEquals(points(n), lines(i)%a)) then
         alreadyDone = .true.
       endif
     enddo
     if( .not. alreadyDone ) then
       points(pointNumber) = lines(i)%a
       pointNumber = pointNumber + 1
     endif
     alreadyDone = .false.
     do n = 0,pointNumber-1
       if (PointEquals(points(n), lines(i)%b)) then
         alreadyDone = .true.
       endif
     enddo
     if( .not. alreadyDone ) then
       points(pointNumber) = lines(i)%b
       pointNumber = pointNumber + 1
     endif
   enddo
! ------------------------------------------------------
!  remove unconnected line segments...
! ------------------------------------------------------
   mask(0:lineIndex-1) = .true.
   do r = 0,lineIndex-1
     do p = 0,lineIndex-1
       if( p /= r ) then
         if( PointEquals(lines(r)%a,lines(p)%a) .or. &
             PointEquals(lines(r)%a,lines(p)%b) ) then
           exit
         endif
       endif
     enddo
     if( p /= lineIndex ) then
       do p = 0,lineIndex-1
         if( p /= r ) then
           if( PointEquals(lines(r)%b,lines(p)%a) .or. &
               PointEquals(lines(r)%b,lines(p)%b) ) then
             exit
           endif
         endif
       enddo
       mask(r) = p /= lineIndex
     else
       mask(r) = .false.
     endif
   enddo

   vtxcnt = count( mask(0:pointNumber-1) )

   if( vtxcnt < 3 ) then
     XsectArea = ZERO
     return 
   endif
        
! ------------------------------------------------------
!  setup vertices of intersection polygon
! ------------------------------------------------------
   cnt = 2
   do p = 0,pointNumber-1
     if( mask(p) ) then
       poly_x(0) = lines(p)%a%x
       poly_y(0) = lines(p)%a%y
       poly_x(1) = lines(p)%b%x
       poly_y(1) = lines(p)%b%y
       matchPnt  = lines(p)%b
       mask(p) = .false.
       exit
     endif
   enddo

   do while( cnt < vtxcnt )
     do p = 0,pointNumber-1
       if( mask(p) ) then
         if( PointEquals(lines(p)%a,matchPnt) ) then
           poly_x(cnt) = lines(p)%b%x
           poly_y(cnt) = lines(p)%b%y
           matchPnt  = lines(p)%b
           mask(p) = .false.
           cnt = cnt + 1
           exit
         elseif( PointEquals(lines(p)%b,matchPnt) ) then
           poly_x(cnt) = lines(p)%a%x
           poly_y(cnt) = lines(p)%a%y
           matchPnt  = lines(p)%a
           mask(p) = .false.
           cnt = cnt + 1
           exit
         endif
       endif
     enddo
   enddo

   XsectArea = abs(poly_area(poly_x(0:vtxcnt-1), poly_y(0:vtxcnt-1)))

   end function quadintersectarea

   function LineTest( one,two ) result( A )

   type(LINESEGMENT), intent(in) :: one, two
   real(8)                       :: A

   real(8) :: numer, denom

   numer = (one%b%y - two%b%y)*(two%a%x - two%b%x) - (one%b%x - two%b%x)*(two%a%y - two%b%y)
   denom = (one%a%y - one%b%y)*(two%a%x - two%b%x) - (one%a%x - one%b%x)*(two%a%y - two%b%y)
   A     = -numer/denom

   end function LineTest

   function GetX( line,A ) result(X)
     
   type(LINESEGMENT), intent(in) :: line
   real(8), intent(in)           :: A
   real(8)                       :: X

   X = line%a%x * A + line%b%x * (ONE - A)

   end function GetX

   function GetY( line,A ) result(Y)
     
   type(LINESEGMENT), intent(in) :: line
   real(8), intent(in)           :: A
   real(8)                       :: Y

   Y = line%a%y * A + line%b%y * (ONE - A)

   end function GetY

   function Parallel( one,two ) result(ParallelLines)

   type(LINESEGMENT), intent(in) :: one, two
   logical                       :: ParallelLines

   ParallelLines = ((one%a%y - one%b%y)*(two%a%x - two%b%x) - (one%a%x - one%b%x)*(two%a%y - two%b%y)) == ZERO

   end function Parallel


   function DoubleEquals( one,two ) result(Equal)

   real(8), intent(in)  :: one, two
   logical              :: Equal

   Equal = abs(one - two) < 0.0000000001_8

   endfunction DoubleEquals

   function PointEquals( one,two ) result(Equal)

   type(FPOINT), intent(in) :: one, two
   logical                  :: Equal

   Equal = DoubleEquals(one%x,two%x) .and. DoubleEquals(one%y,two%y)

   end function PointEquals


   function LineEquals( one,two ) result(Equal)

   type(LINESEGMENT), intent(in) :: one, two
   logical                       :: Equal

   Equal  = PointEquals(one%a,two%a) .and. PointEquals(one%b,two%b)
   if( .not. Equal ) then
     Equal = PointEquals(one%a,two%b) .and. PointEquals(one%b,two%a)
   endif

   end function LineEquals

  end module CartQuadXsectArea