# tem_determine_discreteVector Subroutine

## public subroutine tem_determine_discreteVector(vector, compareVector, angle)

Compare the incoming discrete vector against a set of prevailing vectors and return the found closest prevailing integer vector

A set of real numbered unit vectors (compareVector) must be given. The integer numbered vector is compared against those and the closest integer vector overwrites the original vector.

### Arguments

Type IntentOptional Attributes Name
integer, intent(inout) :: vector(:)

The given vector, will be set to the vector in compareVector that is best aligned to it.

real(kind=rk), intent(in) :: compareVector(:,:)

Set of unit vectors to select from against. size is (3, nVectors)

real(kind=rk), intent(out), optional :: angle

angle between vectors

## Source Code

  subroutine tem_determine_discreteVector( vector, compareVector, angle )
! -------------------------------------------------------------------- !
!> The given vector, will be set to the vector in compareVector
!! that is best aligned to it.
integer, intent(inout) :: vector(:)
!> Set of unit vectors to select from against.
!! size is (3, nVectors)
real( kind=rk ), intent(in) :: compareVector(:,:)
!> angle between vectors
real( kind=rk ), intent(out), optional :: angle
! -------------------------------------------------------------------- !
integer :: iDir, bestIndex
real(kind=rk) :: length, rv(3)
real(kind=rk) :: rVectIn(3), angleLoc, dotProduct, maxplen, min_nz_comp
! -------------------------------------------------------------------- !

angleLoc = 0.0_rk

! Only need to do anything for non-zero vectors.
if ( maxval(abs(vector)) > 0 ) then
rv = real(vector, kind=rk)

length = sqrt(rv(1)**2 + rv(2)**2 + rv(3)**2)

! Normalize vector
rVectIn = rv / length

! Keep track of the best fitting index.
bestIndex = 0

! Keep track of the maximal projected length.
! As we are looking at two unit vectors, the smallest
! projected length should be -1 (opposing directions).
! Thus, starting with -2 here as a save comparison that is
! smaller than all possible outcomes from compareVectors.
maxplen = -2.0_rk

! Go over all vectors in the compareVector
do iDir = 1, size(compareVector, 2)
! Go through all prevailing directions and check the minimum
dotProduct = rVectIn(1)*compareVector( 1, iDir ) &
&        + rVectIn(2)*compareVector( 2, iDir ) &
&        + rVectIn(3)*compareVector( 3, iDir )

dotProduct = min( max(dotProduct, -1.0_rk), 1.0_rk )
if ( dotProduct > maxplen ) then
maxplen = dotProduct
bestIndex = iDir
if ( maxplen .feq. 1._rk ) EXIT
end if
end do

! Align the vector to the closest matching prevailing direction from
! the set of compareVector:

! Use the smallest nonzero component to scale the unit vector to
! components close to integer numbers.
min_nz_comp = abs( minval( pack(compareVector(:, bestIndex),     &
&                             abs(compareVector(:, bestIndex)) &
&                             > epsilon(0.0_rk))               &
&                      )                                       )
! Use the nearest integers to represent the projected vector.
vector = nint(compareVector(:, bestIndex)/min_nz_comp)

! Angle between the original vector and the projected direction.
angleLoc = acos(maxplen)
end if ! not zero-vector

if (present(angle)) angle = angleLoc

end subroutine tem_determine_discreteVector