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