module utils

#ifdef INTEL
  use ifport, only : getcwd
#endif
  use const_and_precisions, only : wp_

  implicit none

contains

  pure subroutine locate(array, value, index)
    ! Given an `array`, and a `value` returns the `index`
    ! such that `array(index) < value < array(index+1)`
    !
    ! Notes:
    !   1. `array` must be monotonic, either increasing or decreasing.
    !   2. an index equal to 0 or `size(array)` indicate the value was not found
    ! Source: Numerical Recipes

    ! subroutine arguments
    real(wp_), intent(in)  :: array(:), value
    integer,   intent(out) :: index

    ! local variables
    integer :: jl, ju, jm
    logical :: incr

    jl = 0
    ju = size(array) + 1
    incr = array(size(array)) > array(1)

    do while ((ju - jl) > 1)
      jm = (ju + jl)/2
      if (incr .eqv. (value > array(jm))) then
        jl = jm
      else
        ju = jm
      endif
    end do
    index = jl
  end subroutine locate


  pure subroutine locate_unord(array, value, locs, n, nlocs)
    ! Given an `array` of size `n` and a `value`, finds at most
    ! `n` locations `locs` such that `value` is between
    ! `array(locs(i))` and `array(locs(i+i))`, in whichever order.

    ! subroutine arguments
    real(wp_), intent(in)    :: array(:)
    real(wp_), intent(in)    :: value
    integer,   intent(inout) :: locs(n)
    integer,   intent(in)    :: n
    integer,   intent(out)   :: nlocs

    ! local variables
    integer :: i
    logical :: larger_than_last

    nlocs = 0
    if (size(array) < 2) return

    larger_than_last = value > array(1)
    do i = 2, size(array)
      ! Note: the condition is equivalent to
      ! (array(i-1) < value < array(i))
      ! .or. (array(i) < value < array(i-1))
      if (array(i) < value .neqv. larger_than_last) then
        larger_than_last = value > array(i)
        nlocs = nlocs + 1
        if (nlocs <= n) locs(nlocs) = i - 1
      end if
    end do
  end subroutine locate_unord


  pure function linear_interp(x1, y1, x2, y2, x) result(y)
    ! Computes the linear interpolation between
    ! the points `(x1, y1)` and `(x2, y2)` at `x`.
    !
    ! Note: x1 != x2 is assumed

    ! function arguments
    real(wp_), intent(in) :: x1, y1, x2, y2, x
    real(wp_) :: y

    ! local variables
    real(wp_) :: a

    a = (x2 - x)/(x2 - x1)
    y = a*y1 + (1 - a)*y2
  end function linear_interp


  pure subroutine vmaxmin(x, xmin, xmax, imin, imax)
    ! Computes the maximum and minimum of the array `x`
    ! and optionally their indices

    ! subroutine arguments
    real(wp_), intent(in)            :: x(:)
    real(wp_), intent(out)           :: xmin, xmax
    integer,   intent(out), optional :: imin, imax

    ! local variables
    integer :: i

    if (size(x) < 1) then
      if (present(imin)) imin = 0
      if (present(imax)) imax = 0
      return
    end if

    if (present(imin)) imin = 1
    if (present(imax)) imax = 1
    xmin = x(1)
    xmax = x(1)

    do i = 2, size(x)
      if(x(i) < xmin) then
        xmin = x(i)
        if (present(imin)) imin = i
      else if(x(i) > xmax) then
        xmax = x(i)
        if (present(imax)) imax = i
      end if
    end do
  end subroutine vmaxmin


  pure subroutine sort_pair(p)
    ! Sorts the pair `p` in ascending order

    ! subroutine inputs
    real(wp_), intent(inout) :: p(2)

    ! local variables
    real(wp_) :: temp

    if (p(1) > p(2)) then
      temp = p(1)
      p(1) = p(2)
      p(2) = temp
    end if
  end subroutine sort_pair


  function dirname(filepath) result(directory)
    ! Get the parent `directory` of `filepath`

    ! function arguments
    character(*), intent(in)  :: filepath
    character(:), allocatable :: directory

    ! local variables
    character(255) :: cwd
    integer :: last_sep, err

    last_sep = scan(filepath, '/', back=.true.)
    directory = filepath(1:last_sep)

    ! append the cwd to relative paths
    if (isrelative(filepath)) then
      err = getcwd(cwd)
      directory = trim(cwd) // '/' // directory
    end if
  end function dirname


  function isrelative(filepath)
    ! Check if `filepath` is a relative or an absolute path

    ! function arguments
    character(*), intent(in)  :: filepath
    logical :: isrelative

    isrelative = (filepath(1:1) /= '/')
  end function isrelative


  pure function digits(n)
    ! Returns the number of digits of an integer

    ! function arguments
    integer, intent(in) :: n
    integer             :: digits

    ! How it works:
    !  1. `bit_size` returns the number of bits + 1 (sign)
    !  2. subtracting the number of leading zeros gives the significant bits
    !  3. multiplying by log₁₀(2) gives the approximate number of digits
    !  4. `ceiling` rounds to the next closest integer
    digits = ceiling((bit_size(abs(n)) - leadz(abs(n))) * log10(2.0_wp_))
  end function digits


  pure function diag(v) result(D)
    ! Returns a matrix D with v as main diagonal

    ! function arguments
    real(wp_), intent(in) :: v(:)
    real(wp_)             :: D(size(v), size(v))

    integer :: i

    D = 0
    do concurrent (i = 1:size(v))
      D(i,i) = v(i)
    end do
  end function


  pure function rotate(phi) result(R)
    ! Returns a 2D rotation matrix

    ! function arguments
    real(wp_), intent(in) :: phi
    real(wp_)             :: R(2,2)

    R = reshape([cos(phi), -sin(phi), &
                 sin(phi),  cos(phi)], shape=[2,2], order=[2,1])
  end function rotate


  pure function det(M)
    ! Computes the determinant of a real 3x3 matrix M

    ! function arguments
    real(wp_), intent(in) :: M(3,3)
    real(wp_)             :: det

    associate(a => M(1,1), b => M(1,2), c => M(1,3), &
              d => M(2,1), e => M(2,2), f => M(2,3), &
              g => M(3,1), h => M(3,2), i => M(3,3))
      det = (a*e*i + b*f*g + c*d*h) - (a*f*h + b*d*i + c*e*g)
    end associate
  end function det


  pure function eigenvals(A)
    ! Computes the eigenvalues of a real symmetric 3x3 matrix A
    !
    ! Source: https://doi.org/10.1145%2F355578.366316
    use const_and_precisions, only : pi, one

    ! function arguments
    real(wp_), intent(in) :: A(3,3)
    real(wp_)             :: eigenvals(3)

    ! local variables
    real(wp_) :: B(3,3)
    real(wp_) :: p, q, r, p1, p2
    real(wp_) :: phi
    integer   :: i

    p1 = A(1,2)**2 + A(1,3)**2 + A(2,3)**2

    if (p1 == 0) then
      ! A is diagonal
      eigenvals = [A(1,1), A(2,2), A(3,3)]
      return
    end if

    q  = sum([(A(i,i), i=1,3)])/3  ! trace(A)
    p2 = (A(1,1) - q)**2 + (A(2,2) - q)**2 + (A(3,3) - q)**2 + 2*p1
    p  = sqrt(p2 / 6)
    B  = (A - q * diag([one, one, one])) / p
    r = det(B) / 2

    ! Ensure r ∈ [-1,1]
    if (r <= -1) then
       phi = pi / 3
    else if (r >= 1) then
       phi = 0
    else
       phi = acos(r) / 3
    end if

    eigenvals = q - 2*p*[-cos(phi), cos(phi + pi/3), sin(phi + pi/6)]
  end function eigenvals


  pure function singvals(A)
    ! Computes the singular values of a real 3x3 matrix A

    ! function arguments
    real(wp_), intent(in) :: A(3,3)
    real(wp_)             :: singvals(3)

    ! By definition, the singular values of A are the
    ! eigenvalues value of A†A. These are necessarily
    ! real by the spectral theorem.
    singvals = sqrt(eigenvals(matmul(A, transpose(A))))
  end function singvals

end module utils