tem_intersec_line_line Subroutine

public subroutine tem_intersec_line_line(edge, line, intersects, intersection)

This subroutine calculates the intersection between a line and a line. It gives back the coordinates of the intersection, the multiple of the direction vector of the intersection and the distance of the intersection to the center point of the line.

Arguments

Type IntentOptional Attributes Name
type(tem_line), intent(in) :: edge
type(tem_line), intent(in) :: line
logical, intent(out) :: intersects
type(tem_intersec), intent(out) :: intersection

Called by

proc~~tem_intersec_line_line~~CalledByGraph proc~tem_intersec_line_line tem_intersec_line_line proc~exit_element exit_element proc~exit_element->proc~tem_intersec_line_line

Source Code

  subroutine tem_intersec_line_line( edge, line, intersects, intersection )
    ! ---------------------------------------------------------------------------
    type(tem_line), intent(in) :: edge
    type(tem_line), intent(in) :: line
    type(tem_intersec), intent(out) :: intersection
    logical, intent(out) :: intersects
    ! ---------------------------------------------------------------------------
    real(kind=rk), dimension(3) :: diff_vector, normal
    real(kind=rk), dimension(3) :: enormal
    real(kind=rk) :: alignment
    real(kind=rk) :: dist_line, dist_edge
    ! ---------------------------------------------------------------------------

    ! check whether the two lines intersect
    ! They have to be in a common plane for an
    ! intersection, compute the normal of this
    ! plane
    normal(1) = edge%direction(2)*line%direction(3)                            &
      &       - edge%direction(3)*line%direction(2)
    normal(2) = edge%direction(3)*line%direction(1)                            &
      &       - edge%direction(1)*line%direction(3)
    normal(3) = edge%direction(1)*line%direction(2)                            &
      &       - edge%direction(2)*line%direction(1)

    alignment = normal(1)*normal(1)                                            &
      &       + normal(2)*normal(2)                                            &
      &       + normal(3)*normal(3)

    if ((alignment > 16*tiny(alignment))) then
      ! The lines are not colinear, they might intersect, compute
      ! the distance of parallel planes through the lines, if
      ! this is 0, they actually intersect.
      dist_line = dot_product(normal, line%coordStart)
      dist_edge = dot_product(normal, edge%coordStart)
      intersects = (abs(dist_line - dist_edge) < epsilon(dist_line))

      if (intersects) then
        ! They intersect, get the point of intersection
        diff_vector = edge%coordStart - line%coordStart
        enormal(1) = edge%direction(2)*normal(3)                               &
          &        - edge%direction(3)*normal(2)
        enormal(2) = edge%direction(3)*normal(1)                               &
          &        - edge%direction(1)*normal(3)
        enormal(3) = edge%direction(1)*normal(2)                               &
          &        - edge%direction(2)*normal(1)
        intersection%lambda = dot_product(diff_vector, enormal)                &
          &                 / dot_product(line%direction, enormal)
        intersection%coord = line%coordStart + intersection%lambda             &
          &                * line%direction
        intersection%distance = edge%coordStart - intersection%coord
      end if
    else
      ! Lines are colinear
      intersects = .false.
    end if

  end subroutine tem_intersec_line_line