
module newton_interpolation


contains

    !
    !> \brief Calculates the coefficients of the newtonian interpolation
    !

    subroutine get_newton_interpolation_coefficients( n, xn, fn, c)

        implicit none

        integer, intent(in)  :: n !< Degree of the interpolatory polynomial
        real(8), dimension(0:n), intent(in)  :: xn !< x nodes
        real(8), dimension(0:n), intent(in)  :: fn !< f(x)
        real(8), dimension(0:n), intent(out) :: c  !< Coefficients of the interpolation

        ! Auxiliary variables

        integer :: i, j

        real(8), dimension( 0:n , 0:n ) :: M

        M(:,0) = fn

        do j = 1, n

            do i = j, n

                M(i,j) = ( M(i,j-1) - M(i-1,j-1) ) / ( xn(i) - xn(i-j) )

            end do

        end do

        do j = 0, n

            c(j) = M(j,j)

        end do

    end subroutine get_newton_interpolation_coefficients



    !
    !> \brief Calculates the value of the polynomial at a specified point x.
    !

    real(8) function get_newton_interpolation_f0(n, x, xn, c)
        implicit none
        integer, intent(in) :: n !< Degree of the interpolatory polynomial
        real(8), intent(in) :: x !< Value of x where the interpolation will be carried out
        real(8), dimension(0:n), intent(in) :: xn !< x nodes
        real(8), dimension(0:n), intent(in) :: c  !< Coefficients of the interpolation

        ! Auxiliary variables

        integer :: i, j

        real(8) :: aux1, aux2

        aux1 = 0.d0

        do i = 0, n

            aux2 = 1.d0

            do j = 0, i-1

                aux2 = aux2 * ( x - xn(j) )

            end do

            aux1 = aux1 + c(i) * aux2

        end do

        get_newton_interpolation_f0 = aux1

    end function get_newton_interpolation_f0

    !
    !> \brief Calculates the value of the polynomial first derivative at a specified point x.
    !

    real(8) function get_newton_interpolation_f1(n, x, xn, c)
        implicit none
        integer, intent(in) :: n !< Degree of the interpolatory polynomial
        real(8), intent(in) :: x  !< Value of x where the interpolation will be carried out
        real(8), dimension(0:n), intent(in) :: xn !< x nodes
        real(8), dimension(0:n), intent(in) :: c  !< Coefficients of the interpolation

        ! Auxiliary variables

        integer :: i, j, k

        real(8) :: aux1, aux2, aux3

        aux1 = 0.d0

        do i = 0, n

            aux2 = 0.d0

            do j = 0, i-1

                aux3 = 1.d0

                do k = 0, i-1

                    if ( k /= j ) then

                        aux3 = aux3 * ( x - xn(k) )

                    end if

                end do

                aux2 = aux2 + aux3

            end do

            aux1 = aux1 + c(i) * aux2

        end do

        get_newton_interpolation_f1 = aux1

    end function get_newton_interpolation_f1


    !> \brief Searches in the vector xn(1:n) the value of x (xc) where the interpolation
    !! will be carried out. If xc is found, then a interpolation of f(xc) with a FIRST degree
    !! polynomial is performed based on xn and fn. Depending on the results, the variable res
    !! returns: res = 1, meaning "no problems found", res = -1, meaning "out of range", and
    !! res = -2, meaning "error".
    subroutine get_ni_f0_vec_d1(n,xn,fn,xc,fc,res)
        implicit none
        integer, intent(in) :: n !< Number of nodes in the vector xn
        real(8), dimension(n), intent(in) :: xn !< Vector of the nodes
        real(8), dimension(n), intent(in) :: fn !< Vector of f(xn)
        real(8), intent(in)  :: xc  !< Value of x where the interpolation must be done
        real(8), intent(out) :: fc  !< Interpolation of f(xc)
        integer, intent(out) :: res !< Status of the result

        ! Auxiliary variables

        integer :: kc ! Index for xc

        real(8), dimension(0:1) :: xl ! Local values of xn
        real(8), dimension(0:1) :: fl ! Local values of fn
        real(8), dimension(0:1) ::  c ! Coefficients of the newtonian interpolation

        res = 1

        fc = 0.d0

        call searcher(n, xn, xc, kc)

        if ( kc == -1 ) then

            ! Out of range

            res = -1

            return

        else if ( kc == -2 ) then

            ! It was found an error

            res = -2

            return

        end if

        xl = xn( kc : kc+1 )

        fl = fn( kc : kc+1 )

        call get_newton_interpolation_coefficients( 1, xl, fl, c)

        fc = get_newton_interpolation_f0( 1, xc, xl, c)

    end subroutine


    !> \brief Searches in the vector xn(1:n) the value of x (xc) where the interpolation
    !! will be carried out. If xc is found, then a interpolation of f(xc) with a SECOND degree
    !! polynomial is performed based on xn and fn. Depending on the results, the variable res
    !! returns: res = 1, meaning "no problems found", res = -1, meaning "out of range", and
    !! res = -2, meaning "error".
    subroutine get_ni_f0_vec_d2(n,xn,fn,xc,fc,res)
        implicit none
        integer, intent(in) :: n !< Number of nodes in the vector xn
        real(8), dimension(n), intent(in) :: xn !< Vector of the nodes
        real(8), dimension(n), intent(in) :: fn !< Vector of f(xn)
        real(8), intent(in)  :: xc  !< Value of x where the interpolation must be done
        real(8), intent(out) :: fc  !< Interpolation of f(xc)
        integer, intent(out) :: res !< Status of the result

        ! Auxiliary variables

        integer :: kc ! Index for xc

        real(8), dimension(0:2) :: xl ! Local values of xn
        real(8), dimension(0:2) :: fl ! Local values of fn
        real(8), dimension(0:2) ::  c ! Coefficients of the newtonian interpolation

        res = 1

        fc = 0.d0

        call searcher(n, xn, xc, kc)

        if ( kc == -1 ) then

            ! Out of range

            res = -1

            return

        else if ( kc == -2 ) then

            ! It was found an error

            res = -2

            return

        end if

        if ( kc == 1 ) then

            xl = xn( kc : kc+2 )

            fl = fn( kc : kc+2 )

            call get_newton_interpolation_coefficients( 2, xl, fl, c)

            fc = get_newton_interpolation_f0(2, xc, xl, c)

        else if ( kc == n-1 ) then

            xl = xn( kc-1 : kc+1 )

            fl = fn( kc-1 : kc+1 )

            call get_newton_interpolation_coefficients( 2, xl, fl, c)

            fc = get_newton_interpolation_f0(2, xc, xl, c)

        else

            if ( abs( xc - xn(kc) ) < abs( xc - xn(kc+1) ) ) then

                xl = xn( kc-1 : kc+1 )

                fl = fn( kc-1 : kc+1 )

                call get_newton_interpolation_coefficients( 2, xl, fl, c)

                fc = get_newton_interpolation_f0(2, xc, xl, c)

            else

                xl = xn( kc : kc+2 )

                fl = fn( kc : kc+2 )

                call get_newton_interpolation_coefficients( 2, xl, fl, c)

                fc = get_newton_interpolation_f0(2, xc, xl, c)

            end if

        end if

    end subroutine


    !> \brief Searches for the xc in the vector x(1:n). If xc was found, then
    !! kc is the index in the vector x such that x(kc) <= xc <= x(kc+1), else
    !! kc = -1 if xc is out of [x(1),x(n)] and
    !! kc = -2 if an error was found.
    subroutine searcher(n, x, xc, kc)
        implicit none
        integer, intent(in) :: n                !< Number of nodes
        real(8), dimension(n), intent(in) :: x !< Vector of independent variable
        real(8), intent(in)  :: xc              !< Value of x to be found
        integer, intent(out) :: kc              !< Index in the vector such that x(kc) <= xc <= x(kc+1)

        ! Auxiliary variables

        integer:: i, ki, km, kf

        ki = 1
        kf = n

        if ( ( x(ki) - xc ) * ( x(kf) - xc ) > 0.d0 ) then

            ! xc do not belongs to the interval [x(ki), x(kf)]

            kc = -1

            return

        else

            do i = 1, n

                km = ( ki + kf ) / 2

                if ( ( x(ki) - xc ) * ( x(km) - xc ) > 0.d0 ) then

                    ki = km

                else

                    kf = km

                end if

                if ( kf - ki < 2 ) exit

            end do

            kc = ki

            if ( x(kc) > xc .or. x(kc+1) < xc ) then

                ! It was found an error!

                kc = -2

            end if

        end if

    end subroutine



end module newton_interpolation

