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