angle_between Function

private elemental function angle_between(va_x, va_y, vb_x, vb_y) result(angle)

Compute the angle between to vectors (they should not both be the 0 vector).

This function uses the crossproduct and arctan to find the angle between two 2D vectors. It takes the four components of the vectors as scalar arguments to allow vectorized evaluation of multiple vector pairs at once. If one of the vectors is zero, the returned angle is also zero. Also multiplicities of Pi are ignored and a 0 angle will be returned for them.

Arguments

Type IntentOptional Attributes Name
real(kind=rk), intent(in) :: va_x

The first vector va
real(kind=rk), intent(in) :: va_y

The first vector va
real(kind=rk), intent(in) :: vb_x

The second vector vb
real(kind=rk), intent(in) :: vb_y

The second vector vb

Return Value real(kind=rk)

The angle betweend va and vb


Source Code

  elemental function angle_between(va_x, va_y, vb_x, vb_y) result(angle)
    ! ----------------------------------------------------------------------
    !> The first vector va
    real(kind=rk), intent(in) :: va_x, va_y

    !> The second vector vb
    real(kind=rk), intent(in) :: vb_x, vb_y

    !> The angle betweend va and vb
    real(kind=rk) :: angle
    ! ----------------------------------------------------------------------

    angle = 0.0_rk

    ! Only proceed if both vectors are non-zero.
    ! If one of the vectors is 0, we assume a zero angle.
    if (      (abs(va_x)+abs(va_y) > 2*tiny(va_x)) &
      & .and. (abs(vb_x)+abs(vb_y) > 2*tiny(vb_x)) ) then
      ! The angle is arctan( (va X vb)/(va*vb) ):
      ! Using atan2 here for numerical stability and correct signs.
      angle = atan2( (va_x*vb_y - va_y*vb_x), &
        &            (va_x*vb_x + va_y*vb_y)  )

      ! If the angle is close enough to Pi, consider this as a colinear
      ! 0 angle (avoid multiplicities of Pi in the summed angle).
      if (abs(angle) > tolm*Pi) angle = 0.0_rk
    end if

  end function angle_between