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.
Type | Intent | Optional | 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 |
The angle betweend va and vb
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