obl_math.f90

module obl_math
    !< @brief general mathematics module
    ! --------------------------------------------------------------------------------

    ! TO DO:
    !   -- nearest neighbour interpolation/extrapolation (choice)
    !   -- periodic function support in interpolation

    ! modules used
    ! --------------------------------------------------------------------------------

    ! directives list
    implicit none
    private

    ! public interface
    ! --------------------------------------------------------------------------------
    public :: c_interp_linear
    public :: limit_min_array, limit_max_array
    public :: tma
    public :: is_finite, is_finite_array
    ! --------------------------------------------------------------------------------


    contains

    ! --------------------------------------------------------------------------------
    subroutine c_interp_linear(F_res, x_0, F, x, n) 
        !< @brief linear interpolation
        ! ----------------------------------------------------------------------------
        real, intent(out) :: F_res              !< interpolated value

        real, intent(in) :: x_0                 !< interpolated coordinate
        integer, intent(in) :: n                !< number of steps
        real, dimension(n), intent(in) :: F     !< input values
        real, dimension(n), intent(in)  :: x    !< input coorindates

        integer :: i
        ! ----------------------------------------------------------------------------
        
        if (n <= 0) then
            F_res = 0.0
            return
        endif

        if (x_0 <= x(1)) then
            F_res = F(1)
            return
        endif

        do i=1, n-1 
            if (x_0 <= x(i+1)) then
                F_res = F(i) + (F(i+1) - F(i)) / (x(i+1) - x(i)) * (x_0 - x(i))
                return
            endif
        enddo
        F_res = F(n)

    end subroutine c_interp_linear

    ! --------------------------------------------------------------------------------
    subroutine limit_min_array(F, F_min, n) 
        !< @brief limit function (minimum)
        ! ----------------------------------------------------------------------------
        integer, intent(in) :: n                    !< number of steps
        real, intent(inout), dimension(n) :: F      !< value on current z-step
        
        real, intent(in) :: F_min                   !< minimal value
        
        integer :: k
        ! ----------------------------------------------------------------------------

        do k=1, n
            if ( F(k) <  F_min  ) then
                F(k) = F_min  
            endif
        enddo
    end subroutine limit_min_array

    ! --------------------------------------------------------------------------------
    subroutine limit_max_array(F, F_max, n)
        !< @brief limit function (maximum)
        ! ----------------------------------------------------------------------------
        integer, intent(in) :: n                    !< number of steps
        real, intent(inout), dimension(n) :: F      !< value on current z-step

        real, intent(in) :: F_max                   !< maximum value

        integer :: k
        ! ----------------------------------------------------------------------------

        do k=1, n
            if (F(k) > F_max) then
                F(k) = F_max
            endif
        enddo    
    end subroutine limit_max_array

    ! --------------------------------------------------------------------------------
    subroutine tma(X,a,b,c,f,n)
        !< @brief solver for tridiagonal matrix (Thomas algorithm)
        ! ----------------------------------------------------------------------------
        integer, intent(in) :: n
        real, intent(out), dimension(n) :: X

        real, intent(in), dimension(n) :: a, b, c, f
        
        real, dimension(n) :: alpha, beta
        integer :: k
        ! ----------------------------------------------------------------------------

        alpha(1) = 0.0
        beta(1) = 0.0

        do k=1, n-1
            alpha(k+1) = -c(k) / (a(k)*alpha(k) + b(k))
            beta(k+1) = (f(k) -a (k)*beta(k)) / (a(k)*alpha(k) + b(k))
        end do

        X(n) = (f(n) - a(n)*beta(n)) / (b(n) + a(n)*alpha(n))
        
        do k = n, 2, -1
            X(k-1) = alpha(k)*X(k) + beta(k)         
        end do

    end subroutine tma

    ! --------------------------------------------------------------------------------
    elemental function is_finite(value)
        !> @brief check if value is finite
        ! ----------------------------------------------------------------------------
        use ieee_arithmetic
        implicit none
        logical :: is_finite

        real, intent(in) :: value
        ! ----------------------------------------------------------------------------

        is_finite = ieee_is_finite(value)
    end function is_finite
    
    function is_finite_array(F, n)
        !> @brief check if any value in array is finite
        ! ----------------------------------------------------------------------------
        use ieee_arithmetic
        implicit none
        logical :: is_finite_array
        
        integer, intent(in) :: n
        real, intent(in), dimension(n) :: F
        
        integer :: k
        ! ----------------------------------------------------------------------------

        is_finite_array = .true.
        do k=1, n
            is_finite_array = ieee_is_finite(F(k))
            if (.not.is_finite_array) exit
        enddo
    end function is_finite_array

end module