plane_line_intersection Subroutine

   subroutine plane_line_intersection(norm, V0, P0, P1, I, check, dist)
     use modglobal, only : vec0, eps1
     implicit none
     ! determines the intersection of a plane and a line segment
     ! norm: plane normal
     ! V0: point on the plane
     ! P0: start of line segment
     ! P1: end of line segment
     ! I: intersection point
     ! dist: distance from P0 to intersection point
     ! check: 0 if no intersection
     !        1 if unique intersection
     !        2 if line segment is in the plane
     !        3 if intersection is outside line segment
     real, intent(in),  dimension(3) :: norm, V0, P0, P1
     real, intent(out), dimension(3) :: I
     integer, intent(out) :: check
     real, intent(out) :: dist
     real, dimension(3) :: u, w
     real :: D, N, sI

     I = vec0
     w = P0 - V0
     u = P1 - P0
     D = dot_product(norm, u)
     N =-dot_product(norm, w)

     if (abs(D) < eps1) then ! line orthogonal to plane normal -> segment parallel to plane
       if (abs(N) < eps1) then ! start point is on the plane -> segment lies in the plane
         check = 2
         return
       else
         check = 0
         return
       end if
     end if

     sI = N / D
     I = P0 + sI * u
     dist = norm2(I - P0)

     if ((sI < 0.) .or. (sI > 1.)) then
       check = 3
     else
       check = 1
     end if

   end subroutine plane_line_intersection