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.f90"
  include "dqk15_2d_inner.f90"
  include "dqk31_2d_outer.f90"
  include "dqk31_2d_inner.f90"
  include "dqage.f90"
  include "dqk15.f90"
  include "dqk21.f90"
  include "dqk31.f90"
  include "dqk41.f90"
  include "dqk51.f90"
  include "dqk61.f90"
  include "dqk41_2d_outer.f90"
  include "dqk41_2d_inner.f90"
  include "dqag_2d_outer.f90"
  include "dqpsrt.f90"
  include "dqag.f90"
  include "dqage_2d_outer.f90"
  include "dqage_2d_inner.f90"
  include "dqk51_2d_outer.f90"
  include "dqk51_2d_inner.f90"
  include "dqk61_2d_outer.f90"
  include "dqk21_2d_outer.f90"
  include "dqk61_2d_inner.f90"
  include "dqk21_2d_inner.f90"
  include "dqk15_2d_outer.f90"

  
  
  end module fvn_integ