! $Id$
! -----------------------------------------------------------
! Module spline_module
! -----------------------------------------------------------
! Code area 9: general
! -----------------------------------------------------------

!!****h* Conquest/spline_module
!!  NAME
!! 
!!  PURPOSE
!! 
!!  AUTHOR
!!   D.R.Bowler
!!  CREATION DATE
!!   23/05/01
!!  MODIFICATION HISTORY
!!   31/05/2001 dave
!!    Added RCS Id tag to header, increased ROBODoc documentation
!!   20/06/2001 dave
!!    Added RCS Log tag and cq_abort
!!   15:48, 08/04/2003 drb 
!!    Bug fixes in splint and dsplint
!!***
module spline_module
  
  use global_module, only: area_general

  implicit none

  ! RCS tag for object file identification 
  character(len=80), save, private :: RCSid = "$Id$"

contains

! -----------------------------------------------------------
! Subroutine spline
! -----------------------------------------------------------

!!****f* spline_module/spline *
!!
!!  NAME 
!!   spline
!!  USAGE
!! 
!!  PURPOSE
!!   Given an array y(1:n) containing a function tabulated at 
!!   REGULAR intervals with spacing delta_x, and given values
!!   yp1 and ypn for the first derivative of the function at 
!!   the end points, this subroutine returns an array d2y(1:n) 
!!   containing the second derivatives of the interpolating 
!!   function at the tabulated points. If yp1 and/or ypn are 
!!   equal or larger than 1.0d30, the routine is signaled to set 
!!   the corresponding boundary condition for a natural spline, 
!!   with zero second derivative at that point. This subroutine
!!   is based on subroutine spline(), from Numerical Recipes (2nd edition). 
!!  INPUTS
!!   integer :: n - size of table to be splined
!!   real(double) :: delta_x - spacing of values in table
!!   real(double) :: yp1, ypn - end points giving boundaries
!!   real(double), dimension(n) :: y, y2 - table and second derivative
!!  USES
!!   datatypes, numbers
!!  AUTHOR
!!   E.H.Hernandez
!!  CREATION DATE
!!   29/02/96
!!  MODIFICATION HISTORY
!!   23/05/2001 dave
!!    F90, ROBODoc, indentation
!!   31/05/2001 dave
!!    More ROBODoc comments
!!  SOURCE
!!
  subroutine spline( n, delta_x, y, yp1, ypn, y2 )

    use datatypes
    use numbers
    use GenComms, only: cq_abort
    use memory_module, only: reg_alloc_mem, reg_dealloc_mem, type_dbl

    implicit none

    ! Passed variables
    integer :: n

    real(double) :: delta_x, yp1, ypn
    real(double), dimension(n) :: y, y2

    ! Local variables
    integer :: i, k, stat
    real(double) :: p, qn, sig, un
    real(double), parameter ::huge = 1.0e30_double
    real(double), dimension(:), allocatable :: u, x

    allocate(u(n), x(n), STAT=stat)
    if (stat /= 0) call cq_abort("spline: Error alloc mem: ", n)
    call reg_alloc_mem(area_general, 2*n, type_dbl)

    ! generate the array of x entries
    do i=1, n
       x(i) = delta_x * real(i-1,double)
    end do
    if ( yp1 .gt. huge ) then ! lower boundary condition set to 'natural'
       y2(1) = zero
       u(1) = zero
    else ! or else to have a specified first derivative
       y2(1) = -half
       u(1) = ( three/(x(2) - x(1))) *  &
            ((y(2) - y(1))/(x(2) - x(1)) - yp1 )
    end if

    ! this is the decomposition loop for the tridiagonal algorithm; y2 and
    ! u are used for temporary storage of the decomposed factors
    do i=2, n-1
       sig = (x(i) - x(i-1)) / (x(i+1) - x(i-1))
       p = sig * y2(i-1) + two
       y2(i) = (sig - one) / p
       u(i) = (six * ((y(i+1) - y(i))/(x(i+1) - x(i)) -  &
            (y(i) - y(i-1))/ &
            (x(i) - x(i-1)))/(x(i+1) - x(i-1)) - sig * u(i-1))/p 
    end do

    if (ypn>huge) then ! the upper boundary condition is set to be 'natural'...
       qn = zero
       un = zero
    else ! or else to have a specified first derivative
       qn = half
       un = ( three/(x(n) - x(n-1))) *  &
            (ypn - (y(n) - y(n-1))/(x(n) - x(n-1)))
    end if
    y2(n) = (un - qn * u(n-1)) / (qn * y2(n-1) + one )

    ! this is the backsubstitution loop of the tridiagonal algorithm
    do k=n-1, 1, -1
       y2(k) = y2(k) * y2(k+1) + u(k)
    end do

    deallocate(u, x, STAT=stat)
    if (stat /= 0) call cq_abort("spline: Error dealloc mem")
    call reg_dealloc_mem(area_general, 2*n, type_dbl)

    return
  end subroutine spline
!!***

! -----------------------------------------------------------
! Subroutine spline_nonU
! -----------------------------------------------------------

!!****f* spline_module/spline_nonU *
!!
!!  NAME 
!!   spline_nonU - makes a spline table for non-uniform x array
!!  USAGE
!!   spline_nonU(number, x array, y array,start dy/dx,final dy/dx, d2y/dx2)
!!   spline_nonU(n, y, y, yp1, ypn, y2)
!!  PURPOSE
!!   Given an array y(1:n) containing a function tabulated at 
!!   IRREGULAR intervals given by array x(1:n), and given values
!!   yp1 and ypn for the first derivative of the function at 
!!   the end points, this subroutine returns an array d2y(1:n) 
!!   containing the second derivatives of the interpolating 
!!   function at the tabulated points. If yp1 and/or ypn are 
!!   equal or larger than 1.0d30, the routine is signaled to set 
!!   the corresponding boundary condition for a natural spline, 
!!   with zero second derivative at that point. This subroutine
!!   is based on subroutine spline(), from Numerical Recipes (2nd edition). 
!!  INPUTS
!! 
!! 
!!  USES
!! 
!!  AUTHOR
!!   D.R.Bowler
!!  CREATION DATE
!!   20/11/97
!!  MODIFICATION HISTORY
!!   23/05/2001 dave
!!    F90, ROBODoc, indentation
!!  SOURCE
!!
  subroutine spline_nonU(n, x, y, yp1, ypn, y2)

    use datatypes
    use numbers
    use GenComms, only: cq_abort
    use memory_module, only: reg_alloc_mem, reg_dealloc_mem, type_dbl

    implicit none

    ! Passed variables

    integer :: n
    real(double) :: yp1, ypn
    real(double), dimension(n) :: x, y, y2

    !     Local variables
    integer i, k, stat
    real(double) :: p, qn, sig, un
    real(double), parameter :: huge = 1.0e30_double
    real(double), dimension(:), allocatable :: u

    allocate(u(n), STAT=stat)
    if (stat /= 0) call cq_abort("spline_nonU: Error alloc mem: ", n)
    call reg_alloc_mem(area_general, n, type_dbl)

    ! set up boundary conditions
    if (yp1>huge) then ! lower boundary condition set to 'natural'
       y2(1) = zero
       u(1) = zero
    else ! or else to have a specified first derivative
       y2(1) = -half
       u(1) = ( three/(x(2) - x(1))) * ((y(2) - y(1))/(x(2) - x(1)) - yp1 )
    end if

    ! this is the decomposition loop for the tridiagonal algorithm; y2 and
    ! u are used for temporary storage of the decomposed factors
    do i=2, n-1
       sig = (x(i) - x(i-1)) / (x(i+1) - x(i-1))
       p = sig * y2(i-1) + two
       y2(i) = (sig - one) / p
       u(i) = (six * ((y(i+1) - y(i))/(x(i+1) - x(i)) -  &
            (y(i) - y(i-1))/ &
            (x(i) - x(i-1)))/(x(i+1) - x(i-1)) - sig * u(i-1))/p 
    end do

    if (ypn>huge) then ! the upper boundary condition is set to be 'natural'...
       qn = zero
       un = zero
    else ! or else to have a specified first derivative
       qn = half
       un = ( three/(x(n) - x(n-1))) *  &
            (ypn - (y(n) - y(n-1))/(x(n) - x(n-1)))
    end if
    y2(n) = (un - qn * u(n-1)) / (qn * y2(n-1) + one )

    ! this is the backsubstitution loop of the tridiagonal algorithm
    do k=n-1, 1, -1
       y2(k) = y2(k) * y2(k+1) + u(k)
    end do

    deallocate(u, STAT=stat)
    if (stat /= 0) call cq_abort("spline_nonU: Error dealloc mem")
    call reg_dealloc_mem(area_general, n, type_dbl)

    return
  end subroutine spline_nonU
!!***

! -----------------------------------------------------------
! Subroutine splint_nonU
! -----------------------------------------------------------

!!****f* spline_module/splint_nonU *
!!
!!  NAME 
!!   splint_nonU
!!  USAGE
!!   splint_nonU(xarray, yarray, d2y/dx2 array, size of array, x value, y, flag)
!!   splint_nonU(xa,ya,y2a,n,x,y, flag)
!!  PURPOSE
!!   NR routine for interpolation using spline table set up by spline
!!   Takes a spline table with non-uniform x spacings and, for the 
!!   specified value x, finds the y value.
!!  INPUTS
!!   integer :: n (size of table)
!!   real(double) :: x,y (x value specified, y value to be found)
!!   real(double), dimension(n) :: xa, ya, y2a (x, y, d2y/dx2 arrays)
!!   logical :: flag (flag to warn if given x is out of tabulated range)
!!  USES
!!   datatypes
!!  AUTHOR
!!   D.R.Bowler
!!  CREATION DATE
!!   20/11/97
!!  MODIFICATION HISTORY
!!   23/05/2001 dave
!!    F90, ROBODoc, indentation
!!   20/06/2001 dave
!!    Added cq_abort
!!   15/04/02 lkd
!!    Changed title from splint to splint_nonU for consistency with rest of module
!!   15/04/02 lkd
!!    Added "use numbers" statement
!!   17/04/2002 lkd
!!    Added logical flag to return false if all OK, true if the given x is out of range
!! 
!!  SOURCE
!!
  subroutine splint_nonU(xa,ya,y2a,n,x,y,flag)

    use datatypes
    use numbers
    use GenComms, ONLY: cq_abort

    ! Passed variables
    integer, intent(in) :: n
    real(double), intent(in) :: x
    real(double), intent(out) :: y
    real(double), dimension(n), intent(in) :: xa,y2a,ya
    logical, intent(out) :: flag

    ! Local variables
    integer :: k,khi,klo
    real(double) :: a,b,h

    ! initialize warning flag
    flag = .false.
    if (x .gt. xa(n) .or. x .lt. xa(1)) then
       flag = .true. ! the given x lies outside the tabulated values, warn the calling routine
    end if

    ! Search through the table and bracket x value
    klo = 1
    khi = n
    do while((khi-klo)>1)
       k=(khi+klo)/2
       if(xa(k).gt.x) then ! Increment high point
          khi=k
       else                ! Increment low point
          klo=k
       end if
    end do
    h=xa(khi)-xa(klo)      ! Find x spacing
    if(h==0) then
       call cq_abort('splint_nonU: bisection failed, h=0')
    end if
    ! Create value of y
    a = (xa(khi)-x)/h
    b = (x-xa(klo))/h
    y=a*ya(klo)+b*ya(khi)+ &
         ( (a*a*a-a)*y2a(klo)+(b*b*b-b)*y2a(khi) )*(h*h)/six
    return
  end subroutine splint_nonU
!!***

! -----------------------------------------------------------
! Subroutine splint
! -----------------------------------------------------------

!!****f* spline_module/splint *
!!
!!  NAME 
!!   splint
!!  USAGE
!!   splint(xarray, yarray, d2y/dx2 array, size of array, x value, y)
!!   splint(xa,ya,y2a,n,x,y)
!!  PURPOSE
!!   NR routine for interpolation using spline table set up by spline
!!   Takes a spline table with uniform x spacing delta_x and, for the 
!!   specified value x, finds the y value.
!!  INPUTS
!!   integer :: n (size of table)
!!   real(double) :: x, y, delta_x (x value specified, y value to be found, x-value spacing)
!!   real(double), dimension(n) :: ya, y2a (y, d2y/dx2 arrays)
!!   logical :: flag (flag to warn if given x is out of tabulated range)
!!  USES
!!   datatypes, numbers, GenComms
!!  AUTHOR
!!   L.K.Dash
!!  CREATION DATE
!!   15/04/02
!!  MODIFICATION HISTORY
!!   15:47, 08/04/2003 drb 
!!    Fixed major bug: array references were all one element off
!!   10:00, 28/08/2003 drb 
!!    Fixed ANOTHER major bug: I'd forgotten that the array reference was used to create the x position so needs one taken off it !
!!   08:56, 2003/11/20 dave
!!    Added return to prevent array overrun
!!  SOURCE
!!
  subroutine splint(delta_x,ya,y2a,n,x,y,flag)

    use datatypes
    use numbers
    use GenComms, ONLY: cq_abort

    ! Passed variables
    integer, intent(in) :: n
    real(double), intent(in) :: x
    real(double), intent(out) :: y
    real(double), intent(in) :: delta_x
    real(double), dimension(n), intent(in) :: y2a,ya
    logical, intent(out) :: flag

    ! Local variables
    integer :: khi,klo
    real(double) :: a,b,xlo,xhi

    ! initialize warning flag
    flag = .false.
    if (x>real(n-1,double)*delta_x .or. x<0) then
       flag = .true. ! the given x lies outside the tabulated values, warn the calling routine
       y = 0.0_double ! Avoid array overruns
       return
    end if

    ! bracket x value   
    klo = aint(x/delta_x) + 1 ! integer division gives lower bound, add 1 for array value
    khi = klo + 1
    if(khi>n) call cq_abort("Error in splint: index beyond range ",khi,n)
    ! Remember that we added to get array value and subtract it off again !
    xlo =delta_x*real(klo-1,double)
    xhi =delta_x*real(khi-1,double)
    
    ! Create value of y
    a = (xhi-x)/delta_x
    b = (x-xlo)/delta_x
    y=a*ya(klo)+b*ya(khi)+ &
         ( (a*a*a-a)*y2a(klo)+(b*b*b-b)*y2a(khi) )*(delta_x*delta_x)/six
    return
  end subroutine splint
!!***

! -----------------------------------------------------------
! Subroutine dsplint
! -----------------------------------------------------------

!!****f* spline_module/dsplint *
!!
!!  NAME 
!!   dsplint
!!  USAGE
!!   dsplint(xarray, yarray, d2y/dx2 array, size of array, x value, dy)
!!   dsplint(xa,ya,y2a,n,x,dy)
!!  PURPOSE
!!   NR routine for interpolation using spline table set up by spline
!!   Takes a spline table with uniform x spacing delta_x and, for the 
!!   specified value x, finds the derivative value.
!!  INPUTS
!!   integer :: n (size of table)
!!   real(double) :: x, y, delta_x (x value specified, y value to be found, x-value spacing)
!!   real(double), dimension(n) :: ya, y2a (y, d2y/dx2 arrays)
!!   logical :: flag (flag to warn if given x is out of tabulated range)
!!  USES
!!   datatypes, numbers, GenComms
!!  AUTHOR
!!   D.R.Bowler
!!  CREATION DATE
!!   07:50, 2003/03/19 dave
!!  MODIFICATION HISTORY
!!   15:47, 08/04/2003 drb 
!!    Fixed major bug: array references were all one element off (carried over from splint)
!!   11:04, 03/09/2003 drb 
!!    Added xlo/xhi correction that was added to splint (should have been done earlier !)
!!   08:56, 2003/11/20 dave
!!    Added return to prevent array overrun
!!   08:01, 2003/12/04 dave
!!    Added return of spline-interpolated value as well as derivative
!!  SOURCE
!!
  subroutine dsplint(delta_x,ya,y2a,n,x,y,dy,flag)

    use datatypes
    use numbers
    use GenComms, ONLY: cq_abort

    ! Passed variables
    integer, intent(in) :: n
    real(double), intent(in) :: x
    real(double), intent(out) :: y,dy
    real(double), intent(in) :: delta_x
    real(double), dimension(n), intent(in) :: y2a,ya
    logical, intent(out) :: flag

    ! Local variables
    integer :: khi,klo
    real(double) :: a,b,da,db,dc,dd, xlo, xhi

    ! initialize warning flag
    flag = .false.
    if (x>(n-1)*delta_x.OR.x<0) then
       flag = .true. ! the given x lies outside the tabulated values, warn the calling routine
       y = 0.0_double
       dy = 0.0_double ! Avoid array overruns
       return
    end if

    ! bracket x value   
    klo = aint(x/delta_x) + 1 ! integer division gives lower bound, add 1 for array value
    khi = klo + 1
    ! Remember that we added to get array value and subtract it off again !
    xlo =delta_x*real(klo-1,double)
    xhi =delta_x*real(khi-1,double)

    ! Create value of y
    a = (xhi-x)/delta_x
    b = (x-xlo)/delta_x
    da = -one/delta_x
    db =  one/delta_x
    dc = -delta_x*(three*a*a-one)/six
    dd =  delta_x*(three*b*b-one)/six
    y=a*ya(klo)+b*ya(khi)+ &
         ( (a*a*a-a)*y2a(klo)+(b*b*b-b)*y2a(khi) )*(delta_x*delta_x)/six
    dy =  da*ya(klo)+db*ya(khi)+dc*y2a(klo)+dd*y2a(khi)
    return
  end subroutine dsplint
!!***

end module spline_module



