!****************************************************************************** ! Complementary error function ! From GNU scientific library 0.6, translated from c to fortran. !****************************************************************************** function erfc11(x) result(erfc) use precision implicit none real(wp), intent(in) :: x real(wp) :: erfc real(wp) :: ax,t,num,den integer :: i real(wp), parameter :: p(0:5) = (/ & 0.0000000297886562639399288862e+08_wp, & 0.0000007409740605964741794425e+07_wp, & 0.0000061602098531096305440906e+06_wp, & 0.0005019049726784267463450058e+04_wp, & 0.1275366644729965952479585264e+01_wp, & 0.5641895835477550741253201704e+00_wp /) real(wp), parameter :: q(0:6) = (/ & 0.0000000033690752069827527677e+09_wp, & 0.0000009608965327192787870698e+07_wp, & 0.0001708144074746600431571095e+05_wp, & 0.0120489519278551290360340491e+03_wp, & 0.9396034016235054150430579648e+01_wp, & 0.2260528520767326969591866945e+01_wp, & 0.1000000000000000000000000000e+01_wp /) real(wp), parameter :: c1(20) = (/ & +0.106073416421769980345174155056e+01_wp, & -0.042582445804381043569204735291e+01_wp, & +0.004955262679620434040357683080e+01_wp, & +0.000449293488768382749558001242e+01_wp, & -0.000129194104658496953494224761e+01_wp, & -0.000001836389292149396270416979e+01_wp, & +0.000002211114704099526291538556e+01_wp, & -0.523337485234257134673693179020e-06_wp, & -0.278184788833537885382530989578e-06_wp, & +0.141158092748813114560316684249e-07_wp, & +0.272571296330561699984539141865e-08_wp, & -0.206343904872070629406401492476e-09_wp, & -0.214273991996785367924201401812e-10_wp, & +0.222990255539358204580285098119e-11_wp, & +0.136250074650698280575807934155e-12_wp, & -0.195144010922293091898995913038e-13_wp, & -0.685627169231704599442806370690e-15_wp, & +0.144506492869699938239521607493e-15_wp, & +0.245935306460536488037576200030e-17_wp, & -0.929599561220523396007359328540e-18_wp /) real(wp), parameter :: c2(25) = (/ & +0.044045832024338111077637466616e+01_wp, & -0.143958836762168335790826895326e+00_wp, & +0.044786499817939267247056666937e+00_wp, & -0.013343124200271211203618353102e+00_wp, & +0.003824682739750469767692372556e+00_wp, & -0.001058699227195126547306482530e+00_wp, & +0.000283859419210073742736310108e+00_wp, & -0.000073906170662206760483959432e+00_wp, & +0.000018725312521489179015872934e+00_wp, & -0.462530981164919445131297264430e-05_wp, & +0.111558657244432857487884006422e-05_wp, & -0.263098662650834130067808832725e-06_wp, & +0.607462122724551777372119408710e-07_wp, & -0.137460865539865444777251011793e-07_wp, & +0.305157051905475145520096717210e-08_wp, & -0.665174789720310713757307724790e-09_wp, & +0.142483346273207784489792999706e-09_wp, & -0.300141127395323902092018744545e-10_wp, & +0.622171792645348091472914001250e-11_wp, & -0.126994639225668496876152836555e-11_wp, & +0.255385883033257575402681845385e-12_wp, & -0.506258237507038698392265499770e-13_wp, & +0.989705409478327321641264227110e-14_wp, & -0.190685978789192181051961024995e-14_wp, & +0.350826648032737849245113757340e-15_wp /) real(wp), parameter :: c3(20) = (/ & +0.111684990123545698684297865808e+01_wp, & +0.003736240359381998520654927536e+00_wp, & -0.000916623948045470238763619870e+00_wp, & +0.000199094325044940833965078819e+00_wp, & -0.000040276384918650072591781859e+00_wp, & +0.776515264697061049477127605790e-05_wp, & -0.144464794206689070402099225301e-05_wp, & +0.261311930343463958393485241947e-06_wp, & -0.461833026634844152345304095560e-07_wp, & +0.800253111512943601598732144340e-08_wp, & -0.136291114862793031395712122089e-08_wp, & +0.228570483090160869607683087722e-09_wp, & -0.378022521563251805044056974560e-10_wp, & +0.617253683874528285729910462130e-11_wp, & -0.996019290955316888445830597430e-12_wp, & +0.158953143706980770269506726000e-12_wp, & -0.251045971047162509999527428316e-13_wp, & +0.392607828989125810013581287560e-14_wp, & -0.607970619384160374392535453420e-15_wp, & +0.912600607264794717315507477670e-16_wp /) ax = abs(x) if (ax == 0.0_wp) then stop 'GSL_ERFC: zero argument' elseif (ax <= 1.0_wp) then t = 2.0_wp*ax - 1.0_wp erfc = evalcs(t,c1) elseif (ax <= 5.0_wp) then t = 0.5_wp*(ax-3.0_wp) erfc = evalcs(t,c2) * exp(-x**2) elseif (ax < 10.0_wp) then t = (2.0_wp*x - 15.0_wp)/5.0_wp erfc = evalcs(t,c3) * exp(-x**2)/x else num = p(5) do i=4,0,-1 num = x*num + p(i) enddo den = q(6) do i=5,0,-1 den = x*den + q(i) enddo erfc = num/den * exp(-x**2) endif if (x < 0.0_wp) erfc = 2.0_wp - erfc contains !******************************************************************************* function evalcs(x,cs) implicit none real(wp), intent(in) :: x,cs(:) real(wp) :: evalcs real(wp) :: b0,b1,b2,twox integer :: i b0 = 0.0_wp b1 = 0.0_wp twox = 2.0_wp*x do i=size(cs),1,-1 b2 = b1 b1 = b0 b0 = twox*b1 - b2 + cs(i) enddo evalcs = (b0-b2)/2.0_wp end function evalcs !******************************************************************************* end function erfc11