This is the 8-point formula from the original SLATEC routine DGAUS8.
replaced coefficients with high-precision ones from: http://processingjs.nihongoresources.com/bezierinfo/legendre-gauss-values.php
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(integration_class_1d), | intent(inout) | :: | me | |||
| real(kind=wp), | intent(in) | :: | x | |||
| real(kind=wp), | intent(in) | :: | h |
function g8(me, x, h) result(f)
implicit none
class(integration_class_1d),intent(inout) :: me
real(wp), intent(in) :: x
real(wp), intent(in) :: h
real(wp) :: f
!> abscissae:
real(wp),parameter :: x1 = 0.18343464249564980493947614236018398066675781291297378231718847&
&369920447422154211411606822371112335374526765876&
&428676660891960125238768656837885699951606635681&
&044755516171385019663858107642055323708826547494&
&928123149612477646193635627706457164566131594051&
&34052985058171969174306064445289638150514997832_wp
real(wp),parameter :: x2 = 0.52553240991632898581773904918924634904196424312039285775085709&
&927245482076856127252396140019363198206190968292&
&482526085071087937666387799398053953036682536311&
&190182730324023600607174700061279014795875767562&
&412888953366196435283308256242634705401842246036&
&88817537938539658502113876953598879150514997832_wp
real(wp),parameter :: x3 = 0.79666647741362673959155393647583043683717173161596483207017029&
&503921730567647309214715192729572593901919745345&
&309730926536564949170108596027725620746216896761&
&539350162903423256455826342053015458560600957273&
&426035574157612651404288519573419337108037227831&
&36113628137267630651413319993338002150514997832_wp
real(wp),parameter :: x4 = 0.96028985649753623168356086856947299042823523430145203827163977&
&737242489774341928443943895926331226831042439281&
&729417621023895815521712854793736422049096997004&
&339826183266373468087812635533469278673596634808&
&705975425476039293185338665681328688426134748962&
&8923208763998895240977248938732425615051499783203_wp
!> weights:
real(wp),parameter :: w1 = 0.36268378337836198296515044927719561219414603989433054052482306&
&756668673472390667732436604208482850955025876992&
&629670655292582155698951738449955760078620768427&
&783503828625463057710075533732697147148942683287&
&804318227790778467229655355481996014024877675059&
&28976560993309027632737537826127502150514997832_wp
real(wp),parameter :: w2 = 0.31370664587788728733796220198660131326032899900273493769026394&
&507495627194217349696169807623392855604942757464&
&107780861624724683226556160568906242764697589946&
&225031187765625594632872220215204316264677947216&
&038226012952768986525097231851579983531560624197&
&51736972560423953923732838789657919150514997832_wp
real(wp),parameter :: w3 = 0.22238103445337447054435599442624088443013087005124956472590928&
&929361681457044904085365314237719792784215926610&
&121221812311143757985257224193818266745320905779&
&086132895368404027893986488760043856972021574820&
&632532471955902286315706513199655897335454406059&
&52819880671616779621183704306688233150514997832_wp
real(wp),parameter :: w4 = 0.10122853629037625915253135430996219011539409105168495705900369&
&806474017876347078486028273930404500655815438933&
&141326670771549403089234876787319730411360735846&
&905332088240507319763065757292054679614357794675&
&524923287300550259929540899466768105108107294683&
&66466585774650346143712142008566866150514997832_wp
f = h * ( w1*( me%fun(x-x1*h) + me%fun(x+x1*h)) + &
w2*( me%fun(x-x2*h) + me%fun(x+x2*h)) + &
w3*( me%fun(x-x3*h) + me%fun(x+x3*h)) + &
w4*( me%fun(x-x4*h) + me%fun(x+x4*h)) )
end function g8