module fvn_integ
use fvn_common
implicit none

! Gauss legendre
interface fvn_gauss_legendre
    module procedure fvn_d_gauss_legendre
end interface fvn_gauss_legendre

! Simple Gauss Legendre integration
interface fvn_gl_integ
    module procedure fvn_d_gl_integ
end interface fvn_gl_integ

! Adaptative Gauss Kronrod integration f(x)
interface fvn_integ_1_gk
    module procedure fvn_d_integ_1_gk
end interface fvn_integ_1_gk

! Adaptative Gauss Kronrod integration f(x,y)
interface fvn_integ_2_gk
    module procedure fvn_d_integ_2_gk
end interface fvn_integ_2_gk


contains
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!   Integration
!
!   Only double precision coded atm
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


subroutine fvn_d_gauss_legendre(n,qx,qw)
!
! This routine compute the n Gauss Legendre abscissas and weights
! Adapted from Numerical Recipes routine gauleg
!
! n (in) : number of points
! qx(out) : abscissas
! qw(out) : weights
!
implicit none
double precision,parameter :: pi=3.141592653589793d0
integer, intent(in) :: n
double precision, intent(out) :: qx(n),qw(n)

integer :: m,i,j
double precision :: z,z1,p1,p2,p3,pp
m=(n+1)/2
do i=1,m
    z=cos(pi*(dble(i)-0.25d0)/(dble(n)+0.5d0))
iloop:  do 
            p1=1.d0
            p2=0.d0
            do j=1,n
                p3=p2
                p2=p1
                p1=((2.d0*dble(j)-1.d0)*z*p2-(dble(j)-1.d0)*p3)/dble(j)
            end do
            pp=dble(n)*(z*p1-p2)/(z*z-1.d0)
            z1=z
            z=z1-p1/pp
            if (dabs(z-z1)<=epsilon(z)) then
                exit iloop
            end if
        end do iloop
    qx(i)=-z
    qx(n+1-i)=z
    qw(i)=2.d0/((1.d0-z*z)*pp*pp)
    qw(n+1-i)=qw(i)
end do
end subroutine



subroutine fvn_d_gl_integ(f,a,b,n,res)
!
! This is a simple non adaptative integration routine 
! using n gauss legendre abscissas and weights
!
!   f(in)   : the function to integrate
!   a(in)   : lower bound
!   b(in)   : higher bound
!   n(in)   : number of gauss legendre pairs
!   res(out): the evaluation of the integral
!
double precision,external :: f
double precision, intent(in) :: a,b
integer, intent(in):: n
double precision, intent(out) :: res

double precision, allocatable :: qx(:),qw(:)
double precision :: xm,xr
integer :: i

! First compute n gauss legendre abs and weight
allocate(qx(n))
allocate(qw(n))
call fvn_d_gauss_legendre(n,qx,qw)

xm=0.5d0*(b+a)
xr=0.5d0*(b-a)

res=0.d0

do i=1,n
    res=res+qw(i)*f(xm+xr*qx(i))
end do

res=xr*res

deallocate(qw)
deallocate(qx)

end subroutine

!!!!!!!!!!!!!!!!!!!!!!!!
!
! Simple and double adaptative Gauss Kronrod integration based on
! a modified version of quadpack ( http://www.netlib.org/quadpack
!
! Common parameters :
!
!       key (in)
!       epsabs
!       epsrel
!
!
!!!!!!!!!!!!!!!!!!!!!!!!

subroutine fvn_d_integ_1_gk(f,a,b,epsabs,epsrel,key,res,abserr,ier,limit)
!
! Evaluate the integral of function f(x) between a and b
!
! f(in) : the function
! a(in) : lower bound
! b(in) : higher bound
! epsabs(in) : desired absolute error
! epsrel(in) : desired relative error
! key(in) : gauss kronrod rule
!                     1:   7 - 15 points
!                     2:  10 - 21 points
!                     3:  15 - 31 points
!                     4:  20 - 41 points
!                     5:  25 - 51 points
!                     6:  30 - 61 points
!
! limit(in) : maximum number of subintervals in the partition of the 
!               given integration interval (a,b). A value of 500 will give the same
!               behaviour as the imsl routine dqdag
!
! res(out) : estimated integral value
! abserr(out) : estimated absolute error
! ier(out) : error flag from quadpack routines
!               0 : no error
!               1 : maximum number of subdivisions allowed
!                   has been achieved. one can allow more
!                   subdivisions by increasing the value of
!                   limit (and taking the according dimension
!                   adjustments into account). however, if
!                   this yield no improvement it is advised
!                   to analyze the integrand in order to
!                   determine the integration difficulaties.
!                   if the position of a local difficulty can
!                   be determined (i.e.singularity,
!                   discontinuity within the interval) one
!                   will probably gain from splitting up the
!                   interval at this point and calling the
!                   integrator on the subranges. if possible,
!                   an appropriate special-purpose integrator
!                   should be used which is designed for
!                   handling the type of difficulty involved.
!               2 : the occurrence of roundoff error is
!                   detected, which prevents the requested
!                   tolerance from being achieved.
!               3 : extremely bad integrand behaviour occurs
!                   at some points of the integration
!                   interval.
!               6 : the input is invalid, because
!                   (epsabs.le.0 and
!                   epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
!                   or limit.lt.1 or lenw.lt.limit*4.
!                   result, abserr, neval, last are set
!                   to zero.
!                   except when lenw is invalid, iwork(1),
!                   work(limit*2+1) and work(limit*3+1) are
!                   set to zero, work(1) is set to a and
!                   work(limit+1) to b.

implicit none
double precision, external :: f
double precision, intent(in) :: a,b,epsabs,epsrel
integer, intent(in) :: key
integer, intent(in),optional :: limit
double precision, intent(out) :: res,abserr
integer, intent(out) :: ier

double precision, allocatable :: work(:)
integer, allocatable :: iwork(:)
integer :: lenw,neval,last
integer :: limitw

! imsl value for limit is 500
limitw=500
if (present(limit)) limitw=limit

lenw=limitw*4

allocate(iwork(limitw))
allocate(work(lenw))

call dqag(f,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limitw,lenw,last,iwork,work)

deallocate(work)
deallocate(iwork)

end subroutine



subroutine fvn_d_integ_2_gk(f,a,b,g,h,epsabs,epsrel,key,res,abserr,ier,limit)
!
! Evaluate the double integral of function f(x,y) for x between a and b and y between g(x) and h(x)
!
! f(in) : the function
! a(in) : lower bound
! b(in) : higher bound
! g(in) : external function describing lower bound for y
! h(in) : external function describing higher bound for y
! epsabs(in) : desired absolute error
! epsrel(in) : desired relative error
! key(in) : gauss kronrod rule
!                     1:   7 - 15 points
!                     2:  10 - 21 points
!                     3:  15 - 31 points
!                     4:  20 - 41 points
!                     5:  25 - 51 points
!                     6:  30 - 61 points
!
! limit(in) : maximum number of subintervals in the partition of the 
!               given integration interval (a,b). A value of 500 will give the same
!               behaviour as the imsl routine dqdag
!
! res(out) : estimated integral value
! abserr(out) : estimated absolute error
! ier(out) : error flag from quadpack routines
!               0 : no error
!               1 : maximum number of subdivisions allowed
!                   has been achieved. one can allow more
!                   subdivisions by increasing the value of
!                   limit (and taking the according dimension
!                   adjustments into account). however, if
!                   this yield no improvement it is advised
!                   to analyze the integrand in order to
!                   determine the integration difficulaties.
!                   if the position of a local difficulty can
!                   be determined (i.e.singularity,
!                   discontinuity within the interval) one
!                   will probably gain from splitting up the
!                   interval at this point and calling the
!                   integrator on the subranges. if possible,
!                   an appropriate special-purpose integrator
!                   should be used which is designed for
!                   handling the type of difficulty involved.
!               2 : the occurrence of roundoff error is
!                   detected, which prevents the requested
!                   tolerance from being achieved.
!               3 : extremely bad integrand behaviour occurs
!                   at some points of the integration
!                   interval.
!               6 : the input is invalid, because
!                   (epsabs.le.0 and
!                   epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
!                   or limit.lt.1 or lenw.lt.limit*4.
!                   result, abserr, neval, last are set
!                   to zero.
!                   except when lenw is invalid, iwork(1),
!                   work(limit*2+1) and work(limit*3+1) are
!                   set to zero, work(1) is set to a and
!                   work(limit+1) to b.

implicit none
double precision, external:: f,g,h
double precision, intent(in) :: a,b,epsabs,epsrel
integer, intent(in) :: key
integer, intent(in), optional :: limit
integer, intent(out) :: ier
double precision, intent(out) :: res,abserr


double precision, allocatable :: work(:)
integer :: limitw
integer, allocatable :: iwork(:)
integer :: lenw,neval,last

! imsl value for limit is 500
limitw=500
if (present(limit)) limitw=limit

lenw=limitw*4
allocate(work(lenw))
allocate(iwork(limitw))

call dqag_2d_outer(f,a,b,g,h,epsabs,epsrel,key,res,abserr,neval,ier,limitw,lenw,last,iwork,work)

deallocate(iwork)
deallocate(work)
end subroutine



subroutine fvn_d_integ_2_inner_gk(f,x,a,b,epsabs,epsrel,key,res,abserr,ier,limit)
!
! Evaluate the single integral of function f(x,y) for y between a and b with a
! given x value
!
! This function is used for the evaluation of the double integral fvn_d_integ_2_gk
!
! f(in) : the function
! x(in) : x
! a(in) : lower bound
! b(in) : higher bound
! epsabs(in) : desired absolute error
! epsrel(in) : desired relative error
! key(in) : gauss kronrod rule
!                     1:   7 - 15 points
!                     2:  10 - 21 points
!                     3:  15 - 31 points
!                     4:  20 - 41 points
!                     5:  25 - 51 points
!                     6:  30 - 61 points
!
! limit(in) : maximum number of subintervals in the partition of the 
!               given integration interval (a,b). A value of 500 will give the same
!               behaviour as the imsl routine dqdag
!
! res(out) : estimated integral value
! abserr(out) : estimated absolute error
! ier(out) : error flag from quadpack routines
!               0 : no error
!               1 : maximum number of subdivisions allowed
!                   has been achieved. one can allow more
!                   subdivisions by increasing the value of
!                   limit (and taking the according dimension
!                   adjustments into account). however, if
!                   this yield no improvement it is advised
!                   to analyze the integrand in order to
!                   determine the integration difficulaties.
!                   if the position of a local difficulty can
!                   be determined (i.e.singularity,
!                   discontinuity within the interval) one
!                   will probably gain from splitting up the
!                   interval at this point and calling the
!                   integrator on the subranges. if possible,
!                   an appropriate special-purpose integrator
!                   should be used which is designed for
!                   handling the type of difficulty involved.
!               2 : the occurrence of roundoff error is
!                   detected, which prevents the requested
!                   tolerance from being achieved.
!               3 : extremely bad integrand behaviour occurs
!                   at some points of the integration
!                   interval.
!               6 : the input is invalid, because
!                   (epsabs.le.0 and
!                   epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
!                   or limit.lt.1 or lenw.lt.limit*4.
!                   result, abserr, neval, last are set
!                   to zero.
!                   except when lenw is invalid, iwork(1),
!                   work(limit*2+1) and work(limit*3+1) are
!                   set to zero, work(1) is set to a and
!                   work(limit+1) to b.

implicit none
double precision, external:: f
double precision, intent(in) :: x,a,b,epsabs,epsrel
integer, intent(in) :: key
integer, intent(in),optional :: limit
integer, intent(out) :: ier
double precision, intent(out) :: res,abserr


double precision, allocatable :: work(:)
integer :: limitw
integer, allocatable :: iwork(:)
integer :: lenw,neval,last

! imsl value for limit is 500
limitw=500
if (present(limit)) limitw=limit

lenw=limitw*4
allocate(work(lenw))
allocate(iwork(limitw))

call dqag_2d_inner(f,x,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limitw,lenw,last,iwork,work)

deallocate(iwork)
deallocate(work)
end subroutine


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Include the modified quadpack files
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
include "dqag_2d_inner.f"
include "dqk15_2d_inner.f"
include "dqk31_2d_outer.f"
include "dqk31_2d_inner.f"
include "dqage.f"
include "dqk15.f"
include "dqk21.f"
include "dqk31.f"
include "dqk41.f"
include "dqk51.f"
include "dqk61.f"
include "dqk41_2d_outer.f"
include "dqk41_2d_inner.f"
include "dqag_2d_outer.f"
include "dqpsrt.f"
include "dqag.f"
include "dqage_2d_outer.f"
include "dqage_2d_inner.f"
include "dqk51_2d_outer.f"
include "dqk51_2d_inner.f"
include "dqk61_2d_outer.f"
include "dqk21_2d_outer.f"
include "dqk61_2d_inner.f"
include "dqk21_2d_inner.f"
include "dqk15_2d_outer.f"


end module fvn_integ