02_equispaced_interpolation_linear_1D.f90

Go to the documentation of this file.

!----------------------------------------------------------------------------------------
!
! This file is part of EFTCAMB.
!
! Copyright (C) 2013-2016 by the EFTCAMB authors
!
! The EFTCAMB code is free software;
! You can use it, redistribute it, and/or modify it under the terms
! of the GNU General Public License as published by the Free Software Foundation;
! either version 3 of the License, or (at your option) any later version.
! The full text of the license can be found in the file eftcamb/LICENSE at
! the top level of the EFTCAMB distribution.
!
!----------------------------------------------------------------------------------------



!----------------------------------------------------------------------------------------


module equispaced_linear_interpolation_1d

    use precision
    use amlutils
    use eftcamb_mixed_algorithms

    implicit none

    private

    public equispaced_linear_interpolate_function_1d

    !----------------------------------------------------------------------------------------
    type :: equispaced_linear_interpolate_function_1d

        ! options:
        integer                             :: num_points

        ! parameters:
        real(dl)                            :: x_initial
        real(dl)                            :: x_final
        real(dl)                            :: null_value
        logical                             :: has_null_value
        real(dl)                            :: grid_width

        ! arrays with the values:
        real(dl), allocatable, dimension(:) :: x
        real(dl), allocatable, dimension(:) :: y
        real(dl), allocatable, dimension(:) :: yp
        real(dl), allocatable, dimension(:) :: ypp
        real(dl), allocatable, dimension(:) :: yppp
        real(dl), allocatable, dimension(:) :: yint

    contains

        procedure :: initialize             => equispacedlinearintepolatefunction1dinit               
        procedure :: precompute             => equispacedlinearintepolatefunction1dprecompute         
        procedure :: value                  => equispacedlinearintepolatefunction1dvalue              
        procedure :: first_derivative       => equispacedlinearintepolatefunction1dfirstderivative    
        procedure :: second_derivative      => equispacedlinearintepolatefunction1dsecondderivative   
        procedure :: third_derivative       => equispacedlinearintepolatefunction1dthirdderivative    
        procedure :: integral               => equispacedlinearintepolatefunction1dintegral           
        procedure :: initialize_derivatives => equispacedlinearintepolatefunction1dinitderivatives    

    end type equispaced_linear_interpolate_function_1d

    ! ---------------------------------------------------------------------------------------------

contains

    ! ---------------------------------------------------------------------------------------------
    subroutine equispacedlinearintepolatefunction1dinit( self, num_points, x_initial, x_final, null_value )

        implicit none

        class(equispaced_linear_interpolate_function_1d) :: self

        integer , intent(in)                             :: num_points
        real(dl), intent(in)                             :: x_initial
        real(dl), intent(in)                             :: x_final
        real(dl), intent(in), optional                   :: null_value

        integer  :: i

        ! initialize the null value:
        if ( present(null_value) ) then
            self%has_null_value = .true.
            self%null_value     = null_value
        else
            self%has_null_value = .false.
            self%null_value     = 0._dl
        end if

        ! allocate the x vector:
        if ( allocated(self%x) ) deallocate( self%x )

        self%num_points = num_points

        allocate( self%x( self%num_points ) )

        ! store initial and final times:
        self%x_initial  = x_initial
        self%x_final    = x_final

        ! compute the width of the interpolation:
        self%grid_width = ( self%x_final -self%x_initial )/REAL( self%num_points -1 )

        ! fill in the equispaced x array:
        do i=1, self%num_points
            self%x(i)   = self%x_initial + REAL(i-1)*self%grid_width
        end do

        ! allocate the other vectors:
        if ( allocated(self%y)    ) deallocate( self%y    )
        if ( allocated(self%yp)   ) deallocate( self%yp   )
        if ( allocated(self%ypp)  ) deallocate( self%ypp  )
        if ( allocated(self%yppp) ) deallocate( self%yppp )
        if ( allocated(self%yint) ) deallocate( self%yint )

        allocate( self%y( self%num_points )    )
        allocate( self%yp( self%num_points )   )
        allocate( self%ypp( self%num_points )  )
        allocate( self%yppp( self%num_points ) )
        allocate( self%yint( self%num_points ) )

    end subroutine equispacedlinearintepolatefunction1dinit

    ! ---------------------------------------------------------------------------------------------
    subroutine equispacedlinearintepolatefunction1dprecompute( self, x, ind, mu )

        implicit none

        class(equispaced_linear_interpolate_function_1d)  :: self
        real(dl), intent(in)                              :: x
        integer , intent(out)                             :: ind
        real(dl), intent(out)                             :: mu

        real(dl) :: x1, x2

        ! compute the interpolation index:
        ind = int( ( x-self%x_initial)/self%grid_width ) +1
        ! store the x values:
        x1  = self%x(ind)
        x2  = self%x(ind+1)
        ! compute the linear interpolation coefficient:
        mu  = (x-x1)/(x2-x1)

    end subroutine equispacedlinearintepolatefunction1dprecompute

    ! ---------------------------------------------------------------------------------------------
    function equispacedlinearintepolatefunction1dvalue( self, x, index, coeff )

        implicit none

        class(equispaced_linear_interpolate_function_1d)  :: self
        real(dl), intent(in)                              :: x
        integer , intent(in), optional                    :: index
        real(dl), intent(in), optional                    :: coeff
        real(dl) :: equispacedlinearintepolatefunction1dvalue
        integer  :: ind
        real(dl) :: x1, x2, y1, y2, mu

        ! initialize to null value:
        equispacedlinearintepolatefunction1dvalue = self%null_value
        if ( self%has_null_value ) then
            ! if outside the interpolation range return the null value:
            if ( x <= self%x_initial .or. x >= self%x_final ) return
        else
            ! if below the interpolation range return the first value:
            if ( x <= self%x_initial ) then
                equispacedlinearintepolatefunction1dvalue = self%y(1)
                return
            end if
            ! if above the interpolation range return the first value:
            if ( x >= self%x_final   ) then
                equispacedlinearintepolatefunction1dvalue = self%y(self%num_points)
                return
            end if
        end if

        ! return the index of the point:
        if ( present(index) ) then
            ind = index
        else
            ind = int( ( x-self%x_initial)/self%grid_width ) +1
        end if

        ! get the interpolation coefficient:
        if ( present(coeff) ) then
            mu = coeff
        else
            ! store the x values:
            x1  = self%x(ind)
            x2  = self%x(ind+1)
            ! compute the linear interpolation coefficient:
            mu  = (x-x1)/(x2-x1)
        end if

        ! store the y values:
        y1  = self%y(ind)
        y2  = self%y(ind+1)

        ! compute the linear interpolation:
        equispacedlinearintepolatefunction1dvalue = y1*( 1._dl -mu ) +y2*mu

    end function equispacedlinearintepolatefunction1dvalue

    ! ---------------------------------------------------------------------------------------------
    function equispacedlinearintepolatefunction1dfirstderivative( self, x, index, coeff )

        implicit none

        class(equispaced_linear_interpolate_function_1d)  :: self
        real(dl), intent(in)                              :: x
        integer , intent(in), optional                    :: index
        real(dl), intent(in), optional                    :: coeff
        real(dl) :: equispacedlinearintepolatefunction1dfirstderivative

        integer  :: ind
        real(dl) :: x1, x2, y1, y2, mu

        ! initialize to null value:
        equispacedlinearintepolatefunction1dfirstderivative = self%null_value
        if ( self%has_null_value ) then
            ! if outside the interpolation range return the null value:
            if ( x <= self%x_initial .or. x >= self%x_final ) return
        else
            ! if below the interpolation range return the first value:
            if ( x <= self%x_initial ) then
                equispacedlinearintepolatefunction1dfirstderivative = self%yp(1)
                return
            end if
            ! if above the interpolation range return the first value:
            if ( x >= self%x_final   ) then
                equispacedlinearintepolatefunction1dfirstderivative = self%yp(self%num_points)
                return
            end if
        end if

        ! return the index of the point:
        if ( present(index) ) then
            ind = index
        else
            ind = int( ( x-self%x_initial)/self%grid_width ) +1
        end if

        ! get the interpolation coefficient:
        if ( present(coeff) ) then
            mu = coeff
        else
            ! store the x values:
            x1  = self%x(ind)
            x2  = self%x(ind+1)
            ! compute the linear interpolation coefficient:
            mu  = (x-x1)/(x2-x1)
        end if

        ! store the y values:
        y1  = self%yp(ind)
        y2  = self%yp(ind+1)

        ! compute the linear interpolation:
        equispacedlinearintepolatefunction1dfirstderivative = y1*( 1._dl -mu ) +y2*mu

    end function equispacedlinearintepolatefunction1dfirstderivative

    ! ---------------------------------------------------------------------------------------------
    function equispacedlinearintepolatefunction1dsecondderivative( self, x, index, coeff )

        implicit none

        class(equispaced_linear_interpolate_function_1d)  :: self
        real(dl), intent(in)                              :: x
        integer , intent(in), optional                    :: index
        real(dl), intent(in), optional                    :: coeff
        real(dl) :: equispacedlinearintepolatefunction1dsecondderivative

        integer  :: ind
        real(dl) :: x1, x2, y1, y2, mu

        ! initialize to null value:
        equispacedlinearintepolatefunction1dsecondderivative = self%null_value
        if ( self%has_null_value ) then
            ! if outside the interpolation range return the null value:
            if ( x <= self%x_initial .or. x >= self%x_final ) return
        else
            ! if below the interpolation range return the first value:
            if ( x <= self%x_initial ) then
                equispacedlinearintepolatefunction1dsecondderivative = self%ypp(1)
                return
            end if
            ! if above the interpolation range return the first value:
            if ( x >= self%x_final   ) then
                equispacedlinearintepolatefunction1dsecondderivative = self%ypp(self%num_points)
                return
            end if
        end if

        ! return the index of the point:
        if ( present(index) ) then
            ind = index
        else
            ind = int( ( x-self%x_initial)/self%grid_width ) +1
        end if

        ! get the interpolation coefficient:
        if ( present(coeff) ) then
            mu = coeff
        else
            ! store the x values:
            x1  = self%x(ind)
            x2  = self%x(ind+1)
            ! compute the linear interpolation coefficient:
            mu  = (x-x1)/(x2-x1)
        end if

        ! store the y values:
        y1  = self%ypp(ind)
        y2  = self%ypp(ind+1)

        ! compute the linear interpolation:
        equispacedlinearintepolatefunction1dsecondderivative = y1*( 1._dl -mu ) +y2*mu

    end function equispacedlinearintepolatefunction1dsecondderivative

    ! ---------------------------------------------------------------------------------------------
    function equispacedlinearintepolatefunction1dthirdderivative( self, x, index, coeff )

        implicit none

        class(equispaced_linear_interpolate_function_1d)  :: self
        real(dl), intent(in)                              :: x
        integer , intent(in), optional                    :: index
        real(dl), intent(in), optional                    :: coeff
        real(dl) :: equispacedlinearintepolatefunction1dthirdderivative

        integer  :: ind
        real(dl) :: x1, x2, y1, y2, mu

        ! initialize to null value:
        equispacedlinearintepolatefunction1dthirdderivative = self%null_value
        if ( self%has_null_value ) then
            ! if outside the interpolation range return the null value:
            if ( x <= self%x_initial .or. x >= self%x_final ) return
        else
            ! if below the interpolation range return the first value:
            if ( x <= self%x_initial ) then
                equispacedlinearintepolatefunction1dthirdderivative = self%yppp(1)
                return
            end if
            ! if above the interpolation range return the first value:
            if ( x >= self%x_final   ) then
                equispacedlinearintepolatefunction1dthirdderivative = self%yppp(self%num_points)
                return
            end if
        end if

        ! return the index of the point:
        if ( present(index) ) then
            ind = index
        else
            ind = int( ( x-self%x_initial)/self%grid_width ) +1
        end if

        ! get the interpolation coefficient:
        if ( present(coeff) ) then
            mu = coeff
        else
            ! store the x values:
            x1  = self%x(ind)
            x2  = self%x(ind+1)
            ! compute the linear interpolation coefficient:
            mu  = (x-x1)/(x2-x1)
        end if

        ! store the y values:
        y1  = self%yppp(ind)
        y2  = self%yppp(ind+1)

        ! compute the linear interpolation:
        equispacedlinearintepolatefunction1dthirdderivative = y1*( 1._dl -mu ) +y2*mu

    end function equispacedlinearintepolatefunction1dthirdderivative

    ! ---------------------------------------------------------------------------------------------
    function equispacedlinearintepolatefunction1dintegral( self, x, index, coeff )

        implicit none

        class(equispaced_linear_interpolate_function_1d)  :: self
        real(dl), intent(in)                              :: x
        integer , intent(in), optional                    :: index
        real(dl), intent(in), optional                    :: coeff
        real(dl) :: equispacedlinearintepolatefunction1dintegral

        integer  :: ind
        real(dl) :: x1, x2, y1, y2, mu

        ! initialize to null value:
        equispacedlinearintepolatefunction1dintegral = self%null_value
        if ( self%has_null_value ) then
            ! if outside the interpolation range return the null value:
            if ( x <= self%x_initial .or. x >= self%x_final ) return
        else
            ! if below the interpolation range return the first value:
            if ( x <= self%x_initial ) then
                equispacedlinearintepolatefunction1dintegral = self%yint(1)
                return
            end if
            ! if above the interpolation range return the first value:
            if ( x >= self%x_final   ) then
                equispacedlinearintepolatefunction1dintegral = self%yint(self%num_points)
                return
            end if
        end if

        ! return the index of the point:
        if ( present(index) ) then
            ind = index
        else
            ind = int( ( x-self%x_initial)/self%grid_width ) +1
        end if

        ! get the interpolation coefficient:
        if ( present(coeff) ) then
            mu = coeff
        else
            ! store the x values:
            x1  = self%x(ind)
            x2  = self%x(ind+1)
            ! compute the linear interpolation coefficient:
            mu  = (x-x1)/(x2-x1)
        end if

        ! store the y values:
        y1  = self%yint(ind)
        y2  = self%yint(ind+1)

        ! compute the linear interpolation:
        equispacedlinearintepolatefunction1dintegral = y1*( 1._dl -mu ) +y2*mu

    end function equispacedlinearintepolatefunction1dintegral

    ! ---------------------------------------------------------------------------------------------
    subroutine equispacedlinearintepolatefunction1dinitderivatives( self )

        implicit none

        class(equispaced_linear_interpolate_function_1d)  :: self

        write(*,*) 'ERROR: not yet implemented (IW)'
        stop

    end subroutine equispacedlinearintepolatefunction1dinitderivatives

    ! ---------------------------------------------------------------------------------------------

end module equispaced_linear_interpolation_1d

!----------------------------------------------------------------------------------------