!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!  contact: delft3d.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands
!!
!!  All indications and logos of, and references to registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

      subroutine boombounce( xold, yold, xnew, ynew, nboom, xboom, &
                                    yboom, xcatch, ycatch, catch, xbounce, ybounce, bounce)
!
!     Does a boom catch a particle in it's path from xold, yold to xnew, ynew?
!     And if it does, does it bounce without crossing another part of the boom?

      use precision_part      ! single/double precision
      use timers

      implicit none

      real    ( rp ) :: xold, yold, xnew, ynew
      integer ( ip ) :: nboom
      integer ( ip ) :: npolbounce
      real    ( sp ) :: xboom(nboom), yboom(nboom)
      real    ( sp ) :: xcatch, ycatch
      logical        :: catch
      real    ( sp ) :: xbounce, ybounce
      logical        :: bounce

!     local variables
      real    ( rp ) :: xmin, ymin, xmax, ymax
      real    ( dp ) :: xoldd, yoldd, xnewd, ynewd
      real    ( dp ) :: xboomd(nboom), yboomd(nboom)
      real    ( dp ) :: a1, b1, c1, rl1, rlnear
      real    ( dp ) :: a2(nboom), b2(nboom), c2(nboom), rl2(nboom)
      real    ( dp ) :: xcross(nboom), ycross(nboom)
      real    ( sp ) :: xcrossboom(2),ycrossboom(2)
      real    ( dp ) :: rl, xn, yn, xb, yb
      integer ( 4  ) :: iboom
      integer ( 4  ) :: inear
      real    ( dp ) :: det
      logical        :: boomgap(nboom), notnear

      integer(4) ithndl              ! handle to time this subroutine
      data       ithndl / 0 /
      if ( timon ) call timstrt( "lines_cross", ithndl )

      catch = .false.
      xcatch = 999.999e0
      ycatch = 999.999e0
      bounce = .false.
      xbounce = 999.999e0
      ybounce = 999.999e0
      inear = 0
      rlnear = 1.1D0
      boomgap = .false.
      xcrossboom = 999.999e0
      ycrossboom = 999.999e0

      xoldd = dble(xold)
      yoldd = dble(yold)
      xnewd = dble(xnew)
      ynewd = dble(ynew)
      xboomd = dble(xboom)
      yboomd = dble(yboom)

!     at least one coordinate must have changed, otherwise there is no line between old and new location      
      if ( xold .ne. xnew .or. yold .ne. ynew ) then
!        fill matrix for particle track         
         a1 = ynewd - yoldd
         b1 = xoldd - xnewd
         c1 = xnewd * yoldd - xoldd * ynewd
         xmin = min(xold, xnew)
         xmax = max(xold, xnew)
         ymin = min(yold, ynew)
         ymax = max(yold, ynew)

!        loop over booms
         do iboom = 1, (nboom - 1)
!           skip 'gaps' when both coordinates of either start or end are 999.999e0 (missing value)
            boomgap(iboom) = (xboom(iboom) .eq. 999.999e0 .and. xboom(iboom) .eq. 999.999e0) .or. &
                          (xboom(iboom + 1) .eq. 999.999e0 .and. xboom(iboom + 1) .eq. 999.999e0)
!           when both old and new x- or y-coordinate are less or higher than both the boom x- or y-coordinate
            notnear = (xmax .lt. xboom(iboom) .and. xmax .lt. xboom(iboom+1)) .or. &
                      (xmin .gt. xboom(iboom) .and. xmin .gt. xboom(iboom+1)) .or. &
                      (ymax .lt. yboom(iboom) .and. ymax .lt. yboom(iboom+1)) .or. &
                      (ymin .gt. yboom(iboom) .and. ymin .gt. yboom(iboom+1))            
!           if both are false we might have a colission (sorry for the double denial)
            if (.not. boomgap(iboom) .and. .not. notnear) then
!              only when at least one coordinate is different we have a boom line               
               if (xboomd(iboom) .ne. xboomd(iboom+1) .or. yboomd(iboom) .ne. yboomd(iboom+1) ) then
!                 fill matrix for boom track         
                  a2(iboom) = yboomd(iboom+1) - yboomd(iboom)
                  b2(iboom) = xboomd(iboom) - xboomd(iboom+1)
                  c2(iboom) = xboomd(iboom+1) * yboomd(iboom) - xboomd(iboom) * yboomd(iboom+1)

!                 matrix determinant
                  det =  a1 * b2(iboom) - b1 * a2(iboom)

!                 when the determinant is zero, the lines are exactly parallel, and never cross
                  if ( det .ne. 0.0D+00 ) then
!                    calculate where the lines through the points would cross
                     xcross(iboom) = (b1 / det) * c2(iboom) - (b2(iboom) / det) * c1
                     ycross(iboom) = (a2(iboom) / det) * c1 - (a1 / det) * c2(iboom)

!                    relative distance of crosspoint between old and new location of particle
                     if (xnew .ne. xold) then 
                        rl1 = (xcross(iboom) - xoldd) / (xnewd - xoldd)
                     else
                        rl1 = (ycross(iboom) - yoldd) / (ynewd - yoldd)
                     end if
!                    relative location of crosspoint on the boom
                     if (xboom(iboom+1) .ne. xboom(iboom)) then
                        rl2(iboom) = (xcross(iboom) - xboom(iboom)) / (xboom(iboom+1) - xboom(iboom))
                     else
                        rl2(iboom) = (ycross(iboom) - yboom(iboom)) / (yboom(iboom+1) - yboom(iboom))
                     end if
!                    both relative locations must be between 0.0 and 1.0 to cross
                     if ( (rl1 .ge. 0.0D0 .and. rl1 .le. 1.0D0) .and. (rl2(iboom) .ge. 0.0D0 .and. rl2(iboom) .le. 1.0D0)) then
!                       determine if this is the lowest relative location of the particle path
                        if (rl1 .lt. rlnear) then
                           catch = .true.
                           rlnear = rl1
                           inear = iboom
                        end if
                     end if
                  end if
               end if
            end if
         end do
         if (catch) then
!           return information about collision location
            xcrossboom(1) = xboom(inear)
            xcrossboom(2) = xboom(inear+1)
            ycrossboom(1) = yboom(inear)
            ycrossboom(2) = yboom(inear+1)
            xcatch = sngl(xcross(inear))
            ycatch = sngl(ycross(inear))
            npolbounce = inear

!           bounce without crossing boom?
            rl  = ((xboomd(inear) - xnewd) * b2(inear) + (ynewd - yboomd(inear)) * a2(inear)) / (b2(inear)**2 + a2(inear)**2)
            xn = xboomd(inear) - rl * b2(inear)
            yn = yboomd(inear) + rl * a2(inear)
!           potential location after bounce at boom
            xb = xn - (xnewd - xn)
            yb = yn - (ynewd - yn)
            bounce = .true.

!           fill matrix for particle track         
            a1 = yb - yoldd
            b1 = xoldd - xb
            c1 = xb * yoldd - xoldd * yb
            xmin = min(xold, xb)
            xmax = max(xold, xb)
            ymin = min(yold, yb)
            ymax = max(yold, yb)

!           loop over booms, we saved some data from the previous loop (boomgap, a2, b2, c2)
            do iboom = 1, (nboom - 1)
!              when both old and new x- or y-coordinate are less or higher than both the boom x- or y-coordinate
               notnear = (xmax .lt. xboom(iboom) .and. xmax .lt. xboom(iboom+1)) .or. &
                         (xmin .gt. xboom(iboom) .and. xmin .gt. xboom(iboom+1)) .or. &
                         (ymax .lt. yboom(iboom) .and. ymax .lt. yboom(iboom+1)) .or. &
                         (ymin .gt. yboom(iboom) .and. ymin .gt. yboom(iboom+1))            
!              if both are false we might have a colission (sorry for the double denial)
               if (.not.(iboom.eq.npolbounce) .and. bounce .and. (.not.boomgap(iboom) .and. .not. notnear)) then
!                 matrix determinant
                  det =  a1 * b2(iboom) - b1 * a2(iboom)
!                 when the determinant is zero, the lines are exactly parallel, and never cross
                  if ( det .ne. 0.0D+00 ) then
!                    calculate where the lines through the points would cross
                     xcross(iboom) = (b1 / det) * c2(iboom) - (b2(iboom) / det) * c1
                     ycross(iboom) = (a2(iboom) / det) * c1 - (a1 / det) * c2(iboom)

!                    relative distance of crosspoint between old and new location of particle
                     if (xb .ne. xold) then 
                        rl1 = (xcross(iboom) - xoldd) / (xb - xoldd)
                     else
                        rl1 = (ycross(iboom) - yoldd) / (yb - yoldd)
                     end if
!                    relative location of crosspoint on the boom
                     if (xboom(iboom+1) .ne. xboom(iboom)) then
                        rl2(iboom) = (xcross(iboom) - xboom(iboom)) / (xboom(iboom+1) - xboom(iboom))
                     else
                        rl2(iboom) = (ycross(iboom) - yboom(iboom)) / (yboom(iboom+1) - yboom(iboom))
                     end if
!                    if both relative locations are between 0.0 and 1.0 we have crossed another part of the boom when bouncing,
!                    so boncing is not possible unless we would futher iterate
                     if ( (rl1 .ge. 0.0D0 .and. rl1 .le. 1.0D0) .and. (rl2(iboom) .ge. 0.0D0 .and. rl2(iboom) .le. 1.0D0)) then
                        bounce = .false.
                     end if
                  end if
               end if
            end do

!           return bouncing location
            xbounce = sngl(xb)
            ybounce = sngl(yb)
         end if
      end if
99999 if ( timon ) call timstop ( ithndl )
      return
      end subroutine boombounce      
	  