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 |