!******************************************************************************* ! Complementary Error Function ! Naval Surface Warfare Center, Dahlgren, Virginia, April 1992 !******************************************************************************* function erfc10(x) result(erfc) use precision, only: wp implicit none real(wp), intent(in) :: x real(wp) :: erfc real(wp) :: ax,t,w integer :: i real(wp), parameter :: a(21) = (/ & 0.1283791670955125738961589031215e+00_wp, & -0.3761263890318375246320529677070e+00_wp, & 0.1128379167095512573896158902931e+00_wp, & -0.2686617064513125175943235372542e-01_wp, & 0.5223977625442187842111812447877e-02_wp, & -0.8548327023450852832540164081187e-03_wp, & 0.1205533298178966425020717182498e-03_wp, & -0.1492565035840625090430728526820e-04_wp, & 0.1646211436588924261080723578109e-05_wp, & -0.1636584469123468757408968429674e-06_wp, & 0.1480719281587021715400818627811e-07_wp, & -0.1229055530145120140800510155331e-08_wp, & 0.9422759058437197017313055084212e-10_wp, & -0.6711366740969385085896257227159e-11_wp, & 0.4463222608295664017461758843550e-12_wp, & -0.2783497395542995487275065856998e-13_wp, & 0.1634095572365337143933023780777e-14_wp, & -0.9052845786901123985710019387938e-16_wp, & 0.4708274559689744439341671426731e-17_wp, & -0.2187159356685015949749948252160e-18_wp, & 0.7043407712019701609635599701333e-20_wp /) ax = abs(x) ! abs(x) <= 1 if (ax <= 1.0_wp) then t = x**2 w = a(21) do i=20,1,-1 w = t*w + a(i) enddo erfc = 0.5_wp + (0.5_wp - x*(1.0_wp + w)) return endif ! x < -1 if (x <= 0.0_wp) then erfc = 2.0_wp if (x < -8.3_wp) return t = x**2 erfc = 2.0_wp - exp(-t)*erfc0(ax) return endif ! x > 1 erfc = 0.0_wp if (x > 100.0_wp) return t = x**2 if (t > -exparg(1)) return erfc = exp(-t) * erfc0(x) contains !******************************************************************************* ! Evaluation of exp(x**2)*erfc(x) for x >= 1 ! Written by Alfred H. Morris, Jr. ! Naval Surface Warfare Center, Dahlgren, Virginia, April 1992 !******************************************************************************* function erfc0(x) use precision, only: wp implicit none real(wp), intent(in) :: x real(wp) :: erfc0 real(wp) :: t,u,v,z ! rpinv = 1/sqrt(pi) real(wp), parameter :: rpinv = 0.56418958354775628694807945156077259_wp real(wp), parameter :: p0 = 0.16506148041280876191828601e-03_wp, & p1 = 0.15471455377139313353998665e-03_wp, & p2 = 0.44852548090298868465196794e-04_wp, & p3 = -0.49177280017226285450486205e-05_wp, & p4 = -0.69353602078656412367801676e-05_wp, & p5 = -0.20508667787746282746857743e-05_wp, & p6 = -0.28982842617824971177267380e-06_wp, & p7 = -0.17272433544836633301127174e-07_wp real(wp), parameter :: q0 = 1.00000000000000000000000000e+00_wp, & q1 = 0.16272656776533322859856317e+01_wp, & q2 = 0.12040996037066026106794322e+01_wp, & q3 = 0.52400246352158386907601472e+00_wp, & q4 = 0.14497345252798672362384241e+00_wp, & q5 = 0.25592517111042546492590736e-01_wp, & q6 = 0.26869088293991371028123158e-02_wp, & q7 = 0.13133767840925681614496481e-03_wp real(wp), parameter :: r0 = 0.145589721275038539045668824025e+00_wp, & r1 = -0.273421931495426482902320421863e+00_wp, & r2 = 0.226008066916621506788789064272e+00_wp, & r3 = -0.163571895523923805648814425592e+00_wp, & r4 = 0.102604312032193978662297299832e+00_wp, & r5 = -0.548023266949835519254211506880e-01_wp, & r6 = 0.241432239725390106956523668160e-01_wp, & r7 = -0.822062115403915116036874169600e-02_wp, & r8 = 0.180296241564687154310619200000e-02_wp real(wp), parameter :: a0 = -0.45894433406309678202825375e-03_wp, & a1 = -0.12281298722544724287816236e-01_wp, & a2 = -0.91144359512342900801764781e-01_wp, & a3 = -0.28412489223839285652511367e-01_wp, & a4 = 0.14083827189977123530129812e+01_wp, & a5 = 0.11532175281537044570477189e+01_wp, & a6 = -0.72170903389442152112483632e+01_wp, & a7 = -0.19685597805218214001309225e+01_wp, & a8 = 0.93846891504541841150916038e+01_wp real(wp), parameter :: b0 = 1.00000000000000000000000000e+00_wp, & b1 = 0.25136329960926527692263725e+02_wp, & b2 = 0.15349442087145759184067981e+03_wp, & b3 = -0.29971215958498680905476402e+03_wp, & b4 = -0.33876477506888115226730368e+04_wp, & b5 = 0.28301829314924804988873701e+04_wp, & b6 = 0.22979620942196507068034887e+05_wp, & b7 = -0.24280681522998071562462041e+05_wp, & b8 = -0.36680620673264731899504580e+05_wp, & b9 = 0.42278731622295627627042436e+05_wp, & b10= 0.28834257644413614344549790e+03_wp, & b11= 0.70226293775648358646587341e+03_wp real(wp), parameter :: c0 = -0.7040906288250128001000086e-04_wp, & c1 = -0.3858822461760510359506941e-02_wp, & c2 = -0.7708202127512212359395078e-01_wp, & c3 = -0.6713655014557429480440263e+00_wp, & c4 = -0.2081992124162995545731882e+01_wp, & c5 = 0.2898831421475282558867888e+01_wp, & c6 = 0.2199509380600429331650192e+02_wp, & c7 = 0.2907064664404115316722996e+01_wp, & c8 = -0.4766208741588182425380950e+02_wp real(wp), parameter :: d0 = 1.0000000000000000000000000e+00_wp, & d1 = 0.5238852785508439144747174e+02_wp, & d2 = 0.9646843357714742409535148e+03_wp, & d3 = 0.7007152775135939601804416e+04_wp, & d4 = 0.8515386792259821780601162e+04_wp, & d5 = -0.1002360095177164564992134e+06_wp, & d6 = -0.2065250031331232815791912e+06_wp, & d7 = 0.5695324805290370358175984e+06_wp, & d8 = 0.6589752493461331195697873e+06_wp, & d9 = -0.1192930193156561957631462e+07_wp real(wp), parameter :: e0 = 0.540464821348814822409610122136e+00_wp, & e1 = -0.261515522487415653487049835220e-01_wp, & e2 = -0.288573438386338758794591212600e-02_wp, & e3 = -0.529353396945788057720258856000e-03_wp real(wp), parameter :: s0 = 1.00000000000000000000e+00_wp, & s1 = -0.50000000000000000000e+00_wp, & s2 = 0.75000000000000000000e+00_wp, & s3 = -0.18750000000000000000e+01_wp, & s4 = 0.65625000000000000000e+01_wp, & s5 = -0.29531250000000000000e+02_wp, & s6 = 0.16242187500000000000e+03_wp, & s7 = -0.10557421875000000000e+04_wp, & s8 = 0.79180664062500000000e+04_wp, & s9 = -0.67303564453125000000e+05_wp, & s10 = 0.63938386230468750000e+06_wp, & s11 = -0.67135305541992187500e+07_wp, & s12 = 0.77205601373291015625e+08_wp if (x < 1.0_wp) then stop 'ERFC0: argument x < 1' elseif (x <= 2.0_wp) then u = ((((((p7*x + p6)*x + p5)*x + p4)*x + p3)*x + p2)*x + p1)*x + p0 v = ((((((q7*x + q6)*x + q5)*x + q4)*x + q3)*x + q2)*x + q1)*x + q0 t = (x-3.75_wp) / (x+3.75_wp) erfc0 = (((((((((u/v)*t + r8)*t + r7)*t + r6)*t + r5) & *t + r4)*t + r3)*t + r2)*t + r1)*t + r0 elseif (x <= 4.0_wp) then z = 1.0_wp/(2.5_wp + x*x) u = (((((((a8*z + a7)*z + a6)*z + a5)*z + a4) & *z + a3)*z + a2)*z + a1)*z + a0 v = ((((((((((b11*z + b10)*z + b9)*z + b8)*z + b7)*z + b6) & *z + b5)*z + b4)*z + b3)*z + b2)*z + b1)*z + b0 t = 13.0_wp*z - 1.0_wp erfc0 = ((((u/v)*t + e2)*t + e1)*t + e0) / x elseif (x < 50.0_wp) then z = 1.0_wp/(2.5_wp + x*x) u = (((((((c8*z + c7)*z + c6)*z + c5)*z + c4) & *z + c3)*z + c2)*z + c1)*z + c0 v = ((((((((d9*z + d8)*z + d7)*z + d6)*z + d5) & *z + d4)*z + d3)*z + d2)*z + d1)*z + d0 t = 13.0_wp*z - 1.0_wp erfc0 = (((((u/v)*t + e3)*t + e2)*t + e1)*t + e0) / x else t = 1.0_wp/(x*x) z = (((((((((((s12*t + s11)*t + s10)*t + s9)*t + s8)*t + s7)*t + s6) & *t + s5)*t + s4)*t + s3)*t + s2)*t + s1)*t + s0 erfc0 = rpinv*(z/x) endif end function erfc0 !******************************************************************************* ! exparg(l) = ! if l == 0: the largest positive w for which exp(w) can be computed. ! if l /= 0: the largest negative w for which the computed value of ! exp(w) is nonzero. ! ! Note: only an approximate value for exparg(l) is needed. !******************************************************************************* function exparg(l) use precision, only: wp implicit none integer, intent(in) :: l real(wp) :: exparg real(wp), parameter :: one = 1.0_wp if (l == 0) then exparg = log(huge(one)) else exparg = log(tiny(one)) endif end function exparg !******************************************************************************* end function erfc10