splines.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
!/*----------------------------------------------------------------------------
! This module contains a data type (spline_t) to contain 1D functions,
! along with a series of procedures to manage them (to define them, to obtain
! its values, to operate on them, etc). The internal representation of the
! functions is done through cubic splines, handled by the GSL library. For
! the user of the module, this internal representation is hidden; one just
! works with what are called hereafter "spline functions".
!
! To define a function, one must supply a set {x(i),y(i)} of pairs of values
! -- the abscissa and the value of the function.
!
! [*] DATA TYPE:
!     To define a spline function:
!     type(spline_t) :: f
!
! [1] INITIALIZATION:
!     Before using any function, one should initialize it:
!
!     Interface:
!     subroutine spline_init(spl)
!        type(spline_t), intent(out) :: spl [or spl(:) or spl(:, :)]
!     end subroutine spline_init
!
!     Usage:
!     call spline_init(f)
!
! [2] FINALIZATION:
!     To empty any allocated space, one should finalize the function:
!
!     Interface:
!     subroutine spline_end(spl)
!        type(spline_t), intent(inout) :: spl [or spl(:) or spl(:, :)]
!     end subroutine spline_end
!
!     Usage
!     type(spline_t) :: f
!     call spline_end(f)
!
! [3] TO DEFINE A FUNCTION:
!     To "fill" an initialized function f,use spline_fit
!
!     Interface:
!     subroutine spline_fit(n, x, y, spl, threshold)
!        integer, intent(in) :: nrc
!        real(X), intent(in) :: x(n), y(n)
!        type(spline_t), intent(out) :: spl
!        real(real64), optional, intent(in) :: threshold
!     end subroutine spline_fit
!
!     (X may be 4 or eight, for single or double precision)
!
!     Usage:
!     call spline_fit(n, x, y, f)
!     n is the number of values that are supplied, x the abscissas, and y
!     the value of the function to represent at each point.
!     The optional argument threshold is used to find the radius associated with the spline.
!
! [4] FUNCTION VALUES:
!     To retrieve the value of a function at a given point:
!
!     Interface:
!
!     real(real64) function spline_eval(spl, x)
!       type(spline_t), intent(in) :: spl
!       real(real64), intent(in) :: x
!     end function spline_eval
!
!     Usage:
!     To retrieve the value of function f at point x, and place it into
!     real value val:
!     val = spline_eval(f, x)
!
! [5] SUM:
!     If you have two defined functions, f and g, and an initialized function
!     h, you may sum f and g and place the result onto g. For this purpose,
!     the grid that defined the first operand, f, will be used -- the values
!     of g will be interpolated to get the sum and place it in h.
!
!     Interface:
!     subroutine spline_sum(spl1, spl2, splsum, threshold)
!       type(spline_t), intent(in)  :: spl1, spl2
!       type(spline_t), intent(out) :: splsum
!       real(r8), optional, intent(in) :: threshold
!     end subroutine spline_sum
!
!     Usage:
!     call spline_init(f)
!     call spline_init(g)
!     call spline_init(h)
!     call spline_fit(npointsf, xf, yf, f)
!     call spline_fit(npointsg, xg, yg, g)
!     call spline_sum(f, g, h)
!
! [6] MULTIPLICATION BY A SCALAR
!     You may multiply a given spline-represented spline by a real number:
!
!     Interface:
!     subroutine spline_times(a, spl)
!       type(spline_t), intent(inout)  :: spl
!       real(real64), intent(in) :: a
!       real(real64), optional, intent(in) :: threshold
!     end subroutine spline_times
!
!     Usage:
!     call spline_init(f)
!     call spline_fit(npoints, x, y, f) ! Fill f with y values at x points
!     call spline_times(a, f) ! Now f contains a*y values at x points.
!
! [7] INTEGRAL:
!     Given a defined function, the function spline_integral returns its
!     integral. The interval of integration may or may not be supplied.
!
!     Interface:
!     real(real64) function spline_integral(spl [,a,b])
!       type(spline_integral), intent(in) :: spl
!       real(real64), intent(in), optional :: a, b
!     end function spline_integral
!
! [8] DOT PRODUCT:
!     Given two defined functions f and g, the function spline_dotp returns
!     the value of their dot-product: int {dx f(x)*g(x)}. The mesh used to do
!     so the mesh of the fist-function (note that as a result the definition
!     is no longer conmutative).
!
!     Interface:
!     real(real64) function spline_dotp(spl1, spl2)
!       type(spline_t), intent(in) :: spl1, spl2
!     end function spline_dotp
!
! Note: The following routines, spline_3dft, spline_cut and spline_filter,
! assume that the spline functions are the radial part of a 3 dimensional function with
! spherical symmetry. This is why the Fourier transform of F(\vec{r}) = f(r), is:
!       F(\vec{g}) = f(g) = \frac{4\pi}{g} \int_{0}^{\infty} { r*sin(g*r)*f(r) }
! which coincides with the inverse Fourier transform, except that the inverse Fourier
! transform should be multiplied by a (2*\pi)^{-3} factor.
!
! [9] FOURIER TRANSFORM:
!     If a spline function f is representing the radial part of a spherically
!     symmetric function F(\vec{r}), its Fourier transform is:
!       F(\vec{g}) = f(g) = \frac{4\pi}{g} \int_{0}^{\infty} { r*sin(g*r)*f(r) }
!     It is assumed that f is defined in some interval [0,rmax], and that it is
!     null at rmax and beyond. One may obtain f(g) by using spline_3dft.
!     The result is placed on the spline data type splw. This has to be initialized,
!     and may or may not be filled. If it is filled, the abscissas that it has
!     are used to define the function. If it is not filled, an equally spaced
!     grid is constructed to define the function, in the interval [0, gmax], where
!     gmax has to be supplied by the caller.
!
!     Interface:
!     subroutine spline_3dft(spl, splw, gmax)
!       type(spline_t), intent(in)    :: spl
!       type(spline_t), intent(inout) :: splw
!       real(real64), intent(in), optional :: gmax
!     end subroutine spline_3dft
!
! [10] BESSEL TRANSFORM:
!
!
! [11] CUTTING A FUNCTION:
!     spline_cut multiplies a given function by a cutoff-function, which
!     is defined to be one in [0, cutoff], and \exp\{-beta*(x/cutoff-1)^2\}
!
!     Interface:
!     subroutine spline_cut(spl, cutoff, beta)
!       type(spline_t), intent(in) :: spl
!       real(real64), intent(in) :: cutoff, beta
!     end subroutine spline_cut
!
! [13] PRINTING A FUNCTION:
!     It prints to a file the (x,y) values that were used to define a function.
!     The file is pointed to by its Fortran unit given by argument iunit.
!
!     Interface:
!     subroutine spline_print(spl, iunit)
!       type(spline_t), intent(in) :: spl [ or spl(:) or spl(:, :)]
!       integer, intent(in) :: iunit
!     end subroutine spline_print
!
! [14] DIFFERENTIATE A FUNCTION:
!
!----------------------------------------------------------------------------*/!
module splines_oct_m
  use debug_oct_m
  use global_oct_m
  use iso_c_binding
  use, intrinsic :: iso_fortran_env
  use loct_math_oct_m
  use messages_oct_m
  use profiling_oct_m
 
  implicit none
 
  ! Define which routines can be seen from the outside.
  private
  public ::               &
    spline_t,        & ! [*]
    spline_init,     & ! [1]
    spline_end,      & ! [2]
    spline_fit,      & ! [3]
    spline_eval,     & ! [4]
    spline_eval_vec, & ! [4]
    spline_sum,      & ! [5]
    spline_times,    & ! [6]
    spline_integral, & ! [7]
    spline_dotp,     & ! [8]
    spline_3dft,     & ! [9]
    spline_besselft, & ! [10]
    spline_cut,      & ! [11]
    spline_print,    & ! [13]
    spline_der,      &
    spline_der2,     &
    spline_div,      &
    spline_mult,     &
    spline_force_pos, &
    spline_range_min, &
    spline_range_max, &
    spline_x_threshold, &
    spline_generate_shifted_grid
 
  type spline_t
    private
    real(real64)     :: x_limit(2)
    type(c_ptr)      :: spl, acc
    real(real64), public :: x_threshold
  end type spline_t
 
  interface spline_eval_vec
    module procedure spline_eval8_array
    module procedure spline_evalz_array
  end interface spline_eval_vec
 
  interface spline_integral
    module procedure spline_integral_full
    module procedure spline_integral_limits
  end interface spline_integral
 
  interface spline_init
    module procedure spline_init_0
    module procedure spline_init_1
    module procedure spline_init_2
  end interface spline_init
 
  interface spline_end
    module procedure spline_end_0
    module procedure spline_end_1
    module procedure spline_end_2
  end interface spline_end
 
  interface spline_print
    module procedure spline_print_0
    module procedure spline_print_1
    module procedure spline_print_2
  end interface spline_print
 
  ! These are interfaces to functions defined in oct_gsl_f.c, which in turn
  ! take care of calling the GSL library.
  interface
    subroutine oct_spline_end(spl, acc)
      use iso_c_binding
      implicit none
      type(c_ptr), intent(inout) :: spl, acc
    end subroutine oct_spline_end
 
    subroutine oct_spline_fit(nrc, x, y, spl, acc)
      use iso_c_binding
      use, intrinsic :: iso_fortran_env
      implicit none
      integer,     intent(in) :: nrc
      real(real64),    intent(in) :: x
      real(real64),    intent(in) :: y
      type(c_ptr), intent(inout) :: spl
      type(c_ptr), intent(inout) :: acc
    end subroutine oct_spline_fit
 
    real(real64) function oct_spline_eval(x, spl, acc)
      use iso_c_binding
      use, intrinsic :: iso_fortran_env
      implicit none
      real(real64),    intent(in) :: x
      type(c_ptr), intent(in) :: spl
      type(c_ptr), intent(in) :: acc
    end function oct_spline_eval
 
    subroutine oct_spline_eval_array(nn, xf, spl, acc)
      use iso_c_binding
      use, intrinsic :: iso_fortran_env
      implicit none
      integer,     intent(in)    :: nn
      real(real64),    intent(inout) :: xf
      type(c_ptr), intent(in) :: spl
      type(c_ptr), intent(in) :: acc
    end subroutine oct_spline_eval_array
 
    subroutine oct_spline_eval_arrayz(nn, xf, spl, acc)
      use iso_c_binding
      use, intrinsic :: iso_fortran_env
      implicit none
      integer,     intent(in)    :: nn
      complex(real64), intent(inout) :: xf
      type(c_ptr), intent(in) :: spl
      type(c_ptr), intent(in) :: acc
    end subroutine oct_spline_eval_arrayz
 
    real(real64) function oct_spline_eval_der(x, spl, acc)
      use iso_c_binding
      use, intrinsic :: iso_fortran_env
      implicit none
      real(real64),    intent(in) :: x
      type(c_ptr), intent(in) :: spl
      type(c_ptr), intent(in) :: acc
    end function oct_spline_eval_der
 
    real(real64) function oct_spline_eval_der2(x, spl, acc)
      use iso_c_binding
      use, intrinsic :: iso_fortran_env
      implicit none
      real(real64),    intent(in) :: x
      type(c_ptr), intent(in) :: spl
      type(c_ptr), intent(in) :: acc
    end function oct_spline_eval_der2
 
    integer pure function oct_spline_npoints(spl, acc)
      use iso_c_binding
      implicit none
      type(c_ptr), intent(in) :: spl
      type(c_ptr), intent(in) :: acc
    end function oct_spline_npoints
 
    pure subroutine oct_spline_x(spl, x)
      use iso_c_binding
      use, intrinsic :: iso_fortran_env
      implicit none
      type(c_ptr), intent(in)  :: spl
      real(real64),    intent(out) :: x
    end subroutine oct_spline_x
 
    subroutine oct_spline_y(spl, y)
      use iso_c_binding
      use, intrinsic :: iso_fortran_env
      implicit none
      type(c_ptr), intent(in) :: spl
      real(real64),    intent(out) :: y
    end subroutine oct_spline_y
 
    real(real64) pure function oct_spline_eval_integ(spl, a, b, acc)
      use iso_c_binding
      use, intrinsic :: iso_fortran_env
      implicit none
      type(c_ptr), intent(in) :: spl
      real(real64),    intent(in) :: a
      real(real64),    intent(in) :: b
      type(c_ptr), intent(in) :: acc
    end function oct_spline_eval_integ
 
    real(real64) pure function oct_spline_eval_integ_full(spl, acc)
      use iso_c_binding
      use, intrinsic :: iso_fortran_env
      implicit none
      type(c_ptr), intent(in) :: spl
      type(c_ptr), intent(in) :: acc
    end function oct_spline_eval_integ_full
  end interface
 
  real(real64), parameter :: tol_x = 1e-14_real64 
 
contains
 
  !------------------------------------------------------------
  subroutine spline_init_0(spl)
    type(spline_t), intent(out) :: spl
 
    ! No PUSH SUB, called too often
 
    spl%spl = c_null_ptr
    spl%acc = c_null_ptr
 
    ! deliberately illegal values, for checking
    spl%x_limit(1) = -1_real64
    spl%x_limit(2) = -2_real64
 
    spl%x_threshold = m_zero
 
  end subroutine spline_init_0
 
 
  !------------------------------------------------------------
  subroutine spline_init_1(spl)
    type(spline_t), intent(out) :: spl(:)
 
    integer :: i
 
    push_sub(spline_init_1)
 
    do i = 1, size(spl)
      call spline_init_0(spl(i))
    end do
 
    pop_sub(spline_init_1)
  end subroutine spline_init_1
 
 
  !------------------------------------------------------------
  subroutine spline_init_2(spl)
    type(spline_t), intent(out) :: spl(:, :)
 
    integer :: i, j
 
    push_sub(spline_init_2)
 
    do i = 1, size(spl, 1)
      do j = 1, size(spl, 2)
        call spline_init_0(spl(i, j))
      end do
    end do
 
    pop_sub(spline_init_2)
  end subroutine spline_init_2
 
 
  !------------------------------------------------------------
  subroutine spline_end_0(spl)
    type(spline_t), intent(inout) :: spl
 
    ! No PUSH SUB, called too often.
 
    if (c_associated(spl%spl) .and. c_associated(spl%acc)) then
      call oct_spline_end(spl%spl, spl%acc)
    end if
    spl%spl = c_null_ptr
    spl%acc = c_null_ptr
 
  end subroutine spline_end_0
 
 
  !------------------------------------------------------------
  subroutine spline_end_1(spl)
    type(spline_t), intent(inout) :: spl(:)
 
    integer :: i
 
    push_sub(spline_end_1)
 
    do i = 1, size(spl)
      call spline_end_0(spl(i))
    end do
 
    pop_sub(spline_end_1)
  end subroutine spline_end_1
 
 
  !------------------------------------------------------------
  subroutine spline_end_2(spl)
    type(spline_t), intent(inout) :: spl(:, :)
 
    integer :: i, j
 
    push_sub(spline_end_2)
 
    do i = 1, size(spl, 1)
      do j = 1, size(spl, 2)
        call spline_end_0(spl(i, j))
      end do
    end do
 
    pop_sub(spline_end_2)
  end subroutine spline_end_2
 
 
  !------------------------------------------------------------
  subroutine spline_fit(nrc, rofi, ffit, spl, threshold)
    integer,                intent(in)    :: nrc
    real(real64),           intent(in)    :: rofi(:)
    real(real64),           intent(in)    :: ffit(:)
    type(spline_t),         intent(inout) :: spl
    real(real64), optional, intent(in)    :: threshold
 
    !No PUSH SUB, called too often
 
    assert(nrc > 0)
 
    spl%x_limit(1) = rofi(1)
    spl%x_limit(2) = rofi(nrc)
    call oct_spline_fit(nrc, rofi(1), ffit(1), spl%spl, spl%acc)
 
    if (present(threshold)) then
      if (threshold > m_zero) then
        spl%x_threshold = spline_x_threshold(spl, threshold)
      else
        spl%x_threshold = spl%x_limit(2)
      end if
    else
      spl%x_threshold = m_zero
    end if
 
  end subroutine spline_fit
 
  !------------------------------------------------------------
  real(real64) function spline_eval(spl, x)
    type(spline_t),     intent(in) :: spl
    real(real64),       intent(in) :: x
 
    spline_eval = oct_spline_eval(x, spl%spl, spl%acc)
  end function spline_eval
 
 
  !------------------------------------------------------------
  subroutine spline_eval8_array(spl, nn, xf)
    type(spline_t), intent(in)    :: spl
    integer,        intent(in)    :: nn
    real(real64),   intent(inout) :: xf(:)
 
    assert(nn > 0)
 
    call oct_spline_eval_array(nn, xf(1), spl%spl, spl%acc)
  end subroutine spline_eval8_array
 
 
  !------------------------------------------------------------
  subroutine spline_evalz_array(spl, nn, xf)
    type(spline_t),  intent(in)    :: spl
    integer,         intent(in)    :: nn
    complex(real64), intent(inout) :: xf(:)
 
    assert(nn > 0)
 
    call oct_spline_eval_arrayz(nn, xf(1), spl%spl, spl%acc)
  end subroutine spline_evalz_array
 
  !------------------------------------------------------------
  subroutine spline_sum(spl1, spl2, splsum, threshold)
    type(spline_t), intent(in)  :: spl1
    type(spline_t), intent(in)  :: spl2
    type(spline_t), intent(out) :: splsum
    real(real64),   intent(in)  :: threshold
 
    integer :: npoints, i
    real(real64), allocatable :: x(:), y(:), y2(:)
 
    push_sub(spline_sum)
 
    npoints = oct_spline_npoints(spl1%spl, spl1%acc)
    assert(npoints > 0)
 
    safe_allocate( x(1:npoints))
    safe_allocate( y(1:npoints))
    safe_allocate(y2(1:npoints))
 
    call oct_spline_x(spl1%spl, x(1))
    call oct_spline_y(spl1%spl, y(1))
 
    !$omp parallel do
    do i = 1, npoints
      y2(i) = spline_eval(spl2, x(i)) + y(i)
    end do
 
    call spline_fit(npoints, x, y2, splsum, threshold)
 
    safe_deallocate_a(x)
    safe_deallocate_a(y)
    safe_deallocate_a(y2)
 
    pop_sub(spline_sum)
  end subroutine spline_sum
 
 
  !------------------------------------------------------------
  subroutine spline_times(a, spl, threshold)
    real(real64),   intent(in)     :: a
    type(spline_t), intent(inout)  :: spl
    real(real64),   intent(in)     :: threshold
 
    integer :: npoints, i
    real(real64), allocatable :: x(:), y(:)
 
    push_sub(spline_times)
 
    npoints = oct_spline_npoints(spl%spl, spl%acc)
    assert(npoints > 0)
 
    safe_allocate(x(1:npoints))
    safe_allocate(y(1:npoints))
 
    call oct_spline_x(spl%spl, x(1))
    call oct_spline_y(spl%spl, y(1))
    call oct_spline_end(spl%spl, spl%acc)
    !$omp parallel do simd
    do i = 1, npoints
      y(i) = a*y(i)
    end do
    call spline_fit(npoints, x, y, spl, threshold)
 
    safe_deallocate_a(x)
    safe_deallocate_a(y)
 
    pop_sub(spline_times)
  end subroutine spline_times
 
 
  !------------------------------------------------------------
  real(real64) function spline_integral_full(spl) result(res)
    type(spline_t), intent(in) :: spl
 
    push_sub(spline_integral_full)
 
    res = oct_spline_eval_integ_full(spl%spl, spl%acc)
 
    pop_sub(spline_integral_full)
  end function spline_integral_full
 
 
  !------------------------------------------------------------
  real(real64) pure function spline_integral_limits(spl, a, b) result(res)
    type(spline_t), intent(in) :: spl
    real(real64),   intent(in) :: a
    real(real64),   intent(in) :: b
 
    res = oct_spline_eval_integ(spl%spl, a, b, spl%acc)
  end function spline_integral_limits
 
 
  !------------------------------------------------------------
  real(real64) function spline_dotp(spl1, spl2) result (res)
    type(spline_t), intent(in) :: spl1
    type(spline_t), intent(in) :: spl2
 
    type(spline_t) :: aux
    integer :: npoints, i
    real(real64), allocatable :: x(:), y(:)
 
    push_sub(spline_dotp)
 
    npoints = oct_spline_npoints(spl1%spl, spl1%acc)
    assert(npoints > 0)
 
    safe_allocate(x(1:npoints))
    safe_allocate(y(1:npoints))
 
    call oct_spline_x(spl1%spl, x(1))
    call oct_spline_y(spl1%spl, y(1))
    !$omp parallel do
    do i = 1, npoints
      y(i) = y(i)*oct_spline_eval(x(i), spl2%spl, spl2%acc)
    end do
    call spline_init(aux)
    call spline_fit(npoints, x, y, aux, m_zero)
    res = oct_spline_eval_integ(aux%spl, x(1), x(npoints), aux%acc)
 
    safe_deallocate_a(x)
    safe_deallocate_a(y)
 
    pop_sub(spline_dotp)
  end function spline_dotp
 
 
  !------------------------------------------------------------
  subroutine spline_3dft(spl, splw, threshold, gmax)
    type(spline_t),      intent(in)    :: spl
    type(spline_t),      intent(inout) :: splw
    real(real64),        intent(in)    :: threshold
    real(real64), optional,     intent(in)    :: gmax
 
    type(spline_t) :: aux
    real(real64) :: g, dg
    integer :: np
    integer :: npoints, i, j
    real(real64), allocatable :: x(:), y(:), y2(:), xw(:), yw(:)
 
    push_sub(spline_3dft)
 
    npoints = oct_spline_npoints(spl%spl, spl%acc)
    assert(npoints > 0)
 
    safe_allocate( x(1:npoints))
    safe_allocate( y(1:npoints))
    safe_allocate(y2(1:npoints))
 
    call oct_spline_x(spl%spl, x(1))
    call oct_spline_y(spl%spl, y(1))
 
    ! Check whether splw comes with a defined grid, or else build it.
    if (c_associated(splw%spl)) then
      np = oct_spline_npoints(splw%spl, splw%acc)
      safe_allocate(xw(1:np))
      safe_allocate(yw(1:np))
      call oct_spline_x(splw%spl, xw(1))
      ! But now we need to kill the input:
      call spline_end(splw)
    else
      np = 200 ! hard coded value
      dg = gmax/(np-1)
      safe_allocate(xw(1:np))
      safe_allocate(yw(1:np))
      do i = 1, np
        g = (i-1)*dg
        xw(i) = g
      end do
    end if
 
    ! The first point, xw(1) = 0.0 and it has to be treated separately.
    !$omp parallel do simd
    do j = 1, npoints
      y2(j) = m_four*m_pi*y(j)*x(j)**2
    end do
    call spline_init(aux)
    call spline_fit(npoints, x, y2, aux, m_zero)
    yw(1) = oct_spline_eval_integ_full(aux%spl, aux%acc)
    call spline_end(aux)
 
    do i = 2, np
      !$omp parallel do simd
      do j = 1, npoints
        y2(j) = (m_four*m_pi/xw(i))*y(j)*x(j)*sin(xw(i)*x(j))
      end do
      call spline_init(aux)
      call spline_fit(npoints, x, y2, aux, m_zero)
      yw(i) = oct_spline_eval_integ_full(aux%spl, aux%acc)
      call spline_end(aux)
    end do
 
    call spline_init(splw)
    call spline_fit(np, xw, yw, splw, threshold)
 
    safe_deallocate_a(x)
    safe_deallocate_a(y)
    safe_deallocate_a(y2)
    safe_deallocate_a(xw)
    safe_deallocate_a(yw)
 
    pop_sub(spline_3dft)
  end subroutine spline_3dft
 
 
  !------------------------------------------------------------
  subroutine spline_besselft(spl, splw, l, threshold, gmax)
    type(spline_t),    intent(in)    :: spl
    type(spline_t),    intent(inout) :: splw
    integer,           intent(in)    :: l
    real(real64),   optional, intent(in)    :: threshold
    real(real64),   optional, intent(in)    :: gmax
 
    type(spline_t) :: aux
    real(real64) :: dg
    integer :: np
    integer :: npoints, i, j
    real(real64), allocatable :: x(:), y(:), y2(:), xw(:), yw(:)
 
    push_sub(spline_besselft)
 
    npoints = oct_spline_npoints(spl%spl, spl%acc)
    assert(npoints > 0)
 
    safe_allocate( x(1:npoints))
    safe_allocate( y(1:npoints))
    safe_allocate(y2(1:npoints))
 
    call oct_spline_x(spl%spl, x(1))
    call oct_spline_y(spl%spl, y(1))
 
    ! Check whether splw comes with a defined grid, or else build it.
    if (c_associated(splw%spl)) then
      np = oct_spline_npoints(splw%spl, splw%acc)
      safe_allocate(xw(1:np))
      safe_allocate(yw(1:np))
      call oct_spline_x(splw%spl, xw(1))
      ! But now we need to kill the input:
      call spline_end(splw)
    else
      assert(present(gmax))
      np = 1000 ! hard coded value
      dg = gmax/(np-1)
      safe_allocate(xw(1:np))
      safe_allocate(yw(1:np))
      !$omp parallel do
      do i = 1, np
        xw(i) = (i-1)*dg
      end do
    end if
 
    do i = 1, np
      !$omp parallel do
      do j = 1, npoints
        y2(j) = y(j) * x(j)**2 * loct_sph_bessel(l, x(j)*xw(i))
      end do
      call spline_init(aux)
      call spline_fit(npoints, x, y2, aux)
      yw(i) = sqrt(m_two/m_pi)*oct_spline_eval_integ_full(aux%spl, aux%acc)
 
      call spline_end(aux)
    end do
 
    call spline_init(splw)
    call spline_fit(np, xw, yw, splw, threshold)
 
    safe_deallocate_a(x)
    safe_deallocate_a(y)
    safe_deallocate_a(y2)
    safe_deallocate_a(xw)
    safe_deallocate_a(yw)
 
    pop_sub(spline_besselft)
  end subroutine spline_besselft
 
 
  !------------------------------------------------------------
  subroutine spline_cut(spl, cutoff, beta, threshold)
    type(spline_t), intent(inout) :: spl
    real(real64),   intent(in)    :: cutoff
    real(real64),   intent(in)    :: beta
    real(real64), optional,intent(in)    :: threshold
 
    integer :: npoints, i
    real(real64), allocatable :: x(:), y(:)
    real(real64) :: exp_arg
 
    push_sub(spline_cut)
 
    npoints = oct_spline_npoints(spl%spl, spl%acc)
    assert(npoints > 0)
 
    safe_allocate(x(1:npoints))
    safe_allocate(y(1:npoints))
 
    call oct_spline_x(spl%spl, x(1))
    call oct_spline_y(spl%spl, y(1))
    call oct_spline_end(spl%spl, spl%acc)
 
    do i = npoints, 1, -1
      if (x(i)<cutoff) then
        exit
      end if
 
      exp_arg = -beta*(x(i)/cutoff - m_one)**2
      !To avoid underflows
      if (exp_arg > m_min_exp_arg) then
        y(i) = y(i) * exp(exp_arg)
      else
        y(i) = m_zero
      end if
    end do
    call spline_fit(npoints, x, y, spl, threshold)
 
    safe_deallocate_a(x)
    safe_deallocate_a(y)
 
    pop_sub(spline_cut)
  end subroutine spline_cut
 
 
  !------------------------------------------------------------
  subroutine spline_div(spla, splb, threshold)
    type(spline_t),   intent(inout) :: spla
    type(spline_t),   intent(in)    :: splb
    real(real64), optional,  intent(in)    :: threshold
 
    integer :: npoints, i, new_npoints
    real(real64), allocatable :: x(:), y(:)
    real(real64) :: aa
 
    push_sub(spline_div)
 
    npoints = oct_spline_npoints(spla%spl, spla%acc)
    assert(npoints > 0)
 
    safe_allocate(x(1:npoints))
    safe_allocate(y(1:npoints))
 
    call oct_spline_x(spla%spl, x(1))
    call oct_spline_y(spla%spl, y(1))
    call oct_spline_end(spla%spl, spla%acc)
 
    assert(splb%x_limit(2) >= splb%x_limit(1))
 
    ! As an output, the new spline only contains the new number of points
    ! that goes up to the cutoff of the splb spline
    new_npoints = 0
    do i = npoints, 1, -1
      if (x(i) > splb%x_limit(2)-tol_x) cycle
      aa = spline_eval(splb, x(i))
      y(i) = y(i)/aa
      new_npoints = new_npoints + 1
    end do
 
    call spline_fit(new_npoints, x, y, spla, threshold)
 
    safe_deallocate_a(x)
    safe_deallocate_a(y)
 
    pop_sub(spline_div)
  end subroutine spline_div
 
  !------------------------------------------------------------
  !Force all the values of the spline to be positive
  !------------------------------------------------------------
  subroutine spline_force_pos(spl, threshold)
    type(spline_t),   intent(inout) :: spl
    real(real64),     intent(in)    :: threshold
 
    integer :: npoints, i
    real(real64), allocatable :: x(:), y(:)
 
    push_sub(spline_force_pos)
 
    npoints = oct_spline_npoints(spl%spl, spl%acc)
    assert(npoints > 0)
 
    safe_allocate(x(1:npoints))
    safe_allocate(y(1:npoints))
 
    call oct_spline_x(spl%spl, x(1))
    call oct_spline_y(spl%spl, y(1))
    call oct_spline_end(spl%spl, spl%acc)
 
    !$omp parallel do simd
    do i = npoints, 1, -1
      y(i) = max(y(i),m_zero)
    end do
 
    call spline_fit(npoints, x, y, spl, threshold)
 
    safe_deallocate_a(x)
    safe_deallocate_a(y)
 
    pop_sub(spline_force_pos)
  end subroutine spline_force_pos
 
 
  !------------------------------------------------------------
  subroutine spline_mult(spla, splb, threshold)
    type(spline_t),  intent(inout) :: spla
    type(spline_t),  intent(in)    :: splb
    real(real64),    intent(in)    :: threshold
 
    integer :: npoints, i
    real(real64), allocatable :: x(:), y(:)
    real(real64) :: aa
 
    push_sub(spline_mult)
 
    npoints = oct_spline_npoints(spla%spl, spla%acc)
    if(npoints <= 0) then
      pop_sub(spline_mult)
      return
    end if
 
    safe_allocate(x(1:npoints))
    safe_allocate(y(1:npoints))
 
    call oct_spline_x(spla%spl, x(1))
    call oct_spline_y(spla%spl, y(1))
    call oct_spline_end(spla%spl, spla%acc)
 
    assert(splb%x_limit(2) >= splb%x_limit(1))
 
    !$omp parallel do simd private(aa)
    do i = npoints, 1, -1
      if (x(i) > splb%x_limit(2) - tol_x) then
        aa = m_zero
      else
        aa = spline_eval(splb, x(i))
      end if
      y(i) = y(i)*aa
    end do
 
    call spline_fit(npoints, x, y, spla, threshold)
 
    safe_deallocate_a(x)
    safe_deallocate_a(y)
 
    pop_sub(spline_mult)
  end subroutine spline_mult
 
 
  !------------------------------------------------------------
  subroutine spline_der(spl, dspl, threshold)
    type(spline_t), intent(in)    :: spl
    type(spline_t), intent(inout) :: dspl
    real(real64),   intent(in)    :: threshold
 
    integer :: npoints, i
    real(real64), allocatable :: x(:), y(:)
 
    push_sub(spline_der)
 
    ! Use the grid of dspl if it is present, otherwise use the same one of spl.
    if (.not. c_associated(dspl%spl)) then ! use the grid of spl
      npoints = oct_spline_npoints(spl%spl, spl%acc)
      assert(npoints > 0)
      safe_allocate(x(1:npoints))
      safe_allocate(y(1:npoints))
      call oct_spline_x(spl%spl, x(1))
    else ! use the grid of dspl
      npoints = oct_spline_npoints(dspl%spl, dspl%acc)
      assert(npoints > 0)
      safe_allocate(x(1:npoints))
      safe_allocate(y(1:npoints))
      call oct_spline_x(dspl%spl, x(1))
    end if
    do i = 1, npoints
      y(i) = oct_spline_eval_der(x(i), spl%spl, spl%acc)
    end do
    call spline_fit(npoints, x, y, dspl, threshold)
 
    safe_deallocate_a(x)
    safe_deallocate_a(y)
 
    pop_sub(spline_der)
  end subroutine spline_der
 
 
  !------------------------------------------------------------
  subroutine spline_der2(spl, dspl, threshold, grid)
    type(spline_t),      intent(in)    :: spl
    type(spline_t),      intent(inout) :: dspl
    real(real64),        intent(in)    :: threshold
    real(real64), optional, intent(in) :: grid(:)
 
    integer :: npoints, i
    real(real64), allocatable :: x(:), y(:)
 
    push_sub(spline_der2)
 
    ! Use the grid of dspl if it is present, otherwise use the same one of spl.
    if (.not. c_associated(dspl%spl)) then ! use the grid of spl
      npoints = oct_spline_npoints(spl%spl, spl%acc)
      assert(npoints > 0)
      safe_allocate(x(1:npoints))
      safe_allocate(y(1:npoints))
      if (present(grid)) then
        x = grid(:)
      else
        call oct_spline_x(spl%spl, x(1))
      end if
    else ! use the grid of dspl
      npoints = oct_spline_npoints(dspl%spl, dspl%acc)
      assert(npoints > 0)
      safe_allocate(x(1:npoints))
      safe_allocate(y(1:npoints))
      call oct_spline_x(dspl%spl, x(1))
    end if
    do i = 1, npoints
      y(i) = oct_spline_eval_der2(x(i), spl%spl, spl%acc)
    end do
    call spline_fit(npoints, x, y, dspl, threshold)
 
    safe_deallocate_a(x)
    safe_deallocate_a(y)
 
    pop_sub(spline_der2)
  end subroutine spline_der2
 
 
  !------------------------------------------------------------
  subroutine spline_print_0(spl, iunit)
    type(spline_t), intent(in) :: spl
    integer,        intent(in) :: iunit
 
    integer :: np, i
    real(real64), allocatable :: x(:), y(:)
 
    push_sub(spline_print_0)
 
    np = oct_spline_npoints(spl%spl, spl%acc)
    assert(np > 0)
    safe_allocate(x(1:np))
    safe_allocate(y(1:np))
 
    call oct_spline_x(spl%spl, x(1))
    call oct_spline_y(spl%spl, y(1))
    do i = 1, np
      write(iunit, '(2es18.8E3)') x(i), y(i)
    end do
 
    safe_deallocate_a(x)
    safe_deallocate_a(y)
 
    pop_sub(spline_print_0)
  end subroutine spline_print_0
 
 
  !------------------------------------------------------------
  subroutine spline_print_1(spl, iunit)
    type(spline_t), intent(in) :: spl(:)
    integer,        intent(in) :: iunit
 
    character(len=4)  :: fm
    integer :: np, i, n, j
    real(real64), allocatable :: x(:)
    real(real64) :: x_limit
 
    push_sub(spline_print_1)
 
    n = size(spl)
    if (n <= 0) then
      pop_sub(spline_print_1)
      return
    end if
 
    write(fm,'(i4)') n + 1
    fm = adjustl(fm)
    np = oct_spline_npoints(spl(1)%spl, spl(1)%acc)
    x_limit = spl(1)%x_limit(2)
    safe_allocate(x(1:np))
 
    call oct_spline_x(spl(1)%spl, x(1))
 
    ! Not all splines have the same number of points
    do i = 2, n
      np = min(np, oct_spline_npoints(spl(i)%spl, spl(i)%acc))
      x_limit = min(x_limit, spl(i)%x_limit(2))
    end do
 
    do i = 1, np
      if (x(i) > x_limit - tol_x) cycle
      write(iunit, '('//trim(fm)//'es18.8E3)') x(i), (spline_eval(spl(j), x(i)), j = 1, size(spl))
    end do
 
    safe_deallocate_a(x)
 
    pop_sub(spline_print_1)
  end subroutine spline_print_1
 
 
  !------------------------------------------------------------
  subroutine spline_print_2(spl, iunit)
    type(spline_t), intent(in) :: spl(:, :)
    integer,        intent(in) :: iunit
 
    character(len=4)  :: fm
    integer :: np, i, n1, n2, j, k
    real(real64), allocatable :: x(:)
    real(real64) :: x_limit
 
    push_sub(spline_print_2)
 
    n1 = size(spl, 1)
    n2 = size(spl, 2)
    if (n1 * n2 <= 0) then
      pop_sub(spline_print_2)
      return
    end if
 
    write(fm,'(i4)') n1*n2 + 1
    fm = adjustl(fm)
 
    np = oct_spline_npoints(spl(1, 1)%spl, spl(1, 1)%acc)
    safe_allocate(x(1:np))
 
    call oct_spline_x(spl(1, 1)%spl, x(1))
    x_limit = spl(1,1)%x_limit(2)
 
    ! Not all splines have the same number of points
    do i = 2, n1
      do j = 1, n2
        np = min(np, oct_spline_npoints(spl(i, j)%spl, spl(i, j)%acc))
        x_limit = min(x_limit, spl(i,j)%x_limit(2))
      end do
    end do
 
    do i = 1, np
      if (x(i) > x_limit - tol_x) cycle
      write(iunit, '('//trim(fm)//'es18.8E3)') x(i), &
        ((spline_eval(spl(j, k), x(i)), j = 1, n1), k = 1, n2)
    end do
 
    safe_deallocate_a(x)
 
    pop_sub(spline_print_2)
  end subroutine spline_print_2
 
 
  !------------------------------------------------------------
  real(real64) function spline_x_threshold(spl, threshold) result(r)
    type(spline_t), intent(in) :: spl
    real(real64),   intent(in) :: threshold
 
    integer :: ii, jj
    real(real64), parameter :: dx = 0.01_real64
 
    ! No PUSH SUB, called too often.
    assert(spl%x_limit(2) >= spl%x_limit(1))
 
    call profiling_in('SPLINE_CUTOFF_RADIUS')
 
    jj = floor(spl%x_limit(2)/dx) + 1
    ! Numerical errors can lead to jj to be overestimated by 1
    ! This test prevents this to occur
    do while(dx*(jj-1) > spl%x_limit(2) - m_half * dx)
      jj = jj-1
    end do
 
    do ii = jj, 1, -1
      r = dx*(ii-1)
      if (abs(spline_eval(spl, r)) > threshold) exit
    end do
 
    call profiling_out('SPLINE_CUTOFF_RADIUS')
 
  end function spline_x_threshold
 
  ! -------------------------------------------------------
 
  real(real64) pure function spline_range_min(this) result(range)
    type(spline_t), intent(in) :: this
 
    range = this%x_limit(1)
 
  end function spline_range_min
 
  ! -------------------------------------------------------
 
  real(real64) pure function spline_range_max(this) result(range)
    type(spline_t), intent(in) :: this
 
    range = this%x_limit(2)
 
  end function spline_range_max
 
  ! -------------------------------------------------------
  subroutine spline_generate_shifted_grid(this, x_shifted)
    type(spline_t),               intent(in) :: this
    real(real64), allocatable, intent(inout) :: x_shifted(:)
 
    real(real64), allocatable :: x(:)
    integer :: npoints, i
 
    push_sub(spline_generate_shifted_grid)
 
    ! We shift the grid by half of a spacing to force interpolation of the function
    ! This is useful for getting back the derivative on the same grid
    npoints = oct_spline_npoints(this%spl, this%acc)
    safe_allocate(x(1:npoints))
    safe_allocate(x_shifted(1:npoints))
    call oct_spline_x(this%spl, x(1))
    do i = 1, npoints -1
      x_shifted(i) = m_half*(x(i) + x(i+1))
    end do
    x_shifted(npoints) = x(npoints)
    safe_deallocate_a(x)
 
    pop_sub(spline_generate_shifted_grid)
  end subroutine spline_generate_shifted_grid
 
end module splines_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: