Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2022 The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : ! from Frank Timmes' site, http://www.cococubed.com/code_pages/fermi_dirac.shtml
21 :
22 : !..routine dfermi gets the fermi-dirac functions and their derivatives
23 : !..routine fdfunc1 forms the integrand of the fermi-dirac functions
24 : !..routine fdfunc2 same as fdfunc but with the change of variable z**2=x
25 : !..routine dqleg010 does 10 point gauss-legendre integration 9 fig accuracy
26 : !..routine dqleg020 does 20 point gauss-legendre integration 14 fig accuracy
27 : !..routine dqleg040 does 40 point gauss-legendre integration 18 fig accuracy
28 : !..routine dqleg080 does 80 point gauss-legendre integration 32 fig accuracy
29 : !..routine dqlag010 does 10 point gauss-laguerre integration 9 fig accuracy
30 : !..routine dqlag020 does 20 point gauss-laguerre integration 14 fig accuracy
31 : !..routine dqlag040 does 40 point gauss-laguerre integration 18 fig accuracy
32 : !..routine dqlag080 does 80 point gauss-laguerre integration 32 fig accuracy
33 :
34 : module gauss_fermi
35 :
36 : use const_def, only: dp
37 : use math_lib
38 :
39 : implicit none
40 :
41 : private
42 :
43 : public :: dfermi
44 :
45 : contains
46 :
47 112 : subroutine dfermi(dk,eta,theta,fd,fdeta,fdtheta)
48 : !..this routine computes the fermi-dirac integrals of
49 : !..index dk, with degeneracy parameter eta and relativity parameter theta.
50 : !..input is dk the real(dp) index of the fermi-dirac function,
51 : !..eta the degeneracy parameter, and theta the relativity parameter.
52 : !.. theta = (k * T)/(mass_electron * c^2), k = Boltzmann const.
53 : !..the output is fd is computed by applying three 10-point
54 : !..gauss-legendre and one 10-point gauss-laguerre rules over
55 : !..four appropriate subintervals. the derivative with respect to eta is
56 : !..output in fdeta, and the derivative with respct to theta is in fdtheta.
57 : !..within each subinterval the fd kernel.
58 : !..
59 : !..this routine delivers at least 9, and up to 32, figures of accuracy
60 : !..
61 : !..reference: j.m. aparicio, apjs 117, 632 1998
62 :
63 :
64 : !..declare
65 : ! external fdfunc1, fdfunc2
66 : real(dp) :: dk,eta,theta,fd,fdeta,fdtheta
67 : real(dp) :: d,sg,a1,b1,c1,a2,b2,c2,d2,e2,a3,b3,c3,d3,e3
68 448 : real(dp) :: eta1,xi,xi2,x1,x2,x3,s1,s2,s3,s12,par(3)
69 224 : real(dp) :: res1,dres1,ddres1,res2,dres2,ddres2
70 112 : real(dp) :: res3,dres3,ddres3,res4,dres4,ddres4
71 :
72 : ! parameters defining the location of the breakpoints for the
73 : ! subintervals of integration:
74 : data d / 3.3609d0 /
75 : data sg / 9.1186d-2 /
76 : data a1 / 6.7774d0 /
77 : data b1 / 1.1418d0 /
78 : data c1 / 2.9826d0 /
79 : data a2 / 3.7601d0 /
80 : data b2 / 9.3719d-2 /
81 : data c2 / 2.1063d-2 /
82 : data d2 / 3.1084d1 /
83 : data e2 / 1.0056d0 /
84 : data a3 / 7.5669d0 /
85 : data b3 / 1.1695d0 /
86 : data c3 / 7.5416d-1 /
87 : data d3 / 6.6558d0 /
88 : data e3 /-1.2819d-1 /
89 :
90 :
91 : ! integrand parameters:
92 112 : par(1)=dk
93 112 : par(2)=eta
94 112 : par(3)=theta
95 :
96 :
97 : ! definition of xi:
98 112 : eta1=sg*(eta-d)
99 112 : if (eta1<=5.d1) then
100 112 : xi=log(1.d0+exp(eta1))/sg
101 : else
102 : xi=eta-d
103 : end if
104 112 : xi2=xi*xi
105 :
106 : ! definition of the x_i:
107 112 : x1=(a1 + b1*xi + c1*xi2)/(1.d0+c1*xi)
108 112 : x2=(a2 + b2*xi + c2*d2*xi2)/(1.d0 + e2*xi + c2*xi2)
109 112 : x3=(a3 + b3*xi + c3*d3*xi2)/(1.d0 + e3*xi + c3*xi2)
110 :
111 : ! breakpoints:
112 112 : s1=x1-x2
113 112 : s2=x1
114 112 : s3=x1+x3
115 112 : s12=sqrt(s1)
116 :
117 : ! quadrature integrations:
118 :
119 : ! 9 significant figure accuracy
120 : ! call dqleg010(fdfunc2, 0.d0, s12, res1, dres1, ddres1, par,3)
121 : ! call dqleg010(fdfunc1, s1, s2, res2, dres2, ddres2, par,3)
122 : ! call dqleg010(fdfunc1, s2, s3, res3, dres3, ddres3, par,3)
123 : ! call dqlag010(fdfunc1, s3, 1.d0, res4, dres4, ddres4, par,3)
124 :
125 :
126 : ! 14 significant figure accuracy
127 112 : call dqleg020(fdfunc2, 0.d0, s12, res1, dres1, ddres1, par,3)
128 112 : call dqleg020(fdfunc1, s1, s2, res2, dres2, ddres2, par,3)
129 112 : call dqleg020(fdfunc1, s2, s3, res3, dres3, ddres3, par,3)
130 112 : call dqlag020(fdfunc1, s3, 1.d0, res4, dres4, ddres4, par,3)
131 :
132 :
133 : ! 18 significant figure accuracy
134 : ! call dqleg040(fdfunc2, 0.d0, s12, res1, dres1, ddres1, par,3)
135 : ! call dqleg040(fdfunc1, s1, s2, res2, dres2, ddres2, par,3)
136 : ! call dqleg040(fdfunc1, s2, s3, res3, dres3, ddres3, par,3)
137 : ! call dqlag040(fdfunc1, s3, 1.d0, res4, dres4, ddres4, par,3)
138 :
139 : ! 32 significant figure accuracy
140 : ! call dqleg080(fdfunc2, 0.d0, s12, res1, dres1, ddres1, par,3)
141 : ! call dqleg080(fdfunc1, s1, s2, res2, dres2, ddres2, par,3)
142 : ! call dqleg080(fdfunc1, s2, s3, res3, dres3, ddres3, par,3)
143 : ! call dqlag080(fdfunc1, s3, 1.d0, res4, dres4, ddres4, par,3)
144 :
145 :
146 : !..sum the contributions
147 112 : fd = res1 + res2 + res3 + res4
148 112 : fdeta = dres1 + dres2 + dres3 + dres4
149 112 : fdtheta = ddres1 + ddres2 + ddres3 + ddres4
150 112 : return
151 : end subroutine dfermi
152 :
153 6720 : subroutine fdfunc1(x,par,n,fd,fdeta,fdtheta)
154 : !..
155 : !..forms the fermi-dirac integrand and its derivatives with eta and theta.
156 : !..on input x is the integration variable, par(1) is the real(dp)
157 : !..index, par(2) is the degeneravy parameter, and par(3) is the relativity
158 : !..parameter. on output fd is the integrand, fdeta is the derivative
159 : !..with respect to eta, and fdtheta is the derivative with respect to theta.
160 : !..
161 : !..declare
162 : integer :: n
163 6720 : real(dp) :: x,par(n),dk,eta,theta,fd,fdeta,fdtheta
164 6720 : real(dp) :: factor,dxst,denom,denom2,xdk,xdkp1
165 :
166 : !..initialize
167 6720 : dk = par(1)
168 6720 : eta = par(2)
169 6720 : theta = par(3)
170 6720 : xdk = pow(x,dk)
171 6720 : xdkp1 = x * xdk
172 6720 : dxst = sqrt(1.0d0 + 0.5d0*x*theta)
173 :
174 : ! avoid overflow in the exponentials at large x
175 6720 : if ((x-eta) < 1.0d2) then
176 5080 : factor = exp(x-eta)
177 5080 : denom = factor + 1.0d0
178 5080 : fd = xdk * dxst / denom
179 5080 : fdeta = fd * factor / denom
180 5080 : denom2 = 4.0d0 * dxst * denom
181 5080 : fdtheta = xdkp1 / denom2
182 :
183 : else
184 1640 : factor = exp(eta-x)
185 1640 : fd = xdk * dxst * factor
186 1640 : fdeta = fd
187 1640 : denom2 = 4.0d0 * dxst
188 1640 : fdtheta = xdkp1/denom2 * factor
189 : end if
190 :
191 6720 : return
192 : end subroutine fdfunc1
193 :
194 2240 : subroutine fdfunc2(x,par,n,fd,fdeta,fdtheta)
195 : !..
196 : !..forms the fermi-dirac integrand and its derivatives with eta and theta,
197 : !..when the z**2=x variable change has been made.
198 : !..on input x is the integration variable, par(1) is the real(dp)
199 : !..index, par(2) is the degeneravy parameter, and par(3) is the relativity
200 : !..parameter. on output fd is the integrand, fdeta is the derivative
201 : !..with respect to eta, and fdtheta is the derivative with respect to theta.
202 : !..
203 : !..declare
204 : integer :: n
205 2240 : real(dp) :: x,par(n),dk,eta,theta,fd,fdeta,fdtheta
206 2240 : real(dp) :: factor,dxst,denom,denom2,xdk,xdkp1,xsq
207 :
208 2240 : dk = par(1)
209 2240 : eta = par(2)
210 2240 : theta = par(3)
211 2240 : xsq = x * x
212 2240 : xdk = pow(x, 2.0d0 * dk + 1.0d0)
213 2240 : xdkp1 = xsq * xdk
214 2240 : dxst = sqrt(1.0d0 + 0.5d0 * xsq * theta)
215 :
216 : ! avoid an overflow in the denominator at large x:
217 2240 : if ((xsq-eta) < 1.d2) then
218 1760 : factor = exp(xsq - eta)
219 1760 : denom = factor + 1.0d0
220 1760 : fd = 2.0d0 * xdk * dxst/denom
221 1760 : fdeta = fd * factor/denom
222 1760 : denom2 = 4.0d0 * dxst * denom
223 1760 : fdtheta = 2.0d0 * xdkp1/denom2
224 :
225 : else
226 480 : factor = exp(eta - xsq)
227 480 : fd = 2.0d0 * xdk * dxst * factor
228 480 : fdeta = fd
229 480 : denom2 = 4.0d0 * dxst
230 480 : fdtheta = 2.0d0 * xdkp1/denom2 * factor
231 : end if
232 :
233 2240 : return
234 : end subroutine fdfunc2
235 :
236 : subroutine dqleg010(f,a,b,res,dres,ddres,par,n)
237 : !..
238 : !..10 point gauss-legendre rule for the fermi-dirac function and
239 : !..its derivatives with respect to eta and theta.
240 : !..on input f is the name of the subroutine containing the integrand,
241 : !..a is the lower end point of the interval, b is the higher end point,
242 : !..par is an array of constant parameters to be passed to subroutine f,
243 : !..and n is the length of the par array. on output res is the
244 : !..approximation from applying the 10-point gauss-legendre rule,
245 : !..dres is the derivative with respect to eta, and ddres is the
246 : !..derivative with respect to theta.
247 : !..
248 : !..note: since the number of nodes is even, zero is not an abscissa.
249 : !..
250 : !..declare
251 : interface
252 : subroutine f(absc1, par, n, fval1, dfval1, ddfval1)
253 : use const_def, only: dp
254 : implicit none
255 : integer :: n
256 : real(dp) :: absc1, par(n), fval1, dfval1, ddfval1
257 : end subroutine f
258 : end interface
259 : integer :: j, n
260 : real(dp) :: a,b,res,dres,ddres,par(n)
261 : real(dp) :: absc1,absc2,center,hlfrun,wg(5),xg(5)
262 : real(dp) :: fval1,dfval1,ddfval1,fval2,dfval2,ddfval2
263 :
264 : ! the abscissae and weights are given for the interval (-1,1).
265 : ! xg - abscissae of the 20-point gauss-legendre rule
266 : ! for half of the usual run (-1,1), i.e.
267 : ! the positive nodes of the 20-point rule
268 : ! wg - weights of the 20-point gauss rule.
269 : !
270 : ! abscissae and weights were evaluated with 100 decimal digit arithmetic.
271 :
272 : data xg ( 1) / 1.48874338981631210884826001129719984d-1 /
273 : data xg ( 2) / 4.33395394129247190799265943165784162d-1 /
274 : data xg ( 3) / 6.79409568299024406234327365114873575d-1 /
275 : data xg ( 4) / 8.65063366688984510732096688423493048d-1 /
276 : data xg ( 5) / 9.73906528517171720077964012084452053d-1 /
277 :
278 : data wg ( 1) / 2.95524224714752870173892994651338329d-1 /
279 : data wg ( 2) / 2.69266719309996355091226921569469352d-1 /
280 : data wg ( 3) / 2.19086362515982043995534934228163192d-1 /
281 : data wg ( 4) / 1.49451349150580593145776339657697332d-1 /
282 : data wg ( 5) / 6.66713443086881375935688098933317928d-2 /
283 :
284 :
285 : ! list of major variables
286 : ! -----------------------
287 : !
288 : ! absc - abscissa
289 : ! fval* - function value
290 : ! res - res of the 10-point gauss formula
291 :
292 : center = 0.5d0 * (a+b)
293 : hlfrun = 0.5d0 * (b-a)
294 : res = 0.0d0
295 : dres = 0.0d0
296 : ddres = 0.0d0
297 : do j=1,5
298 : absc1 = center + hlfrun*xg(j)
299 : absc2 = center - hlfrun*xg(j)
300 : call f(absc1, par, n, fval1, dfval1, ddfval1)
301 : call f(absc2, par, n, fval2, dfval2, ddfval2)
302 : res = res + (fval1 + fval2)*wg(j)
303 : dres = dres + (dfval1 + dfval2)*wg(j)
304 : ddres = ddres + (ddfval1 + ddfval2)*wg(j)
305 : end do
306 : res = res * hlfrun
307 : dres = dres * hlfrun
308 : ddres = ddres * hlfrun
309 : return
310 : end subroutine dqleg010
311 :
312 336 : subroutine dqleg020(f,a,b,res,dres,ddres,par,n)
313 : !..
314 : !..20 point gauss-legendre rule for the fermi-dirac function and
315 : !..its derivatives with respect to eta and theta.
316 : !..on input f is the name of the subroutine containing the integrand,
317 : !..a is the lower end point of the interval, b is the higher end point,
318 : !..par is an array of constant parameters to be passed to subroutine f,
319 : !..and n is the length of the par array. on output res is the
320 : !..approximation from applying the 20-point gauss-legendre rule,
321 : !..dres is the derivative with respect to eta, and ddres is the
322 : !..derivative with respect to theta.
323 : !..
324 : !..note: since the number of nodes is even, zero is not an abscissa.
325 : !..
326 : !..declare
327 : interface
328 : subroutine f(absc1, par, n, fval1, dfval1, ddfval1)
329 : use const_def, only: dp
330 : implicit none
331 : integer :: n
332 : real(dp) :: absc1, par(n), fval1, dfval1, ddfval1
333 : end subroutine f
334 : end interface
335 : integer :: j,n
336 : real(dp) :: a,b,res,dres,ddres,par(n)
337 336 : real(dp) :: absc1,absc2,center,hlfrun,wg(10),xg(10)
338 336 : real(dp) :: fval1,dfval1,ddfval1,fval2,dfval2,ddfval2
339 :
340 :
341 : ! the abscissae and weights are given for the interval (-1,1).
342 : ! xg - abscissae of the 20-point gauss-legendre rule
343 : ! for half of the usual run (-1,1), i.e.
344 : ! the positive nodes of the 20-point rule
345 : ! wg - weights of the 20-point gauss rule.
346 : !
347 : ! abscissae and weights were evaluated with 100 decimal digit arithmetic.
348 :
349 :
350 : data xg ( 1) / 7.65265211334973337546404093988382110d-2 /
351 : data xg ( 2) / 2.27785851141645078080496195368574624d-1 /
352 : data xg ( 3) / 3.73706088715419560672548177024927237d-1 /
353 : data xg ( 4) / 5.10867001950827098004364050955250998d-1 /
354 : data xg ( 5) / 6.36053680726515025452836696226285936d-1 /
355 : data xg ( 6) / 7.46331906460150792614305070355641590d-1 /
356 : data xg ( 7) / 8.39116971822218823394529061701520685d-1 /
357 : data xg ( 8) / 9.12234428251325905867752441203298113d-1 /
358 : data xg ( 9) / 9.63971927277913791267666131197277221d-1 /
359 : data xg ( 10) / 9.93128599185094924786122388471320278d-1 /
360 :
361 : data wg ( 1) / 1.52753387130725850698084331955097593d-1 /
362 : data wg ( 2) / 1.49172986472603746787828737001969436d-1 /
363 : data wg ( 3) / 1.42096109318382051329298325067164933d-1 /
364 : data wg ( 4) / 1.31688638449176626898494499748163134d-1 /
365 : data wg ( 5) / 1.18194531961518417312377377711382287d-1 /
366 : data wg ( 6) / 1.01930119817240435036750135480349876d-1 /
367 : data wg ( 7) / 8.32767415767047487247581432220462061d-2 /
368 : data wg ( 8) / 6.26720483341090635695065351870416063d-2 /
369 : data wg ( 9) / 4.06014298003869413310399522749321098d-2 /
370 : data wg ( 10) / 1.76140071391521183118619623518528163d-2 /
371 :
372 :
373 : ! list of major variables
374 : ! -----------------------
375 : !
376 : ! absc - abscissa
377 : ! fval* - function value
378 : ! res - res of the 20-point gauss formula
379 :
380 :
381 336 : center = 0.5d0 * (a+b)
382 336 : hlfrun = 0.5d0 * (b-a)
383 336 : res = 0.0d0
384 336 : dres = 0.0d0
385 336 : ddres = 0.0d0
386 3696 : do j=1,10
387 3360 : absc1 = center + hlfrun*xg(j)
388 3360 : absc2 = center - hlfrun*xg(j)
389 3360 : call f(absc1, par, n, fval1, dfval1, ddfval1)
390 3360 : call f(absc2, par, n, fval2, dfval2, ddfval2)
391 3360 : res = res + (fval1 + fval2)*wg(j)
392 3360 : dres = dres + (dfval1 + dfval2)*wg(j)
393 3696 : ddres = ddres + (ddfval1 + ddfval2)*wg(j)
394 : end do
395 336 : res = res * hlfrun
396 336 : dres = dres * hlfrun
397 336 : ddres = ddres * hlfrun
398 336 : return
399 : end subroutine dqleg020
400 :
401 : subroutine dqleg040(f,a,b,res,dres,ddres,par,n)
402 : !..
403 : !..40 point gauss-legendre rule for the fermi-dirac function and
404 : !..its derivatives with respect to eta and theta.
405 : !..on input f is the name of the subroutine containing the integrand,
406 : !..a is the lower end point of the interval, b is the higher end point,
407 : !..par is an array of constant parameters to be passed to subroutine f,
408 : !..and n is the length of the par array. on output res is the
409 : !..approximation from applying the 40-point gauss-legendre rule,
410 : !..dres is the derivative with respect to eta, and ddres is the
411 : !..derivative with respect to theta.
412 : !..
413 : !..note: since the number of nodes is even, zero is not an abscissa.
414 : !..
415 : !..declare
416 : interface
417 : subroutine f(absc1, par, n, fval1, dfval1, ddfval1)
418 : use const_def, only: dp
419 : implicit none
420 : integer :: n
421 : real(dp) :: absc1, par(n), fval1, dfval1, ddfval1
422 : end subroutine f
423 : end interface
424 : integer :: j,n
425 : real(dp) :: a,b,res,dres,ddres,par(n)
426 : real(dp) :: absc1,absc2,center,hlfrun,wg(20),xg(20)
427 : real(dp) :: fval1,dfval1,ddfval1,fval2,dfval2,ddfval2
428 :
429 :
430 : ! the abscissae and weights are given for the interval (-1,1).
431 : ! xg - abscissae of the 40-point gauss-legendre rule
432 : ! for half of the usual run (-1,1), i.e.
433 : ! the positive nodes of the 40-point rule
434 : ! wg - weights of the 40-point gauss rule.
435 : !
436 : ! abscissae and weights were evaluated with 100 decimal digit arithmetic.
437 :
438 : data xg ( 1) / 3.87724175060508219331934440246232946d-2 /
439 : data xg ( 2) / 1.16084070675255208483451284408024113d-1 /
440 : data xg ( 3) / 1.92697580701371099715516852065149894d-1 /
441 : data xg ( 4) / 2.68152185007253681141184344808596183d-1 /
442 : data xg ( 5) / 3.41994090825758473007492481179194310d-1 /
443 : data xg ( 6) / 4.13779204371605001524879745803713682d-1 /
444 : data xg ( 7) / 4.83075801686178712908566574244823004d-1 /
445 : data xg ( 8) / 5.49467125095128202075931305529517970d-1 /
446 : data xg ( 9) / 6.12553889667980237952612450230694877d-1 /
447 : data xg ( 10) / 6.71956684614179548379354514961494109d-1 /
448 : data xg ( 11) / 7.27318255189927103280996451754930548d-1 /
449 : data xg ( 12) / 7.78305651426519387694971545506494848d-1 /
450 : data xg ( 13) / 8.24612230833311663196320230666098773d-1 /
451 : data xg ( 14) / 8.65959503212259503820781808354619963d-1 /
452 : data xg ( 15) / 9.02098806968874296728253330868493103d-1 /
453 : data xg ( 16) / 9.32812808278676533360852166845205716d-1 /
454 : data xg ( 17) / 9.57916819213791655804540999452759285d-1 /
455 : data xg ( 18) / 9.77259949983774262663370283712903806d-1 /
456 : data xg ( 19) / 9.90726238699457006453054352221372154d-1 /
457 : data xg ( 20) / 9.98237709710559200349622702420586492d-1 /
458 :
459 : data wg ( 1) / 7.75059479784248112637239629583263269d-2 /
460 : data wg ( 2) / 7.70398181642479655883075342838102485d-2 /
461 : data wg ( 3) / 7.61103619006262423715580759224948230d-2 /
462 : data wg ( 4) / 7.47231690579682642001893362613246731d-2 /
463 : data wg ( 5) / 7.28865823958040590605106834425178358d-2 /
464 : data wg ( 6) / 7.06116473912867796954836308552868323d-2 /
465 : data wg ( 7) / 6.79120458152339038256901082319239859d-2 /
466 : data wg ( 8) / 6.48040134566010380745545295667527300d-2 /
467 : data wg ( 9) / 6.13062424929289391665379964083985959d-2 /
468 : data wg ( 10) / 5.74397690993915513666177309104259856d-2 /
469 : data wg ( 11) / 5.32278469839368243549964797722605045d-2 /
470 : data wg ( 12) / 4.86958076350722320614341604481463880d-2 /
471 : data wg ( 13) / 4.38709081856732719916746860417154958d-2 /
472 : data wg ( 14) / 3.87821679744720176399720312904461622d-2 /
473 : data wg ( 15) / 3.34601952825478473926781830864108489d-2 /
474 : data wg ( 16) / 2.79370069800234010984891575077210773d-2 /
475 : data wg ( 17) / 2.22458491941669572615043241842085732d-2 /
476 : data wg ( 18) / 1.64210583819078887128634848823639272d-2 /
477 : data wg ( 19) / 1.04982845311528136147421710672796523d-2 /
478 : data wg ( 20) / 4.52127709853319125847173287818533272d-3 /
479 :
480 :
481 : ! list of major variables
482 : ! -----------------------
483 : !
484 : ! absc - abscissa
485 : ! fval* - function value
486 : ! res - res of the 20-point gauss formula
487 :
488 :
489 : center = 0.5d0 * (a+b)
490 : hlfrun = 0.5d0 * (b-a)
491 : res = 0.0d0
492 : dres = 0.0d0
493 : ddres = 0.0d0
494 : do j=1,20
495 : absc1 = center + hlfrun*xg(j)
496 : absc2 = center - hlfrun*xg(j)
497 : call f(absc1, par, n, fval1, dfval1, ddfval1)
498 : call f(absc2, par, n, fval2, dfval2, ddfval2)
499 : res = res + (fval1 + fval2)*wg(j)
500 : dres = dres + (dfval1 + dfval2)*wg(j)
501 : ddres = ddres + (ddfval1 + ddfval2)*wg(j)
502 : end do
503 : res = res * hlfrun
504 : dres = dres * hlfrun
505 : ddres = ddres * hlfrun
506 : return
507 : end subroutine dqleg040
508 :
509 : subroutine dqleg080(f,a,b,res,dres,ddres,par,n)
510 : !..
511 : !..80 point gauss-legendre rule for the fermi-dirac function and
512 : !..its derivatives with respect to eta and theta.
513 : !..on input f is the name of the subroutine containing the integrand,
514 : !..a is the lower end point of the interval, b is the higher end point,
515 : !..par is an array of constant parameters to be passed to subroutine f,
516 : !..and n is the length of the par array. on output res is the
517 : !..approximation from applying the 80-point gauss-legendre rule,
518 : !..dres is the derivative with respect to eta, and ddres is the
519 : !..derivative with respect to theta.
520 : !..
521 : !..note: since the number of nodes is even, zero is not an abscissa.
522 : !..
523 : !..declare
524 : interface
525 : subroutine f(absc1, par, n, fval1, dfval1, ddfval1)
526 : use const_def, only: dp
527 : implicit none
528 : integer :: n
529 : real(dp) :: absc1, par(n), fval1, dfval1, ddfval1
530 : end subroutine f
531 : end interface
532 : integer :: j,n
533 : real(dp) :: a,b,res,dres,ddres,par(n)
534 : real(dp) :: absc1,absc2,center,hlfrun,wg(40),xg(40)
535 : real(dp) :: fval1,dfval1,ddfval1,fval2,dfval2,ddfval2
536 :
537 :
538 : ! the abscissae and weights are given for the interval (-1,1).
539 : ! xg - abscissae of the 80-point gauss-legendre rule
540 : ! for half of the usual run (-1,1), i.e.
541 : ! the positive nodes of the 80-point rule
542 : ! wg - weights of the 80-point gauss rule.
543 : !
544 : ! abscissae and weights were evaluated with 100 decimal digit arithmetic.
545 :
546 :
547 : data xg ( 1) / 1.95113832567939976543512341074545479d-2 /
548 : data xg ( 2) / 5.85044371524206686289933218834177944d-2 /
549 : data xg ( 3) / 9.74083984415845990632784501049369020d-2 /
550 : data xg ( 4) / 1.36164022809143886559241078000717067d-1 /
551 : data xg ( 5) / 1.74712291832646812559339048011286195d-1 /
552 : data xg ( 6) / 2.12994502857666132572388538666321823d-1 /
553 : data xg ( 7) / 2.50952358392272120493158816035004797d-1 /
554 : data xg ( 8) / 2.88528054884511853109139301434713898d-1 /
555 : data xg ( 9) / 3.25664370747701914619112943627358695d-1 /
556 : data xg ( 10) / 3.62304753499487315619043286358963588d-1 /
557 : data xg ( 11) / 3.98393405881969227024379642517533757d-1 /
558 : data xg ( 12) / 4.33875370831756093062386700363181958d-1 /
559 : data xg ( 13) / 4.68696615170544477036078364935808657d-1 /
560 : data xg ( 14) / 5.02804111888784987593672750367568003d-1 /
561 : data xg ( 15) / 5.36145920897131932019857253125400904d-1 /
562 : data xg ( 16) / 5.68671268122709784725485786624827158d-1 /
563 : data xg ( 17) / 6.00330622829751743154746299164006848d-1 /
564 : data xg ( 18) / 6.31075773046871966247928387289336863d-1 /
565 : data xg ( 19) / 6.60859898986119801735967122844317234d-1 /
566 : data xg ( 20) / 6.89637644342027600771207612438935266d-1 /
567 : data xg ( 21) / 7.17365185362099880254068258293815278d-1 /
568 : data xg ( 22) / 7.44000297583597272316540527930913673d-1 /
569 : data xg ( 23) / 7.69502420135041373865616068749026083d-1 /
570 : data xg ( 24) / 7.93832717504605449948639311738454358d-1 /
571 : data xg ( 25) / 8.16954138681463470371124994012295707d-1 /
572 : data xg ( 26) / 8.38831473580255275616623043902867064d-1 /
573 : data xg ( 27) / 8.59431406663111096977192123491656492d-1 /
574 : data xg ( 28) / 8.78722567678213828703773343639124407d-1 /
575 : data xg ( 29) / 8.96675579438770683194324071967395986d-1 /
576 : data xg ( 30) / 9.13263102571757654164733656150947478d-1 /
577 : data xg ( 31) / 9.28459877172445795953045959075453133d-1 /
578 : data xg ( 32) / 9.42242761309872674752266004500001735d-1 /
579 : data xg ( 33) / 9.54590766343634905493481517021029508d-1 /
580 : data xg ( 34) / 9.65485089043799251452273155671454998d-1 /
581 : data xg ( 35) / 9.74909140585727793385645230069136276d-1 /
582 : data xg ( 36) / 9.82848572738629070418288027709116473d-1 /
583 : data xg ( 37) / 9.89291302499755531026503167136631385d-1 /
584 : data xg ( 38) / 9.94227540965688277892063503664911698d-1 /
585 : data xg ( 39) / 9.97649864398237688899494208183122985d-1 /
586 : data xg ( 40) / 9.99553822651630629880080499094567184d-1 /
587 :
588 : data wg ( 1) / 3.90178136563066548112804392527540483d-2 /
589 : data wg ( 2) / 3.89583959627695311986255247722608223d-2 /
590 : data wg ( 3) / 3.88396510590519689317741826687871658d-2 /
591 : data wg ( 4) / 3.86617597740764633270771102671566912d-2 /
592 : data wg ( 5) / 3.84249930069594231852124363294901384d-2 /
593 : data wg ( 6) / 3.81297113144776383442067915657362019d-2 /
594 : data wg ( 7) / 3.77763643620013974897749764263210547d-2 /
595 : data wg ( 8) / 3.73654902387304900267053770578386691d-2 /
596 : data wg ( 9) / 3.68977146382760088391509965734052192d-2 /
597 : data wg ( 10) / 3.63737499058359780439649910465228136d-2 /
598 : data wg ( 11) / 3.57943939534160546028615888161544542d-2 /
599 : data wg ( 12) / 3.51605290447475934955265923886968812d-2 /
600 : data wg ( 13) / 3.44731204517539287943642267310298320d-2 /
601 : data wg ( 14) / 3.37332149846115228166751630642387284d-2 /
602 : data wg ( 15) / 3.29419393976454013828361809019595361d-2 /
603 : data wg ( 16) / 3.21004986734877731480564902872506960d-2 /
604 : data wg ( 17) / 3.12101741881147016424428667206035518d-2 /
605 : data wg ( 18) / 3.02723217595579806612200100909011747d-2 /
606 : data wg ( 19) / 2.92883695832678476927675860195791396d-2 /
607 : data wg ( 20) / 2.82598160572768623967531979650145302d-2 /
608 : data wg ( 21) / 2.71882275004863806744187066805442598d-2 /
609 : data wg ( 22) / 2.60752357675651179029687436002692871d-2 /
610 : data wg ( 23) / 2.49225357641154911051178470032198023d-2 /
611 : data wg ( 24) / 2.37318828659301012931925246135684162d-2 /
612 : data wg ( 25) / 2.25050902463324619262215896861687390d-2 /
613 : data wg ( 26) / 2.12440261157820063887107372506131285d-2 /
614 : data wg ( 27) / 1.99506108781419989288919287151135633d-2 /
615 : data wg ( 28) / 1.86268142082990314287354141521572090d-2 /
616 : data wg ( 29) / 1.72746520562693063585842071312909998d-2 /
617 : data wg ( 30) / 1.58961835837256880449029092291785257d-2 /
618 : data wg ( 31) / 1.44935080405090761169620745834605500d-2 /
619 : data wg ( 32) / 1.30687615924013392937868258970563403d-2 /
620 : data wg ( 33) / 1.16241141207978269164667699954326348d-2 /
621 : data wg ( 34) / 1.01617660411030645208318503524069436d-2 /
622 : data wg ( 35) / 8.68394526926085842640945220403428135d-3 /
623 : data wg ( 36) / 7.19290476811731275267557086795650747d-3 /
624 : data wg ( 37) / 5.69092245140319864926910711716201847d-3 /
625 : data wg ( 38) / 4.18031312469489523673930420168135132d-3 /
626 : data wg ( 39) / 2.66353358951268166929353583166845546d-3 /
627 : data wg ( 40) / 1.14495000318694153454417194131563611d-3 /
628 :
629 :
630 : ! list of major variables
631 : ! -----------------------
632 : !
633 : ! absc - abscissa
634 : ! fval* - function value
635 : ! res - res of the 20-point gauss formula
636 :
637 :
638 : center = 0.5d0 * (a+b)
639 : hlfrun = 0.5d0 * (b-a)
640 : res = 0.0d0
641 : dres = 0.0d0
642 : ddres = 0.0d0
643 : do j=1,40
644 : absc1 = center + hlfrun*xg(j)
645 : absc2 = center - hlfrun*xg(j)
646 : call f(absc1, par, n, fval1, dfval1, ddfval1)
647 : call f(absc2, par, n, fval2, dfval2, ddfval2)
648 : res = res + (fval1 + fval2)*wg(j)
649 : dres = dres + (dfval1 + dfval2)*wg(j)
650 : ddres = ddres + (ddfval1 + ddfval2)*wg(j)
651 : end do
652 : res = res * hlfrun
653 : dres = dres * hlfrun
654 : ddres = ddres * hlfrun
655 : return
656 : end subroutine dqleg080
657 :
658 : subroutine dqlag010(f,a,b,res,dres,ddres,par,n)
659 : !..
660 : !..10 point gauss-laguerre rule for the fermi-dirac function.
661 : !..on input f is the external function defining the integrand
662 : !..f(x)=g(x)*w(x), where w(x) is the gaussian weight
663 : !..w(x)=exp(-(x-a)/b) and g(x) a smooth function,
664 : !..a is the lower end point of the interval, b is the higher end point,
665 : !..par is an array of constant parameters to be passed to the function f,
666 : !..and n is the length of the par array. on output res is the
667 : !..approximation from applying the 10-point gauss-laguerre rule.
668 : !..since the number of nodes is even, zero is not an abscissa.
669 : !..
670 : !..declare
671 : interface
672 : subroutine f(absc, par, n, fval, dfval, ddfval)
673 : use const_def, only: dp
674 : implicit none
675 : integer :: n
676 : real(dp) :: absc, par(n), fval, dfval, ddfval
677 : end subroutine f
678 : end interface
679 : integer :: j,n
680 : real(dp) :: a,b,res,dres,ddres,par(n)
681 : real(dp) :: absc,wg(10),xg(10),fval,dfval,ddfval
682 :
683 :
684 : ! the abscissae and weights are given for the interval (0,+inf).
685 : ! xg - abscissae of the 10-point gauss-laguerre rule
686 : ! wg - weights of the 10-point gauss rule. since f yet
687 : ! includes the weight function, the values in wg
688 : ! are actually exp(xg) times the standard
689 : ! gauss-laguerre weights
690 : !
691 : ! abscissae and weights were evaluated with 100 decimal digit arithmetic.
692 :
693 : data xg ( 1) / 1.37793470540492430830772505652711188d-1 /
694 : data xg ( 2) / 7.29454549503170498160373121676078781d-1 /
695 : data xg ( 3) / 1.80834290174031604823292007575060883d0 /
696 : data xg ( 4) / 3.40143369785489951448253222140839067d0 /
697 : data xg ( 5) / 5.55249614006380363241755848686876285d0 /
698 : data xg ( 6) / 8.33015274676449670023876719727452218d0 /
699 : data xg ( 7) / 1.18437858379000655649185389191416139d1 /
700 : data xg ( 8) / 1.62792578313781020995326539358336223d1 /
701 : data xg ( 9) / 2.19965858119807619512770901955944939d1 /
702 : data xg ( 10) / 2.99206970122738915599087933407991951d1 /
703 :
704 : data wg ( 1) / 3.54009738606996308762226891442067608d-1 /
705 : data wg ( 2) / 8.31902301043580738109829658127849577d-1 /
706 : data wg ( 3) / 1.33028856174932817875279219439399369d0 /
707 : data wg ( 4) / 1.86306390311113098976398873548246693d0 /
708 : data wg ( 5) / 2.45025555808301016607269373165752256d0 /
709 : data wg ( 6) / 3.12276415513518249615081826331455472d0 /
710 : data wg ( 7) / 3.93415269556152109865581245924823077d0 /
711 : data wg ( 8) / 4.99241487219302310201148565243315445d0 /
712 : data wg ( 9) / 6.57220248513080297518766871037611234d0 /
713 : data wg ( 10) / 9.78469584037463069477008663871859813d0 /
714 :
715 :
716 : ! list of major variables
717 : ! -----------------------
718 : ! absc - abscissa
719 : ! fval* - function value
720 : ! res - res of the 10-point gauss formula
721 :
722 : res = 0.0d0
723 : dres = 0.0d0
724 : ddres = 0.0d0
725 : do j=1,10
726 : absc = a+b*xg(j)
727 : call f(absc, par, n, fval, dfval, ddfval)
728 : res = res + fval*wg(j)
729 : dres = dres + dfval*wg(j)
730 : ddres = ddres + ddfval*wg(j)
731 : end do
732 : res = res*b
733 : dres = dres*b
734 : ddres = ddres*b
735 : return
736 : end subroutine dqlag010
737 :
738 112 : subroutine dqlag020(f,a,b,res,dres,ddres,par,n)
739 : !..
740 : !..20 point gauss-laguerre rule for the fermi-dirac function.
741 : !..on input f is the external function defining the integrand
742 : !..f(x)=g(x)*w(x), where w(x) is the gaussian weight
743 : !..w(x)=dexp(-(x-a)/b) and g(x) a smooth function,
744 : !..a is the lower end point of the interval, b is the higher end point,
745 : !..par is an array of constant parameters to be passed to the function f,
746 : !..and n is the length of the par array. on output res is the
747 : !..approximation from applying the 20-point gauss-laguerre rule.
748 : !..since the number of nodes is even, zero is not an abscissa.
749 : !..
750 : !..declare
751 : interface
752 : subroutine f(absc, par, n, fval, dfval, ddfval)
753 : use const_def, only: dp
754 : implicit none
755 : integer :: n
756 : real(dp) :: absc, par(n), fval, dfval, ddfval
757 : end subroutine f
758 : end interface
759 : integer :: j,n
760 : real(dp) :: a,b,res,dres,ddres,par(n)
761 112 : real(dp) :: absc,wg(20),xg(20),fval,dfval,ddfval
762 :
763 :
764 : ! the abscissae and weights are given for the interval (0,+inf).
765 : ! xg - abscissae of the 20-point gauss-laguerre rule
766 : ! wg - weights of the 20-point gauss rule. since f yet
767 : ! includes the weight function, the values in wg
768 : ! are actually exp(xg) times the standard
769 : ! gauss-laguerre weights
770 : !
771 : ! abscissae and weights were evaluated with 100 decimal digit arithmetic.
772 :
773 : data xg ( 1) / 7.05398896919887533666890045842150958d-2 /
774 : data xg ( 2) / 3.72126818001611443794241388761146636d-1 /
775 : data xg ( 3) / 9.16582102483273564667716277074183187d-1 /
776 : data xg ( 4) / 1.70730653102834388068768966741305070d0 /
777 : data xg ( 5) / 2.74919925530943212964503046049481338d0 /
778 : data xg ( 6) / 4.04892531385088692237495336913333219d0 /
779 : data xg ( 7) / 5.61517497086161651410453988565189234d0 /
780 : data xg ( 8) / 7.45901745367106330976886021837181759d0 /
781 : data xg ( 9) / 9.59439286958109677247367273428279837d0 /
782 : data xg ( 10) / 1.20388025469643163096234092988655158d1 /
783 : data xg ( 11) / 1.48142934426307399785126797100479756d1 /
784 : data xg ( 12) / 1.79488955205193760173657909926125096d1 /
785 : data xg ( 13) / 2.14787882402850109757351703695946692d1 /
786 : data xg ( 14) / 2.54517027931869055035186774846415418d1 /
787 : data xg ( 15) / 2.99325546317006120067136561351658232d1 /
788 : data xg ( 16) / 3.50134342404790000062849359066881395d1 /
789 : data xg ( 17) / 4.08330570567285710620295677078075526d1 /
790 : data xg ( 18) / 4.76199940473465021399416271528511211d1 /
791 : data xg ( 19) / 5.58107957500638988907507734444972356d1 /
792 : data xg ( 20) / 6.65244165256157538186403187914606659d1 /
793 :
794 : data wg ( 1) / 1.81080062418989255451675405913110644d-1 /
795 : data wg ( 2) / 4.22556767878563974520344172566458197d-1 /
796 : data wg ( 3) / 6.66909546701848150373482114992515927d-1 /
797 : data wg ( 4) / 9.15352372783073672670604684771868067d-1 /
798 : data wg ( 5) / 1.16953970719554597380147822239577476d0 /
799 : data wg ( 6) / 1.43135498592820598636844994891514331d0 /
800 : data wg ( 7) / 1.70298113798502272402533261633206720d0 /
801 : data wg ( 8) / 1.98701589079274721410921839275129020d0 /
802 : data wg ( 9) / 2.28663578125343078546222854681495651d0 /
803 : data wg ( 10) / 2.60583472755383333269498950954033323d0 /
804 : data wg ( 11) / 2.94978373421395086600235416827285951d0 /
805 : data wg ( 12) / 3.32539578200931955236951937421751118d0 /
806 : data wg ( 13) / 3.74225547058981092111707293265377811d0 /
807 : data wg ( 14) / 4.21423671025188041986808063782478746d0 /
808 : data wg ( 15) / 4.76251846149020929695292197839096371d0 /
809 : data wg ( 16) / 5.42172604424557430380308297989981779d0 /
810 : data wg ( 17) / 6.25401235693242129289518490300707542d0 /
811 : data wg ( 18) / 7.38731438905443455194030019196464791d0 /
812 : data wg ( 19) / 9.15132873098747960794348242552950528d0 /
813 : data wg ( 20) / 1.28933886459399966710262871287485278d1 /
814 :
815 :
816 : ! list of major variables
817 : ! -----------------------
818 : ! absc - abscissa
819 : ! fval* - function value
820 : ! res - res of the 20-point gauss formula
821 :
822 112 : res = 0.0d0
823 112 : dres = 0.0d0
824 112 : ddres = 0.0d0
825 2352 : do j=1,20
826 2240 : absc = a+b*xg(j)
827 2240 : call f(absc, par, n, fval, dfval, ddfval)
828 2240 : res = res + fval*wg(j)
829 2240 : dres = dres + dfval*wg(j)
830 2352 : ddres = ddres + ddfval*wg(j)
831 : end do
832 112 : res = res*b
833 112 : dres = dres*b
834 112 : ddres = ddres*b
835 112 : return
836 : end subroutine dqlag020
837 :
838 : subroutine dqlag040(f,a,b,res,dres,ddres,par,n)
839 : !..
840 : !..20 point gauss-laguerre rule for the fermi-dirac function.
841 : !..on input f is the external function defining the integrand
842 : !..f(x)=g(x)*w(x), where w(x) is the gaussian weight
843 : !..w(x)=dexp(-(x-a)/b) and g(x) a smooth function,
844 : !..a is the lower end point of the interval, b is the higher end point,
845 : !..par is an array of constant parameters to be passed to the function f,
846 : !..and n is the length of the par array. on output res is the
847 : !..approximation from applying the 20-point gauss-laguerre rule.
848 : !..since the number of nodes is even, zero is not an abscissa.
849 : !..
850 : !..declare
851 : interface
852 : subroutine f(absc, par, n, fval, dfval, ddfval)
853 : use const_def, only: dp
854 : implicit none
855 : integer :: n
856 : real(dp) :: absc, par(n), fval, dfval, ddfval
857 : end subroutine f
858 : end interface
859 : integer :: j,n
860 : real(dp) :: a,b,res,dres,ddres,par(n)
861 : real(dp) :: absc,wg(40),xg(40),fval,dfval,ddfval
862 :
863 :
864 : ! the abscissae and weights are given for the interval (0,+inf).
865 : ! xg - abscissae of the 20-point gauss-laguerre rule
866 : ! wg - weights of the 20-point gauss rule. since f yet
867 : ! includes the weight function, the values in wg
868 : ! are actually exp(xg) times the standard
869 : ! gauss-laguerre weights
870 : !
871 : ! abscissae and weights were evaluated with 100 decimal digit arithmetic.
872 :
873 : data xg ( 1) / 3.57003943088883851220844712866008554d-2 /
874 : data xg ( 2) / 1.88162283158698516003589346219095913d-1 /
875 : data xg ( 3) / 4.62694281314576453564937524561190364d-1 /
876 : data xg ( 4) / 8.59772963972934922257272224688722412d-1 /
877 : data xg ( 5) / 1.38001082052733718649800032959526559d0 /
878 : data xg ( 6) / 2.02420913592282673344206600280013075d0 /
879 : data xg ( 7) / 2.79336935350681645765351448602664039d0 /
880 : data xg ( 8) / 3.68870267790827020959152635190868698d0 /
881 : data xg ( 9) / 4.71164114655497269361872283627747369d0 /
882 : data xg ( 10) / 5.86385087834371811427316423799582987d0 /
883 : data xg ( 11) / 7.14724790810228825068569195197942362d0 /
884 : data xg ( 12) / 8.56401701758616376271852204208813232d0 /
885 : data xg ( 13) / 1.01166340484519394068496296563952448d1 /
886 : data xg ( 14) / 1.18078922940045848428415867043606304d1 /
887 : data xg ( 15) / 1.36409337125370872283716763606501202d1 /
888 : data xg ( 16) / 1.56192858933390738372019636521880145d1 /
889 : data xg ( 17) / 1.77469059500956630425738774954243772d1 /
890 : data xg ( 18) / 2.00282328345748905296126148101751172d1 /
891 : data xg ( 19) / 2.24682499834984183513717862289945366d1 /
892 : data xg ( 20) / 2.50725607724262037943960862094009769d1 /
893 : data xg ( 21) / 2.78474800091688627207517041404557997d1 /
894 : data xg ( 22) / 3.08001457394454627007543851961911114d1 /
895 : data xg ( 23) / 3.39386570849137196090988585862819990d1 /
896 : data xg ( 24) / 3.72722458804760043283207609906074207d1 /
897 : data xg ( 25) / 4.08114928238869204661556755816006426d1 /
898 : data xg ( 26) / 4.45686031753344627071230206344983559d1 /
899 : data xg ( 27) / 4.85577635330599922809620488067067936d1 /
900 : data xg ( 28) / 5.27956111872169329693520211373917638d1 /
901 : data xg ( 29) / 5.73018633233936274950337469958921651d1 /
902 : data xg ( 30) / 6.21001790727751116121681990578989921d1 /
903 : data xg ( 31) / 6.72193709271269987990802775518887054d1 /
904 : data xg ( 32) / 7.26951588476124621175219277242619385d1 /
905 : data xg ( 33) / 7.85728029115713092805438968334812596d1 /
906 : data xg ( 34) / 8.49112311357049845427015647096663186d1 /
907 : data xg ( 35) / 9.17898746712363769923371934806273153d1 /
908 : data xg ( 36) / 9.93208087174468082501090541654868123d1 /
909 : data xg ( 37) / 1.07672440639388272520796767611322664d2 /
910 : data xg ( 38) / 1.17122309512690688807650644123550702d2 /
911 : data xg ( 39) / 1.28201841988255651192541104389631263d2 /
912 : data xg ( 40) / 1.42280044469159997888348835359541764d2 /
913 :
914 : data wg ( 1) / 9.16254711574598973115116980801374830d-2 /
915 : data wg ( 2) / 2.13420584905012080007193367121512341d-1 /
916 : data wg ( 3) / 3.35718116680284673880510701616292191d-1 /
917 : data wg ( 4) / 4.58540935033497560385432380376452497d-1 /
918 : data wg ( 5) / 5.82068165779105168990996365401543283d-1 /
919 : data wg ( 6) / 7.06495216367219392989830015673016682d-1 /
920 : data wg ( 7) / 8.32026903003485238099112947978349523d-1 /
921 : data wg ( 8) / 9.58878198794443111448122679676028906d-1 /
922 : data wg ( 9) / 1.08727616203054971575386933317202661d0 /
923 : data wg ( 10) / 1.21746232797778097895427785066560948d0 /
924 : data wg ( 11) / 1.34969549135676530792393859442394519d0 /
925 : data wg ( 12) / 1.48425492977684671120561178612978719d0 /
926 : data wg ( 13) / 1.62144416281182197802316884316454527d0 /
927 : data wg ( 14) / 1.76159537467676961118424220420981598d0 /
928 : data wg ( 15) / 1.90507466589479967668299320597279371d0 /
929 : data wg ( 16) / 2.05228834726171671760199582272947454d0 /
930 : data wg ( 17) / 2.20369055324509588909828344328140570d0 /
931 : data wg ( 18) / 2.35979253852320332354037375378901497d0 /
932 : data wg ( 19) / 2.52117414037643299165313690287422820d0 /
933 : data wg ( 20) / 2.68849805540884226415950544706374659d0 /
934 : data wg ( 21) / 2.86252781321044881203476395983104311d0 /
935 : data wg ( 22) / 3.04415066531151710041043967954333670d0 /
936 : data wg ( 23) / 3.23440709726353194177490239428867111d0 /
937 : data wg ( 24) / 3.43452939842774809220398481891602464d0 /
938 : data wg ( 25) / 3.64599282499408907238965646699490434d0 /
939 : data wg ( 26) / 3.87058459721651656808475320213444338d0 /
940 : data wg ( 27) / 4.11049868043282265583582247263951577d0 /
941 : data wg ( 28) / 4.36846872325406347450808338272945025d0 /
942 : data wg ( 29) / 4.64795898407446688299303399883883991d0 /
943 : data wg ( 30) / 4.95344611240989326218696150785562721d0 /
944 : data wg ( 31) / 5.29084840590073657468737365718858968d0 /
945 : data wg ( 32) / 5.66820460903297677000730529023263795d0 /
946 : data wg ( 33) / 6.09679641474342030593376010859198806d0 /
947 : data wg ( 34) / 6.59310886103999953794429664206294899d0 /
948 : data wg ( 35) / 7.18249599553689315064429801626699574d0 /
949 : data wg ( 36) / 7.90666631138422877369310742310586595d0 /
950 : data wg ( 37) / 8.84089249281034652079125595063026792d0 /
951 : data wg ( 38) / 1.01408992656211694839094600306940468d1 /
952 : data wg ( 39) / 1.22100212992046038985226485875881108d1 /
953 : data wg ( 40) / 1.67055206420242974052468774398573553d1 /
954 :
955 :
956 : ! list of major variables
957 : ! -----------------------
958 : ! absc - abscissa
959 : ! fval* - function value
960 : ! res - res of the 20-point gauss formula
961 :
962 :
963 : res = 0.0d0
964 : dres = 0.0d0
965 : ddres = 0.0d0
966 : do j=1,40
967 : absc = a+b*xg(j)
968 : call f(absc, par, n, fval, dfval, ddfval)
969 : res = res + fval*wg(j)
970 : dres = dres + dfval*wg(j)
971 : ddres = ddres + ddfval*wg(j)
972 : end do
973 : res = res*b
974 : dres = dres*b
975 : ddres = ddres*b
976 : return
977 : end subroutine dqlag040
978 :
979 : subroutine dqlag080(f,a,b,res,dres,ddres,par,n)
980 : !..
981 : !..20 point gauss-laguerre rule for the fermi-dirac function.
982 : !..on input f is the external function defining the integrand
983 : !..f(x)=g(x)*w(x), where w(x) is the gaussian weight
984 : !..w(x)=dexp(-(x-a)/b) and g(x) a smooth function,
985 : !..a is the lower end point of the interval, b is the higher end point,
986 : !..par is an array of constant parameters to be passed to the function f,
987 : !..and n is the length of the par array. on output res is the
988 : !..approximation from applying the 20-point gauss-laguerre rule.
989 : !..since the number of nodes is even, zero is not an abscissa.
990 : !..
991 : !..declare
992 : interface
993 : subroutine f(absc, par, n, fval, dfval, ddfval)
994 : use const_def, only: dp
995 : implicit none
996 : integer :: n
997 : real(dp) :: absc, par(n), fval, dfval, ddfval
998 : end subroutine f
999 : end interface
1000 : integer :: j,n
1001 : real(dp) :: a,b,res,dres,ddres,par(n)
1002 : real(dp) :: absc,wg(80),xg(80),fval,dfval,ddfval
1003 :
1004 :
1005 : ! the abscissae and weights are given for the interval (0,+inf).
1006 : ! xg - abscissae of the 20-point gauss-laguerre rule
1007 : ! wg - weights of the 20-point gauss rule. since f yet
1008 : ! includes the weight function, the values in wg
1009 : ! are actually exp(xg) times the standard
1010 : ! gauss-laguerre weights
1011 : !
1012 : ! abscissae and weights were evaluated with 100 decimal digit arithmetic.
1013 :
1014 : data xg ( 1) / 1.79604233006983655540103192474016803d-2 /
1015 : data xg ( 2) / 9.46399129943539888113902724652172943d-2 /
1016 : data xg ( 3) / 2.32622868125867569207706157216349831d-1 /
1017 : data xg ( 4) / 4.31992547802387480255786172497770411d-1 /
1018 : data xg ( 5) / 6.92828861352021839905702213635446867d-1 /
1019 : data xg ( 6) / 1.01523255618947143744625436859935350d0 /
1020 : data xg ( 7) / 1.39932768784287277414419051430978382d0 /
1021 : data xg ( 8) / 1.84526230383584513811177117769599966d0 /
1022 : data xg ( 9) / 2.35320887160926152447244708016140181d0 /
1023 : data xg ( 10) / 2.92336468655542632483691234259732862d0 /
1024 : data xg ( 11) / 3.55595231404613405944967308324638370d0 /
1025 : data xg ( 12) / 4.25122008230987808316485766448577637d0 /
1026 : data xg ( 13) / 5.00944263362016477243367706818206389d0 /
1027 : data xg ( 14) / 5.83092153860871901982127113295605083d0 /
1028 : data xg ( 15) / 6.71598597785131711156550087635199430d0 /
1029 : data xg ( 16) / 7.66499349489177306073418909047823480d0 /
1030 : data xg ( 17) / 8.67833082516770109543442255542661083d0 /
1031 : data xg ( 18) / 9.75641480574293071316550944617366591d0 /
1032 : data xg ( 19) / 1.08996933712878553774361001021489406d1 /
1033 : data xg ( 20) / 1.21086466423656999007054848698315593d1 /
1034 : data xg ( 21) / 1.33837881127786473701629840603833297d1 /
1035 : data xg ( 22) / 1.47256659435085855393358076261838437d1 /
1036 : data xg ( 23) / 1.61348643716624665791658545428990907d1 /
1037 : data xg ( 24) / 1.76120052438144378598635686943586520d1 /
1038 : data xg ( 25) / 1.91577496842412479221729970205674985d1 /
1039 : data xg ( 26) / 2.07727999097920960924419379010489579d1 /
1040 : data xg ( 27) / 2.24579012045404583114095916950877516d1 /
1041 : data xg ( 28) / 2.42138440689586473771922469392447092d1 /
1042 : data xg ( 29) / 2.60414665601655866929390053565435682d1 /
1043 : data xg ( 30) / 2.79416568418594655558233069293692111d1 /
1044 : data xg ( 31) / 2.99153559649009855011270412115737715d1 /
1045 : data xg ( 32) / 3.19635609022089207107748887542636533d1 /
1046 : data xg ( 33) / 3.40873278647261898749834947342860505d1 /
1047 : data xg ( 34) / 3.62877759287814544588031988436216948d1 /
1048 : data xg ( 35) / 3.85660910092922104582563052172908535d1 /
1049 : data xg ( 36) / 4.09235302180312671999095850595544326d1 /
1050 : data xg ( 37) / 4.33614266517312302957826760468219500d1 /
1051 : data xg ( 38) / 4.58811946612788863456266489974878378d1 /
1052 : data xg ( 39) / 4.84843356608331891358737273353563006d1 /
1053 : data xg ( 40) / 5.11724445446070105959889432334907144d1 /
1054 : data xg ( 41) / 5.39472167895544471206210278787572430d1 /
1055 : data xg ( 42) / 5.68104563346362231341248503244102122d1 /
1056 : data xg ( 43) / 5.97640843421099549427295961277471927d1 /
1057 : data xg ( 44) / 6.28101489639264772036272917590288682d1 /
1058 : data xg ( 45) / 6.59508362574560573434640627160792248d1 /
1059 : data xg ( 46) / 6.91884824202362773741980288648237373d1 /
1060 : data xg ( 47) / 7.25255875442633453588389652616568450d1 /
1061 : data xg ( 48) / 7.59648311278641748269449794974796502d1 /
1062 : data xg ( 49) / 7.95090896290888369620572826259980809d1 /
1063 : data xg ( 50) / 8.31614564010536896630429506875848705d1 /
1064 : data xg ( 51) / 8.69252644196156234481165926040448396d1 /
1065 : data xg ( 52) / 9.08041123009407559518411727820318427d1 /
1066 : data xg ( 53) / 9.48018942159474332072071889138735302d1 /
1067 : data xg ( 54) / 9.89228344469405791648019372738036790d1 /
1068 : data xg ( 55) / 1.03171527508039130233047094167345654d2 /
1069 : data xg ( 56) / 1.07552984977539906327607890798975954d2 /
1070 : data xg ( 57) / 1.12072690484128333623930046166211013d2 /
1071 : data xg ( 58) / 1.16736664673503666318157888130801099d2 /
1072 : data xg ( 59) / 1.21551542490952625566863895752110813d2 /
1073 : data xg ( 60) / 1.26524665796515540341570265431653573d2 /
1074 : data xg ( 61) / 1.31664195252120310870089086308006192d2 /
1075 : data xg ( 62) / 1.36979246686936973947570637289463788d2 /
1076 : data xg ( 63) / 1.42480058912161601930826569200455232d2 /
1077 : data xg ( 64) / 1.48178202455004441818652384836007732d2 /
1078 : data xg ( 65) / 1.54086842281798697859417425265596259d2 /
1079 : data xg ( 66) / 1.60221072870095715935268416893010646d2 /
1080 : data xg ( 67) / 1.66598351934053918744521179733712213d2 /
1081 : data xg ( 68) / 1.73239071334249503830906503775056999d2 /
1082 : data xg ( 69) / 1.80167323049032317982430208997701523d2 /
1083 : data xg ( 70) / 1.87411949676963772390490134588021771d2 /
1084 : data xg ( 71) / 1.95008022441532991450390479600599643d2 /
1085 : data xg ( 72) / 2.02998984195074937824807677823714777d2 /
1086 : data xg ( 73) / 2.11439870494836466691484904695542608d2 /
1087 : data xg ( 74) / 2.20402368151735739654044206677763168d2 /
1088 : data xg ( 75) / 2.29983206075680004348410969675844754d2 /
1089 : data xg ( 76) / 2.40319087055841540417597460479219628d2 /
1090 : data xg ( 77) / 2.51615879330499611167444939310973194d2 /
1091 : data xg ( 78) / 2.64213823883199102097696108691435553d2 /
1092 : data xg ( 79) / 2.78766733046004563652014172530611597d2 /
1093 : data xg ( 80) / 2.96966511995651345758852859155703581d2 /
1094 :
1095 : data wg ( 1) / 4.60931031330609664705251321395510083d-2 /
1096 : data wg ( 2) / 1.07313007783932752564150320304398860d-1 /
1097 : data wg ( 3) / 1.68664429547948111794220457782702406d-1 /
1098 : data wg ( 4) / 2.30088089384940054411257181978193282d-1 /
1099 : data wg ( 5) / 2.91601302502437964832169318772943752d-1 /
1100 : data wg ( 6) / 3.53226753575408236352723125805647046d-1 /
1101 : data wg ( 7) / 4.14988177550940466187197686311280092d-1 /
1102 : data wg ( 8) / 4.76909792302936241314777025418505661d-1 /
1103 : data wg ( 9) / 5.39016218474955374499507656522327912d-1 /
1104 : data wg ( 10) / 6.01332497447190529086765248840739512d-1 /
1105 : data wg ( 11) / 6.63884136396680571849442240727299214d-1 /
1106 : data wg ( 12) / 7.26697163614156688973567296249140514d-1 /
1107 : data wg ( 13) / 7.89798189428428531349793078398788294d-1 /
1108 : data wg ( 14) / 8.53214471438152298354598162431362968d-1 /
1109 : data wg ( 15) / 9.16973983833892698590342900031553302d-1 /
1110 : data wg ( 16) / 9.81105491004005747195060155984218607d-1 /
1111 : data wg ( 17) / 1.04563862580654218147568445663176029d0 /
1112 : data wg ( 18) / 1.11060397300025890771124763259729371d0 /
1113 : data wg ( 19) / 1.17603315841226175056651076519208666d0 /
1114 : data wg ( 20) / 1.24195894449809359279351761817871338d0 /
1115 : data wg ( 21) / 1.30841533303134064261188542845954645d0 /
1116 : data wg ( 22) / 1.37543767574892843813155917093490796d0 /
1117 : data wg ( 23) / 1.44306279387849270398312417207247308d0 /
1118 : data wg ( 24) / 1.51132910758830693847655020559917703d0 /
1119 : data wg ( 25) / 1.58027677653099415830201878723121659d0 /
1120 : data wg ( 26) / 1.64994785280267874116012042819355036d0 /
1121 : data wg ( 27) / 1.72038644781283277182004281452290770d0 /
1122 : data wg ( 28) / 1.79163891476093832891442620527688915d0 /
1123 : data wg ( 29) / 1.86375404864909708435925709028688162d0 /
1124 : data wg ( 30) / 1.93678330603070923513925434327841646d0 /
1125 : data wg ( 31) / 2.01078104701134222912614988175555546d0 /
1126 : data wg ( 32) / 2.08580480238741046429303978512989079d0 /
1127 : data wg ( 33) / 2.16191556924159897378316344048827763d0 /
1128 : data wg ( 34) / 2.23917813882364652373453997447445645d0 /
1129 : data wg ( 35) / 2.31766146114651854068606048043496370d0 /
1130 : data wg ( 36) / 2.39743905144001430514117238638849980d0 /
1131 : data wg ( 37) / 2.47858944444973417756369164455222527d0 /
1132 : data wg ( 38) / 2.56119670357790455335115509222572643d0 /
1133 : data wg ( 39) / 2.64535099306968892850463441000367534d0 /
1134 : data wg ( 40) / 2.73114922289915138861410287131169260d0 /
1135 : data wg ( 41) / 2.81869577775934171703141873747811157d0 /
1136 : data wg ( 42) / 2.90810334368223018934550276777492687d0 /
1137 : data wg ( 43) / 2.99949384839685626832412451829968724d0 /
1138 : data wg ( 44) / 3.09299953469357468116695108353033660d0 /
1139 : data wg ( 45) / 3.18876418994712376429365271501623466d0 /
1140 : data wg ( 46) / 3.28694455975337531998378107012216956d0 /
1141 : data wg ( 47) / 3.38771197960397652334054908762154571d0 /
1142 : data wg ( 48) / 3.49125426598732012281732423782764895d0 /
1143 : data wg ( 49) / 3.59777791769613046096294730174902943d0 /
1144 : data wg ( 50) / 3.70751069001745708341027155659228179d0 /
1145 : data wg ( 51) / 3.82070461965311695152029959430467622d0 /
1146 : data wg ( 52) / 3.93763959771430720676800540657330923d0 /
1147 : data wg ( 53) / 4.05862761338354481597420116187988679d0 /
1148 : data wg ( 54) / 4.18401782381424031850607692334503121d0 /
1149 : data wg ( 55) / 4.31420264929613425820084573217987912d0 /
1150 : data wg ( 56) / 4.44962515053655906604982820155377774d0 /
1151 : data wg ( 57) / 4.59078802263617511042959849148929810d0 /
1152 : data wg ( 58) / 4.73826464598929537394753873505838770d0 /
1153 : data wg ( 59) / 4.89271277966692168696886936743283567d0 /
1154 : data wg ( 60) / 5.05489168534039512820572507135175938d0 /
1155 : data wg ( 61) / 5.22568375594272391089278010166022467d0 /
1156 : data wg ( 62) / 5.40612213379727909853323512340717863d0 /
1157 : data wg ( 63) / 5.59742640184041404016553694158980053d0 /
1158 : data wg ( 64) / 5.80104932137643943530626162455394841d0 /
1159 : data wg ( 65) / 6.01873893878099108768015151514026344d0 /
1160 : data wg ( 66) / 6.25262247491437403092934213480091928d0 /
1161 : data wg ( 67) / 6.50532173517668675787482719663696133d0 /
1162 : data wg ( 68) / 6.78011521200777294201287347980059368d0 /
1163 : data wg ( 69) / 7.08117122025414518776174311916759402d0 /
1164 : data wg ( 70) / 7.41389244615305421421695606226687752d0 /
1165 : data wg ( 71) / 7.78544154841612700386232740339230532d0 /
1166 : data wg ( 72) / 8.20557347814596472333905086100917119d0 /
1167 : data wg ( 73) / 8.68801383996161871469419958058255237d0 /
1168 : data wg ( 74) / 9.25286973415578523923556506201979918d0 /
1169 : data wg ( 75) / 9.93114471840215736008370986534009772d0 /
1170 : data wg ( 76) / 1.07739736414646829405750843522990655d1 /
1171 : data wg ( 77) / 1.18738912465097447081950887710877400d1 /
1172 : data wg ( 78) / 1.34228858497264236139734940154089734d1 /
1173 : data wg ( 79) / 1.59197801616897924449554252200185978d1 /
1174 : data wg ( 80) / 2.14214542964372259537521036186415127d1 /
1175 :
1176 :
1177 : ! list of major variables
1178 : ! -----------------------
1179 : ! absc - abscissa
1180 : ! fval* - function value
1181 : ! res - res of the 20-point gauss formula
1182 :
1183 : res = 0.0d0
1184 : dres = 0.0d0
1185 : ddres = 0.0d0
1186 : do j=1,80
1187 : absc = a+b*xg(j)
1188 : call f(absc, par, n, fval, dfval, ddfval)
1189 : res = res + fval*wg(j)
1190 : dres = dres + dfval*wg(j)
1191 : ddres = ddres + ddfval*wg(j)
1192 : end do
1193 : res = res*b
1194 : dres = dres*b
1195 : ddres = ddres*b
1196 : return
1197 : end subroutine dqlag080
1198 :
1199 : end module gauss_fermi
|
