| 1 |
! ********************************************************************** |
| 2 |
! erfc.F used with permission from Naval Surface Warfare Group |
| 3 |
!****************************************************************************** |
| 4 |
function erfc (x) |
| 5 |
|
| 6 |
use definitions |
| 7 |
|
| 8 |
!! ERFC: the complementary error function |
| 9 |
|
| 10 |
real( kind = DP ) :: erfc |
| 11 |
real( kind = DP ) :: an, ax, c, eps, rpinv, t, x, w |
| 12 |
real( kind = DP ) :: a(21), b(44), e(44) |
| 13 |
real( kind = DP ) :: dpmpar, dcsevl |
| 14 |
! |
| 15 |
! rpinv = 1/sqrt(pi) |
| 16 |
! |
| 17 |
data rpinv /.56418958354775628694807945156077259e0_dp/ |
| 18 |
! |
| 19 |
data a(1) / .1283791670955125738961589031215e+00_dp/, & |
| 20 |
a(2) /-.3761263890318375246320529677070e+00_dp/, & |
| 21 |
a(3) / .1128379167095512573896158902931e+00_dp/, & |
| 22 |
a(4) /-.2686617064513125175943235372542e-01_dp/, & |
| 23 |
a(5) / .5223977625442187842111812447877e-02_dp/, & |
| 24 |
a(6) /-.8548327023450852832540164081187e-03_dp/, & |
| 25 |
a(7) / .1205533298178966425020717182498e-03_dp/, & |
| 26 |
a(8) /-.1492565035840625090430728526820e-04_dp/, & |
| 27 |
a(9) / .1646211436588924261080723578109e-05_dp/, & |
| 28 |
a(10) /-.1636584469123468757408968429674e-06_dp/ |
| 29 |
data a(11) / .1480719281587021715400818627811e-07_dp/, & |
| 30 |
a(12) /-.1229055530145120140800510155331e-08_dp/, & |
| 31 |
a(13) / .9422759058437197017313055084212e-10_dp/, & |
| 32 |
a(14) /-.6711366740969385085896257227159e-11_dp/, & |
| 33 |
a(15) / .4463222608295664017461758843550e-12_dp/, & |
| 34 |
a(16) /-.2783497395542995487275065856998e-13_dp/, & |
| 35 |
a(17) / .1634095572365337143933023780777e-14_dp/, & |
| 36 |
a(18) /-.9052845786901123985710019387938e-16_dp/, & |
| 37 |
a(19) / .4708274559689744439341671426731e-17_dp/, & |
| 38 |
a(20) /-.2187159356685015949749948252160e-18_dp/, & |
| 39 |
a(21) / .7043407712019701609635599701333e-20_dp/ |
| 40 |
! |
| 41 |
data b(1) / .610143081923200417926465815756e+00_dp/, & |
| 42 |
b(2) /-.434841272712577471828182820888e+00_dp/, & |
| 43 |
b(3) / .176351193643605501125840298123e+00_dp/, & |
| 44 |
b(4) /-.607107956092494148600512158250e-01_dp/, & |
| 45 |
b(5) / .177120689956941144861471411910e-01_dp/, & |
| 46 |
b(6) /-.432111938556729381859986496800e-02_dp/, & |
| 47 |
b(7) / .854216676887098678819832055000e-03_dp/, & |
| 48 |
b(8) /-.127155090609162742628893940000e-03_dp/, & |
| 49 |
b(9) / .112481672436711894688470720000e-04_dp/, & |
| 50 |
b(10) / .313063885421820972630152000000e-06_dp/ |
| 51 |
data b(11) /-.270988068537762022009086000000e-06_dp/, & |
| 52 |
b(12) / .307376227014076884409590000000e-07_dp/, & |
| 53 |
b(13) / .251562038481762293731400000000e-08_dp/, & |
| 54 |
b(14) /-.102892992132031912759000000000e-08_dp/, & |
| 55 |
b(15) / .299440521199499393630000000000e-10_dp/, & |
| 56 |
b(16) / .260517896872669362900000000000e-10_dp/, & |
| 57 |
b(17) /-.263483992417196938600000000000e-11_dp/, & |
| 58 |
b(18) /-.643404509890636443000000000000e-12_dp/, & |
| 59 |
b(19) / .112457401801663447000000000000e-12_dp/, & |
| 60 |
b(20) / .172815333899860980000000000000e-13_dp/ |
| 61 |
data b(21) /-.426410169494237500000000000000e-14_dp/, & |
| 62 |
b(22) /-.545371977880191000000000000000e-15_dp/, & |
| 63 |
b(23) / .158697607761671000000000000000e-15_dp/, & |
| 64 |
b(24) / .208998378443340000000000000000e-16_dp/, & |
| 65 |
b(25) /-.590052686940900000000000000000e-17_dp/, & |
| 66 |
b(26) /-.941893387554000000000000000000e-18_dp/, & |
| 67 |
b(27) / .214977356470000000000000000000e-18_dp/, & |
| 68 |
b(28) / .466609850080000000000000000000e-19_dp/, & |
| 69 |
b(29) /-.724301186200000000000000000000e-20_dp/, & |
| 70 |
b(30) /-.238796682400000000000000000000e-20_dp/ |
| 71 |
data b(31) / .191177535000000000000000000000e-21_dp/, & |
| 72 |
b(32) / .120482568000000000000000000000e-21_dp/, & |
| 73 |
b(33) /-.672377000000000000000000000000e-24_dp/, & |
| 74 |
b(34) /-.574799700000000000000000000000e-23_dp/, & |
| 75 |
b(35) /-.428493000000000000000000000000e-24_dp/, & |
| 76 |
b(36) / .244856000000000000000000000000e-24_dp/, & |
| 77 |
b(37) / .437930000000000000000000000000e-25_dp/, & |
| 78 |
b(38) /-.815100000000000000000000000000e-26_dp/, & |
| 79 |
b(39) /-.308900000000000000000000000000e-26_dp/, & |
| 80 |
b(40) / .930000000000000000000000000000e-28_dp/ |
| 81 |
data b(41) / .174000000000000000000000000000e-27_dp/, & |
| 82 |
b(42) / .160000000000000000000000000000e-28_dp/, & |
| 83 |
b(43) /-.800000000000000000000000000000e-29_dp/, & |
| 84 |
b(44) /-.200000000000000000000000000000e-29_dp/ |
| 85 |
! |
| 86 |
data e(1) / .107797785207238315116833591035e+01_dp/, & |
| 87 |
e(2) /-.265598904091486733721465009040e-01_dp/, & |
| 88 |
e(3) /-.148707314669809950960504633300e-02_dp/, & |
| 89 |
e(4) /-.138040145414143859607708920000e-03_dp/, & |
| 90 |
e(5) /-.112803033322874914985073660000e-04_dp/, & |
| 91 |
e(6) /-.117286984274372522405373900000e-05_dp/, & |
| 92 |
e(7) /-.103476150393304615537382000000e-06_dp/, & |
| 93 |
e(8) /-.118991140858924382544470000000e-07_dp/, & |
| 94 |
e(9) /-.101622254498949864047600000000e-08_dp/, & |
| 95 |
e(10) /-.137895716146965692169000000000e-09_dp/ |
| 96 |
data e(11) /-.936961303373730333500000000000e-11_dp/, & |
| 97 |
e(12) /-.191880958395952534900000000000e-11_dp/, & |
| 98 |
e(13) /-.375730172019937070000000000000e-13_dp/, & |
| 99 |
e(14) /-.370537260269833570000000000000e-13_dp/, & |
| 100 |
e(15) / .262756542349037100000000000000e-14_dp/, & |
| 101 |
e(16) /-.112132287643793300000000000000e-14_dp/, & |
| 102 |
e(17) / .184136028922538000000000000000e-15_dp/, & |
| 103 |
e(18) /-.491302565748860000000000000000e-16_dp/, & |
| 104 |
e(19) / .107044551673730000000000000000e-16_dp/, & |
| 105 |
e(20) /-.267189366240500000000000000000e-17_dp/ |
| 106 |
data e(21) / .649326867976000000000000000000e-18_dp/, & |
| 107 |
e(22) /-.165399353183000000000000000000e-18_dp/, & |
| 108 |
e(23) / .426056266040000000000000000000e-19_dp/, & |
| 109 |
e(24) /-.112558407650000000000000000000e-19_dp/, & |
| 110 |
e(25) / .302561744800000000000000000000e-20_dp/, & |
| 111 |
e(26) /-.829042146000000000000000000000e-21_dp/, & |
| 112 |
e(27) / .231049558000000000000000000000e-21_dp/, & |
| 113 |
e(28) /-.654695110000000000000000000000e-22_dp/, & |
| 114 |
e(29) / .188423140000000000000000000000e-22_dp/, & |
| 115 |
e(30) /-.550434100000000000000000000000e-23_dp/ |
| 116 |
data e(31) / .163095000000000000000000000000e-23_dp/, & |
| 117 |
e(32) /-.489860000000000000000000000000e-24_dp/, & |
| 118 |
e(33) / .149054000000000000000000000000e-24_dp/, & |
| 119 |
e(34) /-.459220000000000000000000000000e-25_dp/, & |
| 120 |
e(35) / .143180000000000000000000000000e-25_dp/, & |
| 121 |
e(36) /-.451600000000000000000000000000e-26_dp/, & |
| 122 |
e(37) / .144000000000000000000000000000e-26_dp/, & |
| 123 |
e(38) /-.464000000000000000000000000000e-27_dp/, & |
| 124 |
e(39) / .151000000000000000000000000000e-27_dp/, & |
| 125 |
e(40) /-.500000000000000000000000000000e-28_dp/ |
| 126 |
data e(41) / .170000000000000000000000000000e-28_dp/, & |
| 127 |
e(42) /-.600000000000000000000000000000e-29_dp/, & |
| 128 |
e(43) / .200000000000000000000000000000e-29_dp/, & |
| 129 |
e(44) /-.100000000000000000000000000000e-29_dp/ |
| 130 |
! |
| 131 |
eps = epsilon ( eps ) |
| 132 |
! |
| 133 |
! abs(x) <= 1 |
| 134 |
! |
| 135 |
ax = abs(x) |
| 136 |
if (ax > 1.0_dp) go to 20 |
| 137 |
t = x*x |
| 138 |
w = a(21) |
| 139 |
do i = 1,20 |
| 140 |
k = 21 - i |
| 141 |
w = t*w + a(k) |
| 142 |
end do |
| 143 |
erfc = 0.5_dp + (0.5_dp - x*(1.0_dp + w)) |
| 144 |
return |
| 145 |
! |
| 146 |
! 1 < abs(x) < 2 |
| 147 |
! |
| 148 |
20 if (ax >= 2.0_dp) go to 30 |
| 149 |
n = 44 |
| 150 |
if (eps >= 1.0e-20_dp) n = 30 |
| 151 |
t = (ax - 3.75_dp)/(ax + 3.75_dp) |
| 152 |
erfc = csevl(t, b, n) |
| 153 |
21 erfc = exp(-x*x) * erfc |
| 154 |
if (x < 0.0_dp) erfc = 2.0_dp - erfc |
| 155 |
return |
| 156 |
! |
| 157 |
! 2 < dabs(x) < 12 |
| 158 |
! |
| 159 |
30 if (x < -9.0_dp) go to 60 |
| 160 |
if (x >= 12.0_dp) go to 40 |
| 161 |
n = 44 |
| 162 |
if (eps >= 1.0e-20_dp) n = 25 |
| 163 |
t = (1.0_dp/x)**2 |
| 164 |
w = (10.5_dp*t - 1.0_dp)/(2.5_dp*t + 1.0_dp) |
| 165 |
erfc = csevl(w, e, n) / ax |
| 166 |
go to 21 |
| 167 |
! |
| 168 |
! x >= 12 |
| 169 |
! |
| 170 |
40 if (x > 50.0_dp) go to 70 |
| 171 |
t = (1.0_dp/x)**2 |
| 172 |
an = -0.5_dp |
| 173 |
c = 0.5_dp |
| 174 |
w = 0.0_dp |
| 175 |
50 c = c + 1.0_dp |
| 176 |
an = - c*an*t |
| 177 |
w = w + an |
| 178 |
if (abs(an) > eps) go to 50 |
| 179 |
w = (-0.5_dp + w)*t + 1.0_dp |
| 180 |
erfc = exp(-x*x) * ((rpinv*w)/ax) |
| 181 |
return |
| 182 |
! |
| 183 |
! limit value for large negative x |
| 184 |
! |
| 185 |
60 erfc = 2.0_dp |
| 186 |
return |
| 187 |
! |
| 188 |
! limit value for large positive x |
| 189 |
! |
| 190 |
70 erfc = 0.0_dp |
| 191 |
return |
| 192 |
end function erfc |
| 193 |
|
| 194 |
function csevl (x, a, n) |
| 195 |
|
| 196 |
!! CSEVL: evaluate the n term chebyshev series a at x. |
| 197 |
!! only half of the first coefficient is used. |
| 198 |
|
| 199 |
use definitions |
| 200 |
|
| 201 |
real( kind = DP ) :: csevl |
| 202 |
real( kind = DP ) :: a(n), x, x2, s0, s1, s2 |
| 203 |
|
| 204 |
if (n .le. 1) then |
| 205 |
csevl = 0.5_dp * a(1) |
| 206 |
return |
| 207 |
else |
| 208 |
|
| 209 |
x2 = x + x |
| 210 |
s0 = a(n) |
| 211 |
s1 = 0.0_dp |
| 212 |
do i = 2,n |
| 213 |
s2 = s1 |
| 214 |
s1 = s0 |
| 215 |
k = n - i + 1 |
| 216 |
s0 = x2*s1 - s2 + a(k) |
| 217 |
end do |
| 218 |
csevl = 0.5_dp * (s0 - s2) |
| 219 |
return |
| 220 |
|
| 221 |
endif |
| 222 |
end function csevl |
