module functions

    use parameters

    implicit none

    public
    contains


    ! Functions below this comment are taken from the functions module of atomsk
    !********************************************************
    !  STRUPCASE
    !  This function reads a string of any length
    !  and capitalizes all letters.
    !********************************************************
    FUNCTION StrUpCase (input_string) RESULT (UC_string)
    !
    IMPLICIT NONE
    CHARACTER(*),PARAMETER:: lower_case = 'abcdefghijklmnopqrstuvwxyz'
    CHARACTER(*),PARAMETER:: upper_case = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
    CHARACTER(*),INTENT(IN):: input_string
    CHARACTER(LEN(Input_String)):: UC_string !Upper-Case string that is produced
    INTEGER:: i, n
    !
    IF(LEN(input_string)==0) RETURN
    UC_string = input_string
    ! Loop over string elements
    DO i=1,LEN(UC_string)
      !Find location of letter in lower case constant string
      n = INDEX( lower_case, UC_string(i:i) )
      !If current substring is a lower case letter, make it upper case
      IF(n>0) THEN
        UC_string(i:i) = upper_case(n:n)
      ENDIF
    END DO
    !
    END FUNCTION StrUpCase

    !********************************************************
    !  STRDNCASE
    !  This function reads a string of any length
    !  and transforms all letters to lower case.
    !********************************************************
    FUNCTION StrDnCase (input_string) RESULT (lc_string)
    !
    IMPLICIT NONE
    CHARACTER(*),PARAMETER:: lower_case = 'abcdefghijklmnopqrstuvwxyz'
    CHARACTER(*),PARAMETER:: upper_case = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
    CHARACTER(*),INTENT(IN):: input_string
    CHARACTER(LEN(Input_String)):: lc_string !Lower-Case string that is produced
    INTEGER:: i, n
    !
    IF(LEN(input_string)==0) RETURN
    lc_string = input_string
    ! Loop over string elements
    DO i=1,LEN(lc_string)
      !Find location of letter in upper case constant string
      n = INDEX( upper_case, lc_string(i:i) )
      !If current substring is an upper case letter, make it lower case
      IF(n>0) THEN
        lc_string(i:i) = lower_case(n:n)
      ENDIF
    END DO
!
END FUNCTION StrDnCase

    pure function matrix_normal(a, n)

        integer :: i
        integer, intent(in) :: n
        real(kind = dp), dimension(n) :: v
        real(kind = dp), dimension(n, n),intent(in) :: a
        real(kind = dp), dimension(n,n) :: matrix_normal

        matrix_normal(:, :) = a(:, :)

        do i = 1, n

            v(:) = a(i,:)

            matrix_normal(i, :) = v(:) / norm2(v)

        end do

        return
    end function matrix_normal

    pure function cross_product(a, b)
        !Function which finds the cross product of two vectors

        real(kind = dp), dimension(3), intent(in) :: a, b
        real(kind = dp), dimension(3) :: cross_product

        cross_product(1) = a(2) * b(3) - a(3) * b(2)
        cross_product(2) = a(3) * b(1) - a(1) * b(3)
        cross_product(3) = a(1) * b(2) - a(2) * b(1)

        return
    end function cross_product

    pure function identity_mat(n)
        !Returns the nxn identity matrix

        integer :: i
        integer, intent(in) :: n
        real(kind = dp), dimension(n, n) :: identity_mat

        identity_mat(:, :) = 0.0_dp
        do i = 1, n
          identity_mat(i, i) = 1.0_dp
        end do

     return
     end function identity_mat

     pure function triple_product(a, b, c)
        !triple product between three 3*1 vectors

        real(kind = dp) :: triple_product
        real(kind = dp), dimension(3), intent(in) :: a, b, c
        triple_product = dot_product(a, cross_product(b, c))

     return
     end function triple_product

    function in_bd_lat(v, box_faces, box_norms)
        !This function returns whether the point is within the transformed box boundaries. The transformed 
        !box being the transformed simulation cell in the lattice basis
        
        !Input/output variables
        real(kind=dp), dimension(3), intent(in) :: v !integer lattice position
        real(kind=dp), dimension(3,6), intent(in) :: box_faces !Centroid of all the box faces
        real(kind=dp), dimension(3,6), intent(in) :: box_norms !Box face normals
        logical :: in_bd_lat

        !Other variables
        integer :: i
        real(kind=dp) :: pt_to_face(3)

        in_bd_lat = .true.

        !Check if point is in box bounds, this works by comparing the dot product of the face normal and the 
        !vector from the point to the face. If the dot product is greater than 0 then the point is behind the face
        !if it is equal to zero then the point is on the face, if is less than 0 the the point is in front of the face.
        do i = 1, 6
            pt_to_face(:) = box_faces(:, i) - v 
            if(dot_product(pt_to_face, box_norms(:,i)) <= 0) then 
                in_bd_lat = .false.
                exit
            end if
        end do

        return
    end function in_bd_lat

    function in_block_bd(v, box_bd)
        !This function returns whether a point is within a block in 3d

        !Input/output
        real(kind=dp), dimension(3), intent(in) :: v
        real(kind=dp), dimension(6), intent(in) :: box_bd
        logical :: in_block_bd

        !Other variables
        integer :: i

        in_block_bd = .true.

        do i =1 ,3
            !Check upper bound
            if(v(i) > (box_bd(2*i)+10.0_dp**(-10)) ) then 
                in_block_bd =.false.
                exit
            !Check lower bound
            else if (v(i) < (box_bd(2*i-1)-10.0_dp**(-10))) then 
                in_block_bd = .false.
                exit
            end if
        end do
    end function in_block_bd

    function lcm(a,b)
        !This function returns the smallest least common multiple of two numbers 

        real(kind=dp), intent(in) :: a, b
        real(kind=dp) :: lcm

        integer :: aint, bint, gcd, remainder, placeholder

        !Cast the vector positions to ints. There will be some error associated with this calculation
        aint = a*10**2
        bint = b*10**2

        !Calculate greated common divisor
        gcd = aint
        placeholder = bint
        do while(placeholder /= 0)
            remainder = modulo(gcd, placeholder)
            gcd = placeholder
            placeholder=remainder
        end do
        lcm = real((aint*bint),dp)/(real(gcd,dp))* 10.0_dp**(-2.0_dp)
    end function lcm

    function is_neighbor(rl, rk, r_in, r_out)
        !This function checks to see if two atoms are within a shell with an inner radius r_in and outer radius 
        !r_out
        real(kind=dp), intent(in) :: r_in, r_out
        real(kind=dp), dimension(3), intent(in) :: rl, rk
        logical :: is_neighbor

        !Internal variable
        real(kind=dp) :: rlk

        rlk = norm2(rk - rl)
        is_neighbor=.true.
        if((rlk>r_out).or.(rlk < r_in)) is_neighbor = .false.

        return
    end function is_neighbor
    
    function is_equal(A, B)
        !This function checks if too numbers are equal within a tolerance
        real(kind=dp), intent(in) :: A, B
        logical :: is_equal

        if((A>(B - 10.0_dp**(-10))).and.(A < (B+10.0_dp**(-10)))) then 
            is_equal = .true.
        else
            is_equal = .false.
        end if
    return
    end function is_equal

    pure function unitvec(n,vec)
        integer, intent(in) :: n
        real(kind=dp), intent(in) :: vec(n)
        real(kind=dp) :: unitvec(n)

        unitvec = vec/norm2(vec)
        return
    end function unitvec

end module functions