do_ray_intersect Function

   pure function do_ray_intersect(self, ray_origin, ray_direction) result(do_intersect)
   !< Return true if AABB is intersected by ray from origin and oriented as ray direction vector.
   class(aabb_object), intent(in) :: self          !< AABB box.
   type(vector_R8P),   intent(in) :: ray_origin    !< Ray origin.
   type(vector_R8P),   intent(in) :: ray_direction !< Ray direction.
   logical                        :: do_intersect  !< Test result.
   logical                        :: must_return   !< Flag to check when to return from procedure.
   real(R8P)                      :: tmin, tmax    !< Minimum maximum ray intersections with box slabs.

   do_intersect = .false.
   must_return = .false.
   tmin = 0._R8P
   tmax = MaxR8P
   call check_slab(aabb_min=self%bmin%x, aabb_max=self%bmax%x, &
                   o=ray_origin%x, d=ray_direction%x, must_return=must_return, tmin=tmin, tmax=tmax)
   if (must_return) return
   call check_slab(aabb_min=self%bmin%y, aabb_max=self%bmax%y, &
                   o=ray_origin%y, d=ray_direction%y, must_return=must_return, tmin=tmin, tmax=tmax)
   if (must_return) return
   call check_slab(aabb_min=self%bmin%z, aabb_max=self%bmax%z, &
                   o=ray_origin%z, d=ray_direction%z, must_return=must_return, tmin=tmin, tmax=tmax)
   if (must_return) return
   ! ray intersects all 3 slabs
   do_intersect = .true.
   contains
      pure subroutine check_slab(aabb_min, aabb_max, o, d, must_return, tmin, tmax)
      !< Perform ray intersection check in a direction-split fashion over slabs.
      real(R8P), intent(in)    :: aabb_min     !< Box minimum bound in the current direction.
      real(R8P), intent(in)    :: aabb_max     !< Box maximum bound in the current direction.
      real(R8P), intent(in)    :: o            !< Ray origin in the current direction.
      real(R8P), intent(in)    :: d            !< Ray slope in the current direction.
      logical,   intent(inout) :: must_return  !< Flag to check when to return from procedure.
      real(R8P), intent(inout) :: tmin, tmax   !< Minimum maximum ray intersections with box slabs.
      real(R8P)                :: ood, t1, t2  !< Intersection coefficients.
      real(R8P)                :: tmp          !< Temporary buffer.

      if ((d) < EPS) then
         ! ray is parallel to slab, no hit if origin not within slab
         if ((o < aabb_min).or.(o > aabb_max)) then
            must_return = .true.
            return
         endif
      else
         ! compute intersection t value of ray with near and far plane of slab
         ood = 1._R8P / d
         t1 = (aabb_min - o) * ood
         t2 = (aabb_max - o) * ood
         ! make t1 be intersection with near plane, t2 with far plane
         if (t1 > t2) then
            tmp = t1
            t1 = t2
            t2 = tmp
         endif
         ! compute the intersection of slab intersection intervals
         if (t1 > tmin) tmin = t1
         if (t2 > tmax) tmax = t2
         ! exit with no collision as soon as slab intersection becomes empty
         if (tmin > tmax) then
            must_return = .true.
            return
         endif
      endif
      endsubroutine check_slab
   endfunction do_ray_intersect