polynomials_m.f90 Source File

module polynomials_m
    !! Contains functionality related to polynomials, such as evaluation and
    !! derivates. Polynomials are represented by their coeffcients given as
    !! arrays.

    use precision_m, only: FP

    implicit none
    private

    public :: polyval, polyder

    interface polyval
        procedure :: polyval_1d
        procedure :: polyval_2d
    end interface

contains

    real(FP) pure function polyval_1d(polycoeff, x)
        !! Evaluate a 1D polynomial p given by its coefficients array
        !! at position x
        real(FP), contiguous, dimension(:), intent(in) :: polycoeff
        !! Polynomial coefficients
        real(FP), intent(in) :: x
        !! Point where to evaluate the polynomial
        integer :: n

        polyval_1d = 0.0_FP
        do n = size(polycoeff), 2, -1
            polyval_1d = (polyval_1d + polycoeff(n)) * x
        end do
        polyval_1d = polyval_1d + polycoeff(1)
    end function

    real(FP) pure function polyval_2d(polycoeff, x, y)
        !! Evaluate a 2D polynomial p given by its coefficients c_{n,m} array
        !! at position (x,y)
        !!  p(x,y) = \sum_{n,m} c_{n,m} * x^n * y^m
        !!         = \sum_{m} y^m * ( \sum_{n} c_{n,m} * x^n )
        real(FP), contiguous, dimension(:,:), intent(in) :: polycoeff
        !! Polynomial coefficients
        real(FP), intent(in) :: x
        !! Point x-coordinate where to evaluate the polynomial
        real(FP), intent(in) :: y
        !! Point y-coordinate where to evaluate the polynomial
        real(FP), dimension(1:ubound(polycoeff, 2)) :: polyval_x
        integer :: m

        do m = 1, ubound(polycoeff, 2)
            ! Collapse to single polynomial
            polyval_x(m) = polyval_1d(polycoeff(:,m), x)
        end do
        ! Calculate remaining polynomial
        polyval_2d = polyval_1d(polyval_x, y)
    end function

    pure function polyder(poly) result(dpoly)
        !! Returns the derivative of a 1D polynomial given by its coefficients
        !! array
        real(FP), contiguous, dimension(:), intent(in) :: poly
        !! Polynomial coefficients
        real(FP), dimension(:), allocatable :: dpoly
        !! Polynomial coefficients of the derivative
        integer :: n, degree_orig, degree_der

        degree_orig = size(poly)
        if(degree_orig==1) then
            allocate(dpoly(1), source=0.0_FP)
            return
        endif

        degree_der = degree_orig - 1
        allocate(dpoly(1:degree_der), source=0.0_FP)
        do n = 1, degree_der
            dpoly(n) = n * poly(n + 1)
        end do
    end function

end module