subroutine quad_func(func,a,b,epsabs,epsrel,&
key,&
inf,&
singular_endpoint,&
singular_points,&
cpole,&
alfa,beta,&
omega,&
weight_func,&
verbose,&
strict,&
result)
interface
function func(x)
real(8) :: x
real(8) :: func
end function func
end interface
!optional variables
real(8) :: a
real(8),optional :: b
real(8),optional :: epsabs
real(8),optional :: epsrel
integer,optional :: key !order of GK rule in QAG
integer,optional :: inf !infinite integration limit in QAGI:
! ! -1(-inf:a);1(a:inf);2(-inf:inf)
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 QAWF,QAWO
integer,optional :: weight_func !if(QAWF,QAWO)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
logical,optional :: strict
real(8) :: result
!actual default variables
real(8) :: epsabs_
real(8) :: epsrel_
logical :: verbose_
logical :: strict_
real(8) :: abserr
integer :: Neval
integer :: Ier,i
integer :: Nsingular_points
character(len=20) :: routine_name
!
!
!
epsabs_ = 1d-12 ; if(present(epsabs))epsabs_=epsabs
epsrel_ = 1d-6 ; if(present(epsrel))epsrel_=epsrel
verbose_=.false.; if(present(verbose))verbose_=verbose
strict_ =.true. ; if(present(strict))strict_=strict
!
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(key).AND.present(inf))stop "ERROR in quad: key & inf"
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(singular_endpoint).AND.present(inf))stop "ERROR in quad: singular_endpoint & inf"
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(cpole).AND.present(inf))stop "ERROR in quad: cpole & inf"
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(omega).AND.present(inf))stop "ERROR in quad: omega & inf"
if(present(singular_points).AND.(present(alfa).OR.present(beta)))stop "ERROR in quad: singular_points & alfa,beta"
if(present(singular_points).AND.present(inf))stop "ERROR in quad: singular_points & inf"
if(present(inf).AND.(present(alfa).OR.present(beta)))stop "ERROR in quad: inf & alfa,beta"
!
if(present(b))then
!QNG>QAGS: 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'
else
!QAGI: need inf
!QAWF: need omega + weight_func=1,2
if(present(inf))then
routine_name='QAGI'
elseif(present(omega))then
routine_name='QAWF'
else
stop 'ERROR in quad: b not present but neither inf (QAGI) nor omega (QAWF) are given. stop.'
endif
endif
!
if(verbose_)write(*,"(A,A)")"QUAD: selected procedure =", routine_name
!
result=0d0
!
select case(routine_name)
case ('QNG')
! call QNG(func,a,b,epsabs_,epsrel_,result,abserr,neval,ier)
! call QAG(func,a,b,epsabs_,epsrel_,1,result,abserr,neval,ier)
call QAGS(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
if(strict_)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>2)then
write(*,'(A,I8)') 'Error return code IER =', ier
if(strict_)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>2)then
write(*,'(A,I8)') 'Error return code IER =', ier
if(strict_)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>2)then
write(*,'(A,I8)') 'Error return code IER =', ier
if(strict_)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>2)then
write(*,'(A,I8)') 'Error return code IER =', ier
if(strict_)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
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>2)then
write(*,'(A,I8)') 'Error return code IER =', ier
if(strict_)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>2)then
write(*,'(A,I8)') 'Error return code IER =', ier
if(strict_)stop
endif
endif
!
!
case ('QAGI')
if(inf<-1.OR.inf>2)stop "ERROR in quad: use QAGI, inf !=[-1,1,2]"
call QAGI(func,a,inf,epsabs_,epsrel_,result,abserr,neval,ier)
if(verbose_)then
if(inf==-1)then
write(*,'(A,F14.6)')'Integration interval is (-infty:a)',a
elseif(inf==1)then
write(*,'(A,F14.6)')'Integration interval is (a:infty) ',a
elseif(inf==2)then
write(*,'(A)')'Integration interval is (-infty:infty)'
endif
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>2)then
write(*,'(A,I8)') 'Error return code IER =', ier
if(strict_)stop
endif
endif
!
!
case ('QAWF')
if(weight_func<1.OR.weight_func>2)stop "ERROR in quad: use QAWF, 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 QAWF(func,a,omega,weight_func,epsabs_,result,abserr,neval,ier)
if(verbose_)then
write(*,'(A,2F14.6)')'Endpoints (a,infty) =', a
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>2)then
write(*,'(A,I8)') 'Error return code IER =', ier
if(strict_)stop
endif
endif
!
!
end select
end subroutine quad_func