subroutine quad_sample(fsample,a,b,&
epsabs,epsrel,&
key,&
singular_endpoint,&
singular_points,&
cpole,&
alfa,beta,&
omega,&
weight_func,&
verbose,&
Ninterp,&
result)
real(8),dimension(:) :: fsample
real(8) :: a
real(8) :: b
!optional variables
real(8),optional :: epsabs
real(8),optional :: epsrel
integer,optional :: key !order of GK rule in QAG
logical,optional :: singular_endpoint !T if singular_endpoint exists (QAGS)
real(8),dimension(:),optional :: singular_points !location of singular points in QAGP
real(8),optional :: cpole !location of the pole in QAWC f(x)/x-cpole
real(8),optional :: alfa,beta !parameters for QAWS
real(8),optional :: omega !frequency of the weight funcion in QAWO
integer,optional :: weight_func !if(QAWF)then
! ! weight_func=1,2 -> cos(omega*x),sin(omega*x)
! !if(QAWS)then
! ! weight_func = 1 (x-a)**alfa*(b-x)**beta
! ! weight_func = 2 (x-a)**alfa*(b-x)**beta*log(x-a)
! ! weight_func = 3 (x-a)**alfa*(b-x)**beta*log(b-x)
! ! weight_func = 4 (x-a)**alfa*(b-x)**beta*log(x-a)*log(b-x)
logical,optional :: verbose
integer,optional :: Ninterp
real(8) :: result
!actual default variables
real(8) :: epsabs_
real(8) :: epsrel_
logical :: verbose_
real(8) :: abserr
real(8),dimension(size(fsample)) :: xsample
integer :: Neval
integer :: Ier,i
integer :: Nsingular_points
integer :: Ninterp_,Lin
character(len=20) :: routine_name
!
type finter_type
real(8),allocatable :: X(:)
real(8),allocatable :: F(:)
integer :: Imin=0,Imax=0,N=0
logical :: status=.false.
end type finter_type
type(finter_type) :: Finterp
!
!
epsabs_ = 1d-12 ; if(present(epsabs))epsabs_=epsabs
epsrel_ = 1d-6 ; if(present(epsrel))epsrel_=epsrel
verbose_=.false.; if(present(verbose))verbose_=verbose
Ninterp_=3 ; if(present(Ninterp))Ninterp_=Ninterp
!
!Build the interpolating function to be integrated:
Lin = size(fsample)
xsample = linspace(a,b,Lin)
call init_finter(Finterp,xsample,fsample,Ninterp_)
!
if(present(key).AND.(present(alfa).OR.present(beta)))stop "ERROR in quad: key & alfa,beta"
if(present(key).AND.present(cpole))stop "ERROR in quad: key & cpole"
if(present(key).AND.present(singular_points))stop "ERROR in quad: key & singular_points"
if(present(key).AND.present(singular_endpoint))stop "ERROR in quad: key & singular_endpoint"
if(present(key).AND.present(omega))stop "ERROR in quad: key & omega"
if(present(singular_endpoint).AND.(present(alfa).OR.present(beta)))stop "ERROR in quad: singular_endpoint & alfa,beta"
if(present(singular_endpoint).AND.present(cpole))stop "ERROR in quad: singular_endpoint & cpole"
if(present(singular_endpoint).AND.present(singular_points))stop "ERROR in quad: singular_endpoint & singular_points"
if(present(singular_endpoint).AND.present(omega))stop "ERROR in quad: singular_endpoint & omega"
if(present(cpole).AND.(present(alfa).OR.present(beta)))stop "ERROR in quad: cpole & alfa,beta"
if(present(cpole).AND.present(singular_points))stop "ERROR in quad: cpole & singular_points"
if(present(cpole).AND.present(omega))stop "ERROR in quad: cpole & omega"
if(present(omega).AND.(present(alfa).OR.present(beta)))stop "ERROR in quad: omega & alfa,beta"
if(present(omega).AND.present(singular_points))stop "ERROR in quad: omega & singular_points"
if(present(singular_points).AND.(present(alfa).OR.present(beta)))stop "ERROR in quad: singular_points & alfa,beta"
!
!QNG: no need
!QAGS: need singular_endpoint=T
!QAG : need key=1,6
!QAGP: need singular_points
!QAWC: need cpole
!QAWO: need omega + weight_func=1,2
!QAWS: need alfa + beta + weight_func=1,2,3,4
routine_name='QNG'
if(present(singular_endpoint))routine_name='QAGS'
if(present(key))routine_name='QAG'
if(present(singular_points))routine_name='QAGP'
if(present(cpole))routine_name='QAWC'
if(present(omega))routine_name='QAWO'
if(present(alfa).AND.present(beta))routine_name='QAWS'
!
if(verbose_)write(*,"(A,A)")"QUAD: selected procedure =", routine_name
!
select case(routine_name)
case ('QNG')
call QNG(func,a,b,epsabs_,epsrel_,result,abserr,neval,ier)
if(verbose_)then
write(*,'(A,2F14.6)')'Endpoints (a,b) =', a,b
write(*,'(A,F14.6)') 'Estimated integral =', result
write(*,'(A,F14.6)') 'Estimated integral error =', abserr
write(*,'(A,I8)') 'Error return code IER =', ier
else
if(ier/=0)then
write(*,'(A,I8)') 'Error return code IER =', ier
stop
endif
endif
!
!
case ('QAGS')
call QAGS(func,a,b,epsabs_,epsrel_,result,abserr,neval,ier)
if(verbose_)then
write(*,'(A,2F14.6)')'(Singular) endpoints (a*,b*) =', a,b
write(*,'(A,F14.6)') 'Estimated integral =', result
write(*,'(A,F14.6)') 'Estimated integral error =', abserr
write(*,'(A,I8)') 'Error return code IER =', ier
else
if(ier/=0)then
write(*,'(A,I8)') 'Error return code IER =', ier
stop
endif
endif
!
!
case ('QAG')
if(key<1.OR.key>6)stop "ERROR in quad: use QAG, key !=[1:6]"
call QAG(func,a,b,epsabs_,epsrel_,key,result,abserr,neval,ier)
if(verbose_)then
write(*,'(A,2F14.6)')'Endpoints (a,b) =', a,b
write(*,'(A,I8)') 'Approximation or =', key
write(*,'(A,F14.6)') 'Estimated integral =', result
write(*,'(A,F14.6)') 'Estimated integral error =', abserr
write(*,'(A,I8)') 'Error return code IER =', ier
else
if(ier/=0)then
write(*,'(A,I8)') 'Error return code IER =', ier
stop
endif
endif
!
!
case ('QAGP')
Nsingular_points=size(singular_points)
call QAGP(func,a,b,Nsingular_points,singular_points,epsabs_,epsrel_,result,abserr,neval,ier)
if(verbose_)then
write(*,'(A,2F14.6)')'Endpoints (a,b) =', a,b
write(*,'(A,I8)')'N_singular_points =', Nsingular_points
do i=1,Nsingular_points
write(*,"(A,I6,F14.6)")'I_singular_point =',i,singular_points(i)
enddo
write(*,'(A,F14.6)') 'Estimated integral =', result
write(*,'(A,F14.6)') 'Estimated integral error =', abserr
write(*,'(A,I8)') 'Error return code IER =', ier
else
if(ier/=0)then
write(*,'(A,I8)') 'Error return code IER =', ier
stop
endif
endif
!
!
case ('QAWC')
call QAWC(func,a,b,cpole,epsabs_,epsrel_,result,abserr,neval,ier)
if(verbose_)then
write(*,'(A,2F14.6)')'Endpoints (a,b) =', a,b
write(*,'(A,F14.6)')'C_pole =', cpole
write(*,'(A,F14.6)') 'Estimated integral =', result
write(*,'(A,F14.6)') 'Estimated integral error =', abserr
write(*,'(A,I8)') 'Error return code IER =', ier
else
if(ier/=0)then
write(*,'(A,I8)') 'Error return code IER =', ier
stop
endif
endif
!
!
case ('QAWO')
if(weight_func<1.OR.weight_func>2)stop "ERROR in quad: use QAWO, weight_func !=[1,2]"
if(verbose_)then
if(weight_func==1)then
write(*,'(A)')'W(x) function is cos(omega*x)'
elseif(weight_func==2)then
write(*,'(A)')'W(x) function is sin(omega*x)'
endif
endif
call QAWO(func,a,b,omega,weight_func,epsabs_,epsrel_,result,abserr,neval,ier)
if(verbose_)then
write(*,'(A,2F14.6)')'Endpoints (a,b) =', a,b
write(*,'(A,F14.6)') 'Omega =', omega
write(*,'(A,F14.6)') 'Estimated integral =', result
write(*,'(A,F14.6)') 'Estimated integral error =', abserr
write(*,'(A,I8)') 'Error return code IER =', ier
else
if(ier/=0)then
write(*,'(A,I8)') 'Error return code IER =', ier
stop
endif
endif
!
!
case ('QAWS')
if(weight_func<1.OR.weight_func>4)stop "ERROR in quad: use QAWS, weight_func !=[1,..,4]"
if(verbose_)then
if(weight_func==1)then
write(*,'(A)')'W(x) function is (x-a)**alfa*(b-x)**beta '
elseif(weight_func==2)then
write(*,'(A)')'W(x) function is (x-a)**alfa*(b-x)**beta*log(x-a)'
elseif(weight_func==3)then
write(*,'(A)')'W(x) function is (x-a)**alfa*(b-x)**beta*log(b-x) '
elseif(weight_func==4)then
write(*,'(A)')'W(x) function is (x-a)**alfa*(b-x)**beta*log(x-a)*log(b-x)'
endif
endif
call QAWS(func,a,b,alfa,beta,weight_func,epsabs_,epsrel_,result,abserr,neval,ier)
if(verbose_)then
write(*,'(A,2F14.6)')'Endpoints (a,b) =', a,b
write(*,'(A,2F14.6)')'Parameters (alfa,beta) =', alfa,beta
write(*,'(A,F14.6)') 'Estimated integral =', result
write(*,'(A,F14.6)') 'Estimated integral error =', abserr
write(*,'(A,I8)') 'Error return code IER =', ier
else
if(ier/=0)then
write(*,'(A,I8)') 'Error return code IER =', ier
stop
endif
endif
!
!
end select
!
call delete_finter(finterp)
!
!
contains
!
!
subroutine init_finter(func,Xin,Fin,N)
type(finter_type) :: func
real(8) :: xin(:)
real(8) :: fin(size(xin))
integer :: N,Lin
if(func%status)deallocate(func%x,func%f)
Lin=size(xin)
allocate(func%x(Lin),func%f(Lin))
func%X = Xin
func%F = Fin
func%Imax = Lin
func%Imin = 1
func%N = N
func%status=.true.
end subroutine init_finter
!
subroutine delete_finter(func)
type(finter_type) :: func
if(allocated(func%x))deallocate(func%x)
if(allocated(func%f))deallocate(func%f)
func%imax=0
func%imin=0
func%N =0
func%status=.false.
end subroutine delete_finter
!
function func(x)
real(8) :: x
real(8) :: func
real(8) :: dy
integer :: j,k,k0,k1
integer :: n
N=finterp%N !order of polynomial interpolation
func=0d0
j=locate(finterp%X(finterp%Imin:finterp%Imax),x)
!
k=max(j-(N-1)/2,1)
k0=k
if(k0 < finterp%Imin)k0=finterp%Imin
k1=k+N+1
if(k1 > finterp%Imax)then
k1=finterp%Imax
k0=k1-N-1
endif
call polint(finterp%X(k0:k1),finterp%F(k0:k1),x,func,dy)
end function func
end subroutine quad_sample