(git:374b731)
Loading...
Searching...
No Matches
beta_gamma_psi.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9MODULE beta_gamma_psi
10 ! not tested in the case where dp would stand for single precision
11 ! Routines for the beta function are untested
12
13 USE kinds, ONLY: dp
14#include "./base/base_uses.f90"
15
16 IMPLICIT NONE
17
18 PRIVATE
19
20 PUBLIC :: psi
21
22CONTAINS
23
24! **************************************************************************************************
25!> \brief ...
26!> \param i ...
27!> \return ...
28! **************************************************************************************************
29 FUNCTION ipmpar(i) RESULT(fn_val)
30!-----------------------------------------------------------------------
31
32! IPMPAR PROVIDES THE INTEGER MACHINE CONSTANTS FOR THE COMPUTER
33! THAT IS USED. IT IS ASSUMED THAT THE ARGUMENT I IS AN INTEGER
34! HAVING ONE OF THE VALUES 1-10. IPMPAR(I) HAS THE VALUE ...
35
36! INTEGERS.
37
38! ASSUME INTEGERS ARE REPRESENTED IN THE N-DIGIT, BASE-A FORM
39
40! SIGN ( X(N-1)*A**(N-1) + ... + X(1)*A + X(0) )
41
42! WHERE 0 .LE. X(I) .LT. A FOR I=0,...,N-1.
43
44! IPMPAR(1) = A, THE BASE (radix).
45
46! IPMPAR(2) = N, THE NUMBER OF BASE-A DIGITS (digits).
47
48! IPMPAR(3) = A**N - 1, THE LARGEST MAGNITUDE (huge).
49
50! FLOATING-POINT NUMBERS.
51
52! IT IS ASSUMED THAT THE SINGLE AND DOUBLE PRECISION FLOATING
53! POINT ARITHMETICS HAVE THE SAME BASE, SAY B, AND THAT THE
54! NONZERO NUMBERS ARE REPRESENTED IN THE FORM
55
56! SIGN (B**E) * (X(1)/B + ... + X(M)/B**M)
57
58! WHERE X(I) = 0,1,...,B-1 FOR I=1,...,M,
59! X(1) .GE. 1, AND EMIN .LE. E .LE. EMAX.
60
61! IPMPAR(4) = B, THE BASE.
62
63! SINGLE-PRECISION
64
65! IPMPAR(5) = M, THE NUMBER OF BASE-B DIGITS.
66
67! IPMPAR(6) = EMIN, THE SMALLEST EXPONENT E.
68
69! IPMPAR(7) = EMAX, THE LARGEST EXPONENT E.
70
71! DOUBLE-PRECISION
72
73! IPMPAR(8) = M, THE NUMBER OF BASE-B DIGITS.
74
75! IPMPAR(9) = EMIN, THE SMALLEST EXPONENT E.
76
77! IPMPAR(10) = EMAX, THE LARGEST EXPONENT E.
78
79!-----------------------------------------------------------------------
80
81 INTEGER, INTENT(IN) :: i
82 INTEGER :: fn_val
83
84 SELECT CASE (i)
85 CASE (1)
86 fn_val = radix(i)
87 CASE (2)
88 fn_val = digits(i)
89 CASE (3)
90 fn_val = huge(i)
91 CASE (4)
92 fn_val = radix(1.0)
93 CASE (5)
94 fn_val = digits(1.0)
95 CASE (6)
96 fn_val = minexponent(1.0)
97 CASE (7)
98 fn_val = maxexponent(1.0)
99 CASE (8)
100 fn_val = digits(1.0e0_dp)
101 CASE (9)
102 fn_val = minexponent(1.0e0_dp)
103 CASE (10)
104 fn_val = maxexponent(1.0e0_dp)
105 CASE DEFAULT
106 cpabort("unknown case")
107 END SELECT
108
109 END FUNCTION ipmpar
110
111! **************************************************************************************************
112!> \brief ...
113!> \param i ...
114!> \return ...
115! **************************************************************************************************
116 FUNCTION dpmpar(i) RESULT(fn_val)
117!-----------------------------------------------------------------------
118
119! DPMPAR PROVIDES THE DOUBLE PRECISION MACHINE CONSTANTS FOR
120! THE COMPUTER BEING USED. IT IS ASSUMED THAT THE ARGUMENT
121! I IS AN INTEGER HAVING ONE OF THE VALUES 1, 2, OR 3. IF THE
122! DOUBLE PRECISION ARITHMETIC BEING USED HAS M BASE B DIGITS AND
123! ITS SMALLEST AND LARGEST EXPONENTS ARE EMIN AND EMAX, THEN
124
125! DPMPAR(1) = B**(1 - M), THE MACHINE PRECISION,
126
127! DPMPAR(2) = B**(EMIN - 1), THE SMALLEST MAGNITUDE,
128
129! DPMPAR(3) = B**EMAX*(1 - B**(-M)), THE LARGEST MAGNITUDE.
130!-----------------------------------------------------------------------
131
132 INTEGER, INTENT(IN) :: i
133 REAL(dp) :: fn_val
134
135 REAL(dp), PARAMETER :: one = 1._dp
136
137! Local variable
138
139 SELECT CASE (i)
140 CASE (1)
141 fn_val = epsilon(one)
142 CASE (2)
143 fn_val = tiny(one)
144 CASE (3)
145 fn_val = huge(one)
146 END SELECT
147
148 END FUNCTION dpmpar
149
150! **************************************************************************************************
151!> \brief ...
152!> \param l ...
153!> \return ...
154! **************************************************************************************************
155 FUNCTION dxparg(l) RESULT(fn_val)
156!--------------------------------------------------------------------
157! IF L = 0 THEN DXPARG(L) = THE LARGEST POSITIVE W FOR WHICH
158! DEXP(W) CAN BE COMPUTED.
159!
160! IF L IS NONZERO THEN DXPARG(L) = THE LARGEST NEGATIVE W FOR
161! WHICH THE COMPUTED VALUE OF DEXP(W) IS NONZERO.
162!
163! NOTE... ONLY AN APPROXIMATE VALUE FOR DXPARG(L) IS NEEDED.
164!--------------------------------------------------------------------
165 INTEGER, INTENT(IN) :: l
166 REAL(dp) :: fn_val
167
168 REAL(dp), PARAMETER :: one = 1._dp
169
170! Local variable
171
172 IF (l == 0) THEN
173 fn_val = log(huge(one))
174 ELSE
175 fn_val = log(tiny(one))
176 END IF
177
178 END FUNCTION dxparg
179
180! **************************************************************************************************
181!> \brief ...
182!> \param a ...
183!> \return ...
184! **************************************************************************************************
185 FUNCTION alnrel(a) RESULT(fn_val)
186!-----------------------------------------------------------------------
187! EVALUATION OF THE FUNCTION LN(1 + A)
188!-----------------------------------------------------------------------
189 REAL(dp), INTENT(IN) :: a
190 REAL(dp) :: fn_val
191
192 REAL(dp), PARAMETER :: half = 0.5e0_dp, one = 1.e0_dp, p1 = -.129418923021993e+01_dp, &
193 p2 = .405303492862024e+00_dp, p3 = -.178874546012214e-01_dp, &
194 q1 = -.162752256355323e+01_dp, q2 = .747811014037616e+00_dp, &
195 q3 = -.845104217945565e-01_dp, two = 2.e0_dp, zero = 0.e0_dp
196
197 REAL(dp) :: t, t2, w, x
198
199!--------------------------
200
201 IF (abs(a) <= 0.375e0_dp) THEN
202 t = a/(a + two)
203 t2 = t*t
204 w = (((p3*t2 + p2)*t2 + p1)*t2 + one)/(((q3*t2 + q2)*t2 + q1)*t2 + one)
205 fn_val = two*t*w
206 ELSE
207 x = one + a
208 IF (a < zero) x = (a + half) + half
209 fn_val = log(x)
210 END IF
211
212 END FUNCTION alnrel
213
214! **************************************************************************************************
215!> \brief ...
216!> \param a ...
217!> \return ...
218! **************************************************************************************************
219 FUNCTION gam1(a) RESULT(fn_val)
220!-----------------------------------------------------------------------
221! COMPUTATION OF 1/GAMMA(A+1) - 1 FOR -0.5 <= A <= 1.5
222!-----------------------------------------------------------------------
223 REAL(dp), INTENT(IN) :: a
224 REAL(dp) :: fn_val
225
226 REAL(dp), PARAMETER :: p(7) = (/.577215664901533e+00_dp, -.409078193005776e+00_dp, &
227 -.230975380857675e+00_dp, .597275330452234e-01_dp, .766968181649490e-02_dp, &
228 -.514889771323592e-02_dp, .589597428611429e-03_dp/), q(5) = (/.100000000000000e+01_dp, &
229 .427569613095214e+00_dp, .158451672430138e+00_dp, .261132021441447e-01_dp, &
230 .423244297896961e-02_dp/), r(9) = (/-.422784335098468e+00_dp, -.771330383816272e+00_dp, &
231 -.244757765222226e+00_dp, .118378989872749e+00_dp, .930357293360349e-03_dp, &
232 -.118290993445146e-01_dp, .223047661158249e-02_dp, .266505979058923e-03_dp, &
233 -.132674909766242e-03_dp/)
234 REAL(dp), PARAMETER :: s1 = .273076135303957e+00_dp, s2 = .559398236957378e-01_dp
235
236 REAL(dp) :: bot, d, t, top, w
237
238 t = a
239 d = a - 0.5e0_dp
240 IF (d > 0.0e0_dp) t = d - 0.5e0_dp
241
242 IF (t > 0.e0_dp) THEN
243 top = (((((p(7)*t + p(6))*t + p(5))*t + p(4))*t + p(3))*t + p(2))*t + p(1)
244 bot = (((q(5)*t + q(4))*t + q(3))*t + q(2))*t + 1.0e0_dp
245 w = top/bot
246 IF (d > 0.0e0_dp) THEN
247 fn_val = (t/a)*((w - 0.5e0_dp) - 0.5e0_dp)
248 ELSE
249 fn_val = a*w
250 END IF
251 ELSE IF (t < 0.e0_dp) THEN
252 top = (((((((r(9)*t + r(8))*t + r(7))*t + r(6))*t + r(5))*t &
253 + r(4))*t + r(3))*t + r(2))*t + r(1)
254 bot = (s2*t + s1)*t + 1.0e0_dp
255 w = top/bot
256 IF (d > 0.0e0_dp) THEN
257 fn_val = t*w/a
258 ELSE
259 fn_val = a*((w + 0.5e0_dp) + 0.5e0_dp)
260 END IF
261 ELSE
262 fn_val = 0.0e0_dp
263 END IF
264
265 END FUNCTION gam1
266
267! **************************************************************************************************
268!> \brief ...
269!> \param a ...
270!> \param b ...
271!> \return ...
272! **************************************************************************************************
273 FUNCTION algdiv(a, b) RESULT(fn_val)
274!-----------------------------------------------------------------------
275
276! COMPUTATION OF LN(GAMMA(B)/GAMMA(A+B)) WHEN B >= 8
277
278! --------
279
280! IN THIS ALGORITHM, DEL(X) IS THE FUNCTION DEFINED BY
281! LN(GAMMA(X)) = (X - 0.5)*LN(X) - X + 0.5*LN(2*PI) + DEL(X).
282
283!-----------------------------------------------------------------------
284 REAL(dp), INTENT(IN) :: a, b
285 REAL(dp) :: fn_val
286
287 REAL(dp), PARAMETER :: c0 = .833333333333333e-01_dp, c1 = -.277777777760991e-02_dp, &
288 c2 = .793650666825390e-03_dp, c3 = -.595202931351870e-03_dp, &
289 c4 = .837308034031215e-03_dp, c5 = -.165322962780713e-02_dp
290
291 REAL(dp) :: c, d, h, s11, s3, s5, s7, s9, t, u, v, &
292 w, x, x2
293
294 IF (a > b) THEN
295 h = b/a
296 c = 1.0e0_dp/(1.0e0_dp + h)
297 x = h/(1.0e0_dp + h)
298 d = a + (b - 0.5e0_dp)
299 ELSE
300 h = a/b
301 c = h/(1.0e0_dp + h)
302 x = 1.0e0_dp/(1.0e0_dp + h)
303 d = b + (a - 0.5e0_dp)
304 END IF
305
306! SET SN = (1 - X**N)/(1 - X)
307
308 x2 = x*x
309 s3 = 1.0e0_dp + (x + x2)
310 s5 = 1.0e0_dp + (x + x2*s3)
311 s7 = 1.0e0_dp + (x + x2*s5)
312 s9 = 1.0e0_dp + (x + x2*s7)
313 s11 = 1.0e0_dp + (x + x2*s9)
314
315! SET W = DEL(B) - DEL(A + B)
316
317 t = (1.0e0_dp/b)**2
318 w = ((((c5*s11*t + c4*s9)*t + c3*s7)*t + c2*s5)*t + c1*s3)*t + c0
319 w = w*(c/b)
320
321! COMBINE THE RESULTS
322
323 u = d*alnrel(a/b)
324 v = a*(log(b) - 1.0e0_dp)
325 IF (u > v) THEN
326 fn_val = (w - v) - u
327 ELSE
328 fn_val = (w - u) - v
329 END IF
330
331 END FUNCTION algdiv
332
333! **************************************************************************************************
334!> \brief ...
335!> \param x ...
336!> \return ...
337! **************************************************************************************************
338 FUNCTION rexp(x) RESULT(fn_val)
339!-----------------------------------------------------------------------
340! EVALUATION OF THE FUNCTION EXP(X) - 1
341!-----------------------------------------------------------------------
342 REAL(dp), INTENT(IN) :: x
343 REAL(dp) :: fn_val
344
345 REAL(dp), PARAMETER :: p1 = .914041914819518e-09_dp, p2 = .238082361044469e-01_dp, &
346 q1 = -.499999999085958e+00_dp, q2 = .107141568980644e+00_dp, &
347 q3 = -.119041179760821e-01_dp, q4 = .595130811860248e-03_dp
348
349 REAL(dp) :: e
350
351 IF (abs(x) < 0.15e0_dp) THEN
352 fn_val = x*(((p2*x + p1)*x + 1.0e0_dp)/((((q4*x + q3)*x + q2)*x + q1)*x + 1.0e0_dp))
353 RETURN
354 END IF
355
356 IF (x < 0.0e0_dp) THEN
357 IF (x > -37.0e0_dp) THEN
358 fn_val = (exp(x) - 0.5e0_dp) - 0.5e0_dp
359 ELSE
360 fn_val = -1.0e0_dp
361 END IF
362 ELSE
363 e = exp(x)
364 fn_val = e*(0.5e0_dp + (0.5e0_dp - 1.0e0_dp/e))
365 END IF
366
367 END FUNCTION rexp
368
369! **************************************************************************************************
370!> \brief ...
371!> \param a ...
372!> \param b ...
373!> \param x ...
374!> \param y ...
375!> \param w ...
376!> \param eps ...
377!> \param ierr ...
378! **************************************************************************************************
379 SUBROUTINE bgrat(a, b, x, y, w, eps, ierr)
380!-----------------------------------------------------------------------
381! ASYMPTOTIC EXPANSION FOR IX(A,B) WHEN A IS LARGER THAN B.
382! THE RESULT OF THE EXPANSION IS ADDED TO W. IT IS ASSUMED
383! THAT A <= 15 AND B <= 1. EPS IS THE TOLERANCE USED.
384! IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS.
385!-----------------------------------------------------------------------
386 REAL(dp), INTENT(IN) :: a, b, x, y
387 REAL(dp), INTENT(INOUT) :: w
388 REAL(dp), INTENT(IN) :: eps
389 INTEGER, INTENT(OUT) :: ierr
390
391 REAL(dp), PARAMETER :: half = 0.5e0_dp, one = 1.e0_dp, &
392 quarter = 0.25e0_dp, zero = 0.e0_dp
393
394 INTEGER :: i, n
395 REAL(dp) :: bm1, bp2n, c(30), cn, coef, d(30), dj, &
396 j, l, lnx, n2, nu, q, r, s, sum, t, &
397 t2, tol, u, v, z
398
399 bm1 = (b - half) - half
400 nu = a + half*bm1
401 IF (y > 0.375e0_dp) THEN
402 lnx = log(x)
403 ELSE
404 lnx = alnrel(-y)
405 END IF
406 z = -nu*lnx
407 IF (b*z == zero) THEN
408 ! THE EXPANSION CANNOT BE COMPUTED
409 ierr = 1
410 RETURN
411 END IF
412
413! COMPUTATION OF THE EXPANSION
414! SET R = EXP(-Z)*Z**B/GAMMA(B)
415
416 r = b*(one + gam1(b))*exp(b*log(z))
417 r = r*exp(a*lnx)*exp(half*bm1*lnx)
418 u = algdiv(b, a) + b*log(nu)
419 u = r*exp(-u)
420 IF (u == zero) THEN
421 ! THE EXPANSION CANNOT BE COMPUTED
422 ierr = 1
423 RETURN
424 END IF
425 CALL grat1(b, z, r, q=q, eps=eps)
426
427 tol = 15.0e0_dp*eps
428 v = quarter*(one/nu)**2
429 t2 = quarter*lnx*lnx
430 l = w/u
431 j = q/r
432 sum = j
433 t = one
434 cn = one
435 n2 = zero
436 DO n = 1, 30
437 bp2n = b + n2
438 j = (bp2n*(bp2n + one)*j + (z + bp2n + one)*t)*v
439 n2 = n2 + 2.0e0_dp
440 t = t*t2
441 cn = cn/(n2*(n2 + one))
442 c(n) = cn
443 s = zero
444 IF (.NOT. (n == 1)) THEN
445 coef = b - n
446 DO i = 1, n - 1
447 s = s + coef*c(i)*d(n - i)
448 coef = coef + b
449 END DO
450 END IF
451 d(n) = bm1*cn + s/n
452 dj = d(n)*j
453 sum = sum + dj
454 IF (sum <= zero) THEN
455 ! THE EXPANSION CANNOT BE COMPUTED
456 ierr = 1
457 RETURN
458 END IF
459 IF (abs(dj) <= tol*(sum + l)) EXIT
460 END DO
461
462! ADD THE RESULTS TO W
463
464 ierr = 0
465 w = w + u*sum
466
467 END SUBROUTINE bgrat
468
469! **************************************************************************************************
470!> \brief ...
471!> \param a ...
472!> \param x ...
473!> \param r ...
474!> \param p ...
475!> \param q ...
476!> \param eps ...
477! **************************************************************************************************
478 SUBROUTINE grat1(a, x, r, p, q, eps)
479!-----------------------------------------------------------------------
480! EVALUATION OF P(A,X) AND Q(A,X) WHERE A <= 1 AND
481! THE INPUT ARGUMENT R HAS THE VALUE E**(-X)*X**A/GAMMA(A)
482!-----------------------------------------------------------------------
483 REAL(dp), INTENT(IN) :: a, x, r
484 REAL(dp), INTENT(OUT), OPTIONAL :: p, q
485 REAL(dp), INTENT(IN) :: eps
486
487 REAL(dp), PARAMETER :: half = 0.5e0_dp, one = 1.e0_dp, &
488 quarter = 0.25e0_dp, three = 3.e0_dp, &
489 two = 2.e0_dp, zero = 0.e0_dp
490
491 REAL(dp) :: a2n, a2nm1, an, b2n, b2nm1, c, g, h, j, &
492 l, pp, qq, sum, t, tol, w, z
493
494 IF (a*x == zero) THEN
495 IF (x <= a) THEN
496 pp = zero
497 qq = one
498 ELSE
499 pp = one
500 qq = zero
501 END IF
502 IF (PRESENT(p)) p = pp
503 IF (PRESENT(q)) q = qq
504 RETURN
505 END IF
506 IF (a == half) THEN
507 IF (x < quarter) THEN
508 pp = erf(sqrt(x))
509 qq = half + (half - pp)
510 ELSE
511 qq = erfc(sqrt(x))
512 pp = half + (half - qq)
513 END IF
514 IF (PRESENT(p)) p = pp
515 IF (PRESENT(q)) q = qq
516 RETURN
517 END IF
518 IF (x < 1.1e0_dp) THEN
519 ! TAYLOR SERIES FOR P(A,X)/X**A
520
521 an = three
522 c = x
523 sum = x/(a + three)
524 tol = three*eps/(a + one)
525 an = an + one
526 c = -c*(x/an)
527 t = c/(a + an)
528 sum = sum + t
529 DO WHILE (abs(t) > tol)
530 an = an + one
531 c = -c*(x/an)
532 t = c/(a + an)
533 sum = sum + t
534 END DO
535 j = a*x*((sum/6.0e0_dp - half/(a + two))*x + one/(a + one))
536
537 z = a*log(x)
538 h = gam1(a)
539 g = one + h
540 IF ((x < quarter .AND. z > -.13394e0_dp) .OR. a < x/2.59e0_dp) THEN
541 l = rexp(z)
542 qq = ((half + (half + l))*j - l)*g - h
543 IF (qq <= zero) THEN
544 pp = one
545 qq = zero
546 ELSE
547 pp = half + (half - qq)
548 END IF
549 ELSE
550 w = exp(z)
551 pp = w*g*(half + (half - j))
552 qq = half + (half - pp)
553 END IF
554 ELSE
555 ! CONTINUED FRACTION EXPANSION
556
557 tol = 8.0e0_dp*eps
558 a2nm1 = one
559 a2n = one
560 b2nm1 = x
561 b2n = x + (one - a)
562 c = one
563 DO
564 a2nm1 = x*a2n + c*a2nm1
565 b2nm1 = x*b2n + c*b2nm1
566 c = c + one
567 a2n = a2nm1 + (c - a)*a2n
568 b2n = b2nm1 + (c - a)*b2n
569 a2nm1 = a2nm1/b2n
570 b2nm1 = b2nm1/b2n
571 a2n = a2n/b2n
572 b2n = one
573 IF (abs(a2n - a2nm1/b2nm1) < tol*a2n) EXIT
574 END DO
575
576 qq = r*a2n
577 pp = half + (half - qq)
578 END IF
579
580 IF (PRESENT(p)) p = pp
581 IF (PRESENT(q)) q = qq
582
583 END SUBROUTINE grat1
584
585! **************************************************************************************************
586!> \brief ...
587!> \param mu ...
588!> \param x ...
589!> \return ...
590! **************************************************************************************************
591 FUNCTION esum(mu, x) RESULT(fn_val)
592!-----------------------------------------------------------------------
593! EVALUATION OF EXP(MU + X)
594!-----------------------------------------------------------------------
595 INTEGER, INTENT(IN) :: mu
596 REAL(dp), INTENT(IN) :: x
597 REAL(dp) :: fn_val
598
599 REAL(dp) :: w
600
601 IF (x > 0.0e0_dp) THEN
602 IF (mu > 0 .OR. mu + x < 0.0_dp) THEN
603 w = mu
604 fn_val = exp(w)*exp(x)
605 ELSE
606 w = mu + x
607 fn_val = exp(w)
608 END IF
609 ELSE
610 IF (mu < 0 .OR. mu + x < 0.0_dp) THEN
611 w = mu
612 fn_val = exp(w)*exp(x)
613 ELSE
614 w = mu + x
615 fn_val = exp(w)
616 END IF
617 END IF
618
619 END FUNCTION esum
620
621! **************************************************************************************************
622!> \brief ...
623!> \param x ...
624!> \return ...
625! **************************************************************************************************
626 FUNCTION rlog1(x) RESULT(fn_val)
627!-----------------------------------------------------------------------
628! EVALUATION OF THE FUNCTION X - LN(1 + X)
629!-----------------------------------------------------------------------
630! A = RLOG (0.7)
631! B = RLOG (4/3)
632!------------------------
633 REAL(dp), INTENT(IN) :: x
634 REAL(dp) :: fn_val
635
636 REAL(dp), PARAMETER :: a = .566749439387324e-01_dp, b = .456512608815524e-01_dp, &
637 p0 = .333333333333333e+00_dp, p1 = -.224696413112536e+00_dp, &
638 p2 = .620886815375787e-02_dp, q1 = -.127408923933623e+01_dp, q2 = .354508718369557e+00_dp
639
640 REAL(dp) :: r, t, u, up2, w, w1
641
642 IF (x < -0.39e0_dp .OR. x > 0.57e0_dp) THEN
643 w = (x + 0.5e0_dp) + 0.5e0_dp
644 fn_val = x - log(w)
645 RETURN
646 END IF
647
648! ARGUMENT REDUCTION
649
650 IF (x < -0.18e0_dp) THEN
651 u = (x + 0.3e0_dp)/0.7e0_dp
652 up2 = u + 2.0e0_dp
653 w1 = a - u*0.3e0_dp
654 ELSEIF (x > 0.18e0_dp) THEN
655 t = 0.75e0_dp*x
656 u = t - 0.25e0_dp
657 up2 = t + 1.75e0_dp
658 w1 = b + u/3.0e0_dp
659 ELSE
660 u = x
661 up2 = u + 2.0e0_dp
662 w1 = 0.0e0_dp
663 END IF
664
665! SERIES EXPANSION
666
667 r = u/up2
668 t = r*r
669 w = ((p2*t + p1)*t + p0)/((q2*t + q1)*t + 1.0e0_dp)
670 fn_val = r*(u - 2.0e0_dp*t*w) + w1
671
672 END FUNCTION rlog1
673
674! **************************************************************************************************
675!> \brief ...
676!> \param a ...
677!> \return ...
678! **************************************************************************************************
679 FUNCTION gamln1(a) RESULT(fn_val)
680!-----------------------------------------------------------------------
681! EVALUATION OF LN(GAMMA(1 + A)) FOR -0.2 .LE. A .LE. 1.25
682!-----------------------------------------------------------------------
683 REAL(dp), INTENT(IN) :: a
684 REAL(dp) :: fn_val
685
686 REAL(dp), PARAMETER :: p0 = .577215664901533e+00_dp, p1 = .844203922187225e+00_dp, &
687 p2 = -.168860593646662e+00_dp, p3 = -.780427615533591e+00_dp, &
688 p4 = -.402055799310489e+00_dp, p5 = -.673562214325671e-01_dp, &
689 p6 = -.271935708322958e-02_dp, q1 = .288743195473681e+01_dp, &
690 q2 = .312755088914843e+01_dp, q3 = .156875193295039e+01_dp, q4 = .361951990101499e+00_dp, &
691 q5 = .325038868253937e-01_dp, q6 = .667465618796164e-03_dp, r0 = .422784335098467e+00_dp, &
692 r1 = .848044614534529e+00_dp, r2 = .565221050691933e+00_dp, r3 = .156513060486551e+00_dp, &
693 r4 = .170502484022650e-01_dp, r5 = .497958207639485e-03_dp
694 REAL(dp), PARAMETER :: s1 = .124313399877507e+01_dp, s2 = .548042109832463e+00_dp, &
695 s3 = .101552187439830e+00_dp, s4 = .713309612391000e-02_dp, s5 = .116165475989616e-03_dp
696
697 REAL(dp) :: w, x
698
699 IF (a < 0.6e0_dp) THEN
700 w = ((((((p6*a + p5)*a + p4)*a + p3)*a + p2)*a + p1)*a + p0)/ &
701 ((((((q6*a + q5)*a + q4)*a + q3)*a + q2)*a + q1)*a + 1.0e0_dp)
702 fn_val = -a*w
703 ELSE
704 x = (a - 0.5e0_dp) - 0.5e0_dp
705 w = (((((r5*x + r4)*x + r3)*x + r2)*x + r1)*x + r0)/ &
706 (((((s5*x + s4)*x + s3)*x + s2)*x + s1)*x + 1.0e0_dp)
707 fn_val = x*w
708 END IF
709
710 END FUNCTION gamln1
711
712! **************************************************************************************************
713!> \brief ...
714!> \param xx ...
715!> \return ...
716! **************************************************************************************************
717 FUNCTION psi(xx) RESULT(fn_val)
718!---------------------------------------------------------------------
719
720! EVALUATION OF THE DIGAMMA FUNCTION
721
722! -----------
723
724! PSI(XX) IS ASSIGNED THE VALUE 0 WHEN THE DIGAMMA FUNCTION CANNOT
725! BE COMPUTED.
726
727! THE MAIN COMPUTATION INVOLVES EVALUATION OF RATIONAL CHEBYSHEV
728! APPROXIMATIONS PUBLISHED IN MATH. COMP. 27, 123-127(1973) BY
729! CODY, STRECOK AND THACHER.
730
731!---------------------------------------------------------------------
732! PSI WAS WRITTEN AT ARGONNE NATIONAL LABORATORY FOR THE FUNPACK
733! PACKAGE OF SPECIAL FUNCTION SUBROUTINES. PSI WAS MODIFIED BY
734! A.H. MORRIS (NSWC).
735!---------------------------------------------------------------------
736 REAL(dp), INTENT(IN) :: xx
737 REAL(dp) :: fn_val
738
739 REAL(dp), PARAMETER :: dx0 = 1.461632144968362341262659542325721325e0_dp, p1(7) = (/ &
740 .895385022981970e-02_dp, .477762828042627e+01_dp, .142441585084029e+03_dp, &
741 .118645200713425e+04_dp, .363351846806499e+04_dp, .413810161269013e+04_dp, &
742 .130560269827897e+04_dp/), p2(4) = (/-.212940445131011e+01_dp, -.701677227766759e+01_dp, &
743 -.448616543918019e+01_dp, -.648157123766197e+00_dp/), piov4 = .785398163397448e0_dp, q1(6)&
744 = (/.448452573429826e+02_dp, .520752771467162e+03_dp, .221000799247830e+04_dp, &
745 .364127349079381e+04_dp, .190831076596300e+04_dp, .691091682714533e-05_dp/)
746 REAL(dp), PARAMETER :: q2(4) = (/.322703493791143e+02_dp, .892920700481861e+02_dp, &
747 .546117738103215e+02_dp, .777788548522962e+01_dp/)
748
749 INTEGER :: i, m, n, nq
750 REAL(dp) :: aug, den, sgn, upper, w, x, xmax1, xmx0, &
751 xsmall, z
752
753!---------------------------------------------------------------------
754! PIOV4 = PI/4
755! DX0 = ZERO OF PSI TO EXTENDED PRECISION
756!---------------------------------------------------------------------
757!---------------------------------------------------------------------
758! COEFFICIENTS FOR RATIONAL APPROXIMATION OF
759! PSI(X) / (X - X0), 0.5 <= X <= 3.0
760!---------------------------------------------------------------------
761!---------------------------------------------------------------------
762! COEFFICIENTS FOR RATIONAL APPROXIMATION OF
763! PSI(X) - LN(X) + 1 / (2*X), X > 3.0
764!---------------------------------------------------------------------
765!---------------------------------------------------------------------
766! MACHINE DEPENDENT CONSTANTS ...
767! XMAX1 = THE SMALLEST POSITIVE FLOATING POINT CONSTANT
768! WITH ENTIRELY INTEGER REPRESENTATION. ALSO USED
769! AS NEGATIVE OF LOWER BOUND ON ACCEPTABLE NEGATIVE
770! ARGUMENTS AND AS THE POSITIVE ARGUMENT BEYOND WHICH
771! PSI MAY BE REPRESENTED AS ALOG(X).
772! XSMALL = ABSOLUTE ARGUMENT BELOW WHICH PI*COTAN(PI*X)
773! MAY BE REPRESENTED BY 1/X.
774!---------------------------------------------------------------------
775
776 xmax1 = ipmpar(3)
777 xmax1 = min(xmax1, 1.0e0_dp/dpmpar(1))
778 xsmall = 1.e-9_dp
779!---------------------------------------------------------------------
780 x = xx
781 aug = 0.0e0_dp
782 IF (x < 0.5e0_dp) THEN
783!---------------------------------------------------------------------
784! X .LT. 0.5, USE REFLECTION FORMULA
785! PSI(1-X) = PSI(X) + PI * COTAN(PI*X)
786!---------------------------------------------------------------------
787 IF (abs(x) <= xsmall) THEN
788 IF (x == 0.0e0_dp) THEN
789 ! ERROR RETURN
790 fn_val = 0.0e0_dp
791 RETURN
792 END IF
793!---------------------------------------------------------------------
794! 0 .LT. ABS(X) .LE. XSMALL. USE 1/X AS A SUBSTITUTE
795! FOR PI*COTAN(PI*X)
796!---------------------------------------------------------------------
797 aug = -1.0e0_dp/x
798 x = 1.0e0_dp - x
799 ELSE
800!---------------------------------------------------------------------
801! REDUCTION OF ARGUMENT FOR COTAN
802!---------------------------------------------------------------------
803 w = -x
804 sgn = piov4
805 IF (w <= 0.0e0_dp) THEN
806 w = -w
807 sgn = -sgn
808 END IF
809!---------------------------------------------------------------------
810! MAKE AN ERROR EXIT IF X .LE. -XMAX1
811!---------------------------------------------------------------------
812 IF (w >= xmax1) THEN
813 ! ERROR RETURN
814 fn_val = 0.0e0_dp
815 RETURN
816 END IF
817 nq = int(w)
818 w = w - nq
819 nq = int(w*4.0e0_dp)
820 w = 4.0e0_dp*(w - nq*.25e0_dp)
821!---------------------------------------------------------------------
822! W IS NOW RELATED TO THE FRACTIONAL PART OF 4.0 * X.
823! ADJUST ARGUMENT TO CORRESPOND TO VALUES IN FIRST
824! QUADRANT AND DETERMINE SIGN
825!---------------------------------------------------------------------
826 n = nq/2
827 IF ((n + n) /= nq) w = 1.0e0_dp - w
828 z = piov4*w
829 m = n/2
830 IF ((m + m) /= n) sgn = -sgn
831!---------------------------------------------------------------------
832! DETERMINE FINAL VALUE FOR -PI*COTAN(PI*X)
833!---------------------------------------------------------------------
834 n = (nq + 1)/2
835 m = n/2
836 m = m + m
837 IF (m /= n) THEN
838 aug = sgn*((sin(z)/cos(z))*4.0e0_dp)
839 ELSE
840 !---------------------------------------------------------------------
841 ! CHECK FOR SINGULARITY
842 !---------------------------------------------------------------------
843 IF (z == 0.0e0_dp) THEN
844 ! ERROR RETURN
845 fn_val = 0.0e0_dp
846 RETURN
847 END IF
848!---------------------------------------------------------------------
849! USE COS/SIN AS A SUBSTITUTE FOR COTAN, AND
850! SIN/COS AS A SUBSTITUTE FOR TAN
851!---------------------------------------------------------------------
852 aug = sgn*((cos(z)/sin(z))*4.0e0_dp)
853 END IF
854 x = 1.0e0_dp - x
855 END IF
856 END IF
857 IF (x <= 3.0e0_dp) THEN
858!---------------------------------------------------------------------
859! 0.5 .LE. X .LE. 3.0
860!---------------------------------------------------------------------
861 den = x
862 upper = p1(1)*x
863
864 DO i = 1, 5
865 den = (den + q1(i))*x
866 upper = (upper + p1(i + 1))*x
867 END DO
868
869 den = (upper + p1(7))/(den + q1(6))
870 xmx0 = x - dx0
871 fn_val = den*xmx0 + aug
872 RETURN
873 END IF
874!---------------------------------------------------------------------
875! IF X .GE. XMAX1, PSI = LN(X)
876!---------------------------------------------------------------------
877 IF (x < xmax1) THEN
878!---------------------------------------------------------------------
879! 3.0 .LT. X .LT. XMAX1
880!---------------------------------------------------------------------
881 w = 1.0e0_dp/(x*x)
882 den = w
883 upper = p2(1)*w
884
885 DO i = 1, 3
886 den = (den + q2(i))*w
887 upper = (upper + p2(i + 1))*w
888 END DO
889
890 aug = upper/(den + q2(4)) - 0.5e0_dp/x + aug
891 END IF
892 fn_val = aug + log(x)
893
894 END FUNCTION psi
895
896! **************************************************************************************************
897!> \brief ...
898!> \param a0 ...
899!> \param b0 ...
900!> \return ...
901! **************************************************************************************************
902 FUNCTION betaln(a0, b0) RESULT(fn_val)
903!-----------------------------------------------------------------------
904! EVALUATION OF THE LOGARITHM OF THE BETA FUNCTION
905!-----------------------------------------------------------------------
906! E = 0.5*LN(2*PI)
907!--------------------------
908 REAL(dp), INTENT(IN) :: a0, b0
909 REAL(dp) :: fn_val
910
911 REAL(dp), PARAMETER :: e = .918938533204673e0_dp
912
913 INTEGER :: i, n
914 REAL(dp) :: a, b, c, h, u, v, w, z
915
916!--------------------------
917
918 a = min(a0, b0)
919 b = max(a0, b0)
920!-----------------------------------------------------------------------
921! PROCEDURE WHEN A .GE. 8
922!-----------------------------------------------------------------------
923 IF (a >= 8.0e0_dp) THEN
924 w = bcorr(a, b)
925 h = a/b
926 c = h/(1.0e0_dp + h)
927 u = -(a - 0.5e0_dp)*log(c)
928 v = b*alnrel(h)
929 IF (u > v) THEN
930 fn_val = (((-0.5e0_dp*log(b) + e) + w) - v) - u
931 ELSE
932 fn_val = (((-0.5e0_dp*log(b) + e) + w) - u) - v
933 END IF
934!-----------------------------------------------------------------------
935! PROCEDURE WHEN A .LT. 1
936!-----------------------------------------------------------------------
937 ELSEIF (a < 1.0e0_dp) THEN
938 IF (b < 8.0e0_dp) THEN
939 fn_val = log_gamma(a) + (log_gamma(b) - log_gamma(a + b))
940 ELSE
941 fn_val = log_gamma(a) + algdiv(a, b)
942 END IF
943!-----------------------------------------------------------------------
944! PROCEDURE WHEN 1 .LE. A .LT. 8
945!-----------------------------------------------------------------------
946 ELSEIF (a <= 2.0e0_dp) THEN
947 IF (b <= 2.0e0_dp) THEN
948 fn_val = log_gamma(a) + log_gamma(b) - gsumln(a, b)
949 RETURN
950 END IF
951 w = 0.0e0_dp
952 IF (b < 8.0e0_dp) THEN
953 ! REDUCTION OF B WHEN B .LT. 8
954
955 n = int(b - 1.0e0_dp)
956 z = 1.0e0_dp
957 DO i = 1, n
958 b = b - 1.0e0_dp
959 z = z*(b/(a + b))
960 END DO
961 fn_val = w + log(z) + (log_gamma(a) + (log_gamma(b) - gsumln(a, b)))
962 RETURN
963 END IF
964 fn_val = log_gamma(a) + algdiv(a, b)
965! REDUCTION OF A WHEN B .LE. 1000
966 ELSEIF (b <= 1000.0e0_dp) THEN
967 n = int(a - 1.0e0_dp)
968 w = 1.0e0_dp
969 DO i = 1, n
970 a = a - 1.0e0_dp
971 h = a/b
972 w = w*(h/(1.0e0_dp + h))
973 END DO
974 w = log(w)
975 IF (b >= 8.0e0_dp) THEN
976 fn_val = w + log_gamma(a) + algdiv(a, b)
977 RETURN
978 END IF
979
980 ! REDUCTION OF B WHEN B .LT. 8
981
982 n = int(b - 1.0e0_dp)
983 z = 1.0e0_dp
984 DO i = 1, n
985 b = b - 1.0e0_dp
986 z = z*(b/(a + b))
987 END DO
988 fn_val = w + log(z) + (log_gamma(a) + (log_gamma(b) - gsumln(a, b)))
989 ELSE
990! REDUCTION OF A WHEN B .GT. 1000
991
992 n = int(a - 1.0e0_dp)
993 w = 1.0e0_dp
994 DO i = 1, n
995 a = a - 1.0e0_dp
996 w = w*(a/(1.0e0_dp + a/b))
997 END DO
998 fn_val = (log(w) - n*log(b)) + (log_gamma(a) + algdiv(a, b))
999 END IF
1000
1001 END FUNCTION betaln
1002
1003! **************************************************************************************************
1004!> \brief ...
1005!> \param a ...
1006!> \param b ...
1007!> \return ...
1008! **************************************************************************************************
1009 FUNCTION gsumln(a, b) RESULT(fn_val)
1010!-----------------------------------------------------------------------
1011! EVALUATION OF THE FUNCTION LN(GAMMA(A + B))
1012! FOR 1 .LE. A .LE. 2 AND 1 .LE. B .LE. 2
1013!-----------------------------------------------------------------------
1014 REAL(dp), INTENT(IN) :: a, b
1015 REAL(dp) :: fn_val
1016
1017 REAL(dp) :: x
1018
1019 x = a + b - 2.e0_dp
1020 IF (x <= 0.25e0_dp) THEN
1021 fn_val = gamln1(1.0e0_dp + x)
1022 ELSEIF (x <= 1.25e0_dp) THEN
1023 fn_val = gamln1(x) + alnrel(x)
1024 ELSE
1025 fn_val = gamln1(x - 1.0e0_dp) + log(x*(1.0e0_dp + x))
1026 END IF
1027 END FUNCTION gsumln
1028
1029! **************************************************************************************************
1030!> \brief ...
1031!> \param a0 ...
1032!> \param b0 ...
1033!> \return ...
1034! **************************************************************************************************
1035 FUNCTION bcorr(a0, b0) RESULT(fn_val)
1036!-----------------------------------------------------------------------
1037
1038! EVALUATION OF DEL(A0) + DEL(B0) - DEL(A0 + B0) WHERE
1039! LN(GAMMA(A)) = (A - 0.5)*LN(A) - A + 0.5*LN(2*PI) + DEL(A).
1040! IT IS ASSUMED THAT A0 .GE. 8 AND B0 .GE. 8.
1041
1042!-----------------------------------------------------------------------
1043 REAL(dp), INTENT(IN) :: a0, b0
1044 REAL(dp) :: fn_val
1045
1046 REAL(dp), PARAMETER :: c0 = .833333333333333e-01_dp, c1 = -.277777777760991e-02_dp, &
1047 c2 = .793650666825390e-03_dp, c3 = -.595202931351870e-03_dp, &
1048 c4 = .837308034031215e-03_dp, c5 = -.165322962780713e-02_dp
1049
1050 REAL(dp) :: a, b, c, h, s11, s3, s5, s7, s9, t, w, &
1051 x, x2
1052
1053 a = min(a0, b0)
1054 b = max(a0, b0)
1055
1056 h = a/b
1057 c = h/(1.0e0_dp + h)
1058 x = 1.0e0_dp/(1.0e0_dp + h)
1059 x2 = x*x
1060
1061! SET SN = (1 - X**N)/(1 - X)
1062
1063 s3 = 1.0e0_dp + (x + x2)
1064 s5 = 1.0e0_dp + (x + x2*s3)
1065 s7 = 1.0e0_dp + (x + x2*s5)
1066 s9 = 1.0e0_dp + (x + x2*s7)
1067 s11 = 1.0e0_dp + (x + x2*s9)
1068
1069! SET W = DEL(B) - DEL(A + B)
1070
1071 t = (1.0e0_dp/b)**2
1072 w = ((((c5*s11*t + c4*s9)*t + c3*s7)*t + c2*s5)*t + c1*s3)*t + c0
1073 w = w*(c/b)
1074
1075! COMPUTE DEL(A) + W
1076
1077 t = (1.0e0_dp/a)**2
1078 fn_val = (((((c5*t + c4)*t + c3)*t + c2)*t + c1)*t + c0)/a + w
1079 RETURN
1080 END FUNCTION bcorr
1081
1082! **************************************************************************************************
1083!> \brief ...
1084!> \param a ...
1085!> \param b ...
1086!> \param x ...
1087!> \param y ...
1088!> \param w ...
1089!> \param w1 ...
1090!> \param ierr ...
1091! **************************************************************************************************
1092 SUBROUTINE bratio(a, b, x, y, w, w1, ierr)
1093!-----------------------------------------------------------------------
1094
1095! EVALUATION OF THE INCOMPLETE BETA FUNCTION IX(A,B)
1096
1097! --------------------
1098
1099! IT IS ASSUMED THAT A AND B ARE NONNEGATIVE, AND THAT X <= 1
1100! AND Y = 1 - X. BRATIO ASSIGNS W AND W1 THE VALUES
1101
1102! W = IX(A,B)
1103! W1 = 1 - IX(A,B)
1104
1105! IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS.
1106! IF NO INPUT ERRORS ARE DETECTED THEN IERR IS SET TO 0 AND
1107! W AND W1 ARE COMPUTED. OTHERWISE, IF AN ERROR IS DETECTED,
1108! THEN W AND W1 ARE ASSIGNED THE VALUE 0 AND IERR IS SET TO
1109! ONE OF THE FOLLOWING VALUES ...
1110
1111! IERR = 1 IF A OR B IS NEGATIVE
1112! IERR = 2 IF A = B = 0
1113! IERR = 3 IF X .LT. 0 OR X .GT. 1
1114! IERR = 4 IF Y .LT. 0 OR Y .GT. 1
1115! IERR = 5 IF X + Y .NE. 1
1116! IERR = 6 IF X = A = 0
1117! IERR = 7 IF Y = B = 0
1118
1119!--------------------
1120! WRITTEN BY ALFRED H. MORRIS, JR.
1121! NAVAL SURFACE WARFARE CENTER
1122! DAHLGREN, VIRGINIA
1123! REVISED ... APRIL 1993
1124!-----------------------------------------------------------------------
1125 REAL(dp), INTENT(IN) :: a, b, x, y
1126 REAL(dp), INTENT(OUT) :: w, w1
1127 INTEGER, INTENT(OUT) :: ierr
1128
1129 INTEGER :: ierr1, ind, n
1130 REAL(dp) :: a0, b0, eps, lambda, t, x0, y0, z
1131
1132!-----------------------------------------------------------------------
1133! ****** EPS IS A MACHINE DEPENDENT CONSTANT. EPS IS THE SMALLEST
1134! FLOATING POINT NUMBER FOR WHICH 1.0 + EPS .GT. 1.0
1135
1136 eps = dpmpar(1)
1137
1138 w = 0.0e0_dp
1139 w1 = 0.0e0_dp
1140 IF (a < 0.0e0_dp .OR. b < 0.0e0_dp) THEN
1141 ierr = 1
1142 RETURN
1143 END IF
1144 IF (a == 0.0e0_dp .AND. b == 0.0e0_dp) THEN
1145 ierr = 2
1146 RETURN
1147 END IF
1148 IF (x < 0.0e0_dp .OR. x > 1.0e0_dp) THEN
1149 ierr = 3
1150 RETURN
1151 END IF
1152 IF (y < 0.0e0_dp .OR. y > 1.0e0_dp) THEN
1153 ierr = 4
1154 RETURN
1155 END IF
1156 z = ((x + y) - 0.5e0_dp) - 0.5e0_dp
1157 IF (abs(z) > 3.0e0_dp*eps) THEN
1158 ierr = 5
1159 RETURN
1160 END IF
1161
1162 ierr = 0
1163
1164 IF (x == 0.0e0_dp .OR. a == 0.0e0_dp) THEN
1165 IF (x == 0.0e0_dp .AND. a == 0.0e0_dp) THEN
1166 ierr = 6
1167 ELSE
1168 w = 0.0e0_dp
1169 w1 = 1.0e0_dp
1170 END IF
1171 RETURN
1172 END IF
1173
1174 IF (y == 0.0e0_dp .OR. b == 0.0e0_dp) THEN
1175 IF (y == 0.0e0_dp .AND. b == 0.0e0_dp) THEN
1176 ierr = 7
1177 ELSE
1178 w = 1.0e0_dp
1179 w1 = 0.0e0_dp
1180 END IF
1181 RETURN
1182 END IF
1183
1184 eps = max(eps, 1.e-15_dp)
1185 IF (max(a, b) < 1.e-3_dp*eps) THEN
1186! PROCEDURE FOR A AND B .LT. 1.E-3*EPS
1187 w = b/(a + b)
1188 w1 = a/(a + b)
1189 RETURN
1190 END IF
1191
1192 ind = 0
1193 a0 = a
1194 b0 = b
1195 x0 = x
1196 y0 = y
1197 IF (min(a0, b0) > 1.0e0_dp) GO TO 30
1198
1199! PROCEDURE FOR A0 .LE. 1 OR B0 .LE. 1
1200
1201 IF (x > 0.5e0_dp) THEN
1202 ind = 1
1203 a0 = b
1204 b0 = a
1205 x0 = y
1206 y0 = x
1207 END IF
1208
1209 IF (b0 < min(eps, eps*a0)) GO TO 80
1210 IF (a0 < min(eps, eps*b0) .AND. b0*x0 <= 1.0e0_dp) GO TO 90
1211 IF (max(a0, b0) > 1.0e0_dp) GO TO 20
1212 IF (a0 >= min(0.2e0_dp, b0)) GO TO 100
1213 IF (x0**a0 <= 0.9e0_dp) GO TO 100
1214 IF (x0 >= 0.3e0_dp) GO TO 110
1215 n = 20
1216 GO TO 130
1217
121820 IF (b0 <= 1.0e0_dp) GO TO 100
1219 IF (x0 >= 0.3e0_dp) GO TO 110
1220 IF (x0 < 0.1e0_dp .AND. (x0*b0)**a0 <= 0.7e0_dp) GO TO 100
1221 IF (b0 > 15.0e0_dp) GO TO 131
1222 n = 20
1223 GO TO 130
1224
1225! PROCEDURE FOR A0 .GT. 1 AND B0 .GT. 1
1226
122730 IF (a > b) THEN
1228 lambda = (a + b)*y - b
1229 ELSE
1230 lambda = a - (a + b)*x
1231 END IF
1232
1233 IF (lambda < 0.0e0_dp) THEN
1234 ind = 1
1235 a0 = b
1236 b0 = a
1237 x0 = y
1238 y0 = x
1239 lambda = abs(lambda)
1240 END IF
1241
1242 IF (b0 < 40.0e0_dp .AND. b0*x0 <= 0.7e0_dp) GO TO 100
1243 IF (b0 < 40.0e0_dp) GO TO 140
1244 IF (a0 > b0) GO TO 50
1245 IF (a0 <= 100.0e0_dp) GO TO 120
1246 IF (lambda > 0.03e0_dp*a0) GO TO 120
1247 GO TO 180
124850 IF (b0 <= 100.0e0_dp) GO TO 120
1249 IF (lambda > 0.03e0_dp*b0) GO TO 120
1250 GO TO 180
1251
1252! EVALUATION OF THE APPROPRIATE ALGORITHM
1253
125480 w = fpser(a0, b0, x0, eps)
1255 w1 = 0.5e0_dp + (0.5e0_dp - w)
1256 GO TO 220
1257
125890 w1 = apser(a0, b0, x0, eps)
1259 w = 0.5e0_dp + (0.5e0_dp - w1)
1260 GO TO 220
1261
1262100 w = bpser(a0, b0, x0, eps)
1263 w1 = 0.5e0_dp + (0.5e0_dp - w)
1264 GO TO 220
1265
1266110 w1 = bpser(b0, a0, y0, eps)
1267 w = 0.5e0_dp + (0.5e0_dp - w1)
1268 GO TO 220
1269
1270120 w = bfrac(a0, b0, x0, y0, lambda, 15.0e0_dp*eps)
1271 w1 = 0.5e0_dp + (0.5e0_dp - w)
1272 GO TO 220
1273
1274130 w1 = bup(b0, a0, y0, x0, n, eps)
1275 b0 = b0 + n
1276131 CALL bgrat(b0, a0, y0, x0, w1, eps, ierr1)
1277 IF (ierr1 > 0) cpabort("Error in BGRAT")
1278 w = 0.5e0_dp + (0.5e0_dp - w1)
1279 GO TO 220
1280
1281140 n = int(b0)
1282 b0 = b0 - n
1283 IF (b0 /= 0.0e0_dp) GO TO 141
1284 n = n - 1
1285 b0 = 1.0e0_dp
1286141 w = bup(b0, a0, y0, x0, n, eps)
1287 IF (x0 > 0.7e0_dp) GO TO 150
1288 w = w + bpser(a0, b0, x0, eps)
1289 w1 = 0.5e0_dp + (0.5e0_dp - w)
1290 GO TO 220
1291
1292150 IF (a0 > 15.0e0_dp) GO TO 151
1293 n = 20
1294 w = w + bup(a0, b0, x0, y0, n, eps)
1295 a0 = a0 + n
1296151 CALL bgrat(a0, b0, x0, y0, w, eps, ierr1)
1297 w1 = 0.5e0_dp + (0.5e0_dp - w)
1298 GO TO 220
1299
1300180 w = basym(a0, b0, lambda, 100.0e0_dp*eps)
1301 w1 = 0.5e0_dp + (0.5e0_dp - w)
1302 GO TO 220
1303
1304! TERMINATION OF THE PROCEDURE
1305
1306220 IF (ind == 0) RETURN
1307 t = w
1308 w = w1
1309 w1 = t
1310
1311 END SUBROUTINE bratio
1312
1313! **************************************************************************************************
1314!> \brief ...
1315!> \param a ...
1316!> \param b ...
1317!> \param x ...
1318!> \param eps ...
1319!> \return ...
1320! **************************************************************************************************
1321 FUNCTION fpser(a, b, x, eps) RESULT(fn_val)
1322!-----------------------------------------------------------------------
1323
1324! EVALUATION OF I (A,B)
1325! X
1326
1327! FOR B .LT. MIN(EPS,EPS*A) AND X .LE. 0.5.
1328
1329!-----------------------------------------------------------------------
1330 REAL(dp), INTENT(IN) :: a, b, x, eps
1331 REAL(dp) :: fn_val
1332
1333 REAL(dp) :: an, c, s, t, tol
1334
1335! SET FPSER = X**A
1336
1337 fn_val = 1.0e0_dp
1338 IF (a > 1.e-3_dp*eps) THEN
1339 fn_val = 0.0e0_dp
1340 t = a*log(x)
1341 IF (t < dxparg(1)) RETURN
1342 fn_val = exp(t)
1343 END IF
1344
1345! NOTE THAT 1/B(A,B) = B
1346
1347 fn_val = (b/a)*fn_val
1348 tol = eps/a
1349 an = a + 1.0e0_dp
1350 t = x
1351 s = t/an
1352 an = an + 1.0e0_dp
1353 t = x*t
1354 c = t/an
1355 s = s + c
1356 DO WHILE (abs(c) > tol)
1357 an = an + 1.0e0_dp
1358 t = x*t
1359 c = t/an
1360 s = s + c
1361 END DO
1362
1363 fn_val = fn_val*(1.0e0_dp + a*s)
1364
1365 END FUNCTION fpser
1366
1367! **************************************************************************************************
1368!> \brief ...
1369!> \param a ...
1370!> \param b ...
1371!> \param x ...
1372!> \param eps ...
1373!> \return ...
1374! **************************************************************************************************
1375 FUNCTION apser(a, b, x, eps) RESULT(fn_val)
1376!-----------------------------------------------------------------------
1377! APSER YIELDS THE INCOMPLETE BETA RATIO I(SUB(1-X))(B,A) FOR
1378! A .LE. MIN(EPS,EPS*B), B*X .LE. 1, AND X .LE. 0.5. USED WHEN
1379! A IS VERY SMALL. USE ONLY IF ABOVE INEQUALITIES ARE SATISFIED.
1380!-----------------------------------------------------------------------
1381 REAL(dp), INTENT(IN) :: a, b, x, eps
1382 REAL(dp) :: fn_val
1383
1384 REAL(dp), PARAMETER :: g = .577215664901533e0_dp
1385
1386 REAL(dp) :: aj, bx, c, j, s, t, tol
1387
1388 bx = b*x
1389 t = x - bx
1390 IF (b*eps > 2.e-2_dp) THEN
1391 c = log(bx) + g + t
1392 ELSE
1393 c = log(x) + psi(b) + g + t
1394 END IF
1395
1396 tol = 5.0e0_dp*eps*abs(c)
1397 j = 1.0e0_dp
1398 s = 0.0e0_dp
1399 j = j + 1.0e0_dp
1400 t = t*(x - bx/j)
1401 aj = t/j
1402 s = s + aj
1403 DO WHILE (abs(aj) > tol)
1404 t = t*(x - bx/j)
1405 aj = t/j
1406 s = s + aj
1407 END DO
1408
1409 fn_val = -a*(c + s)
1410
1411 END FUNCTION apser
1412
1413! **************************************************************************************************
1414!> \brief ...
1415!> \param a ...
1416!> \param b ...
1417!> \param x ...
1418!> \param eps ...
1419!> \return ...
1420! **************************************************************************************************
1421 FUNCTION bpser(a, b, x, eps) RESULT(fn_val)
1422!-----------------------------------------------------------------------
1423! POWER SERIES EXPANSION FOR EVALUATING IX(A,B) WHEN B .LE. 1
1424! OR B*X .LE. 0.7. EPS IS THE TOLERANCE USED.
1425!-----------------------------------------------------------------------
1426 REAL(dp), INTENT(IN) :: a, b, x, eps
1427 REAL(dp) :: fn_val
1428
1429 INTEGER :: i, m
1430 REAL(dp) :: a0, apb, b0, c, n, sum, t, tol, u, w, z
1431
1432 fn_val = 0.0e0_dp
1433 IF (x == 0.0e0_dp) RETURN
1434!-----------------------------------------------------------------------
1435! COMPUTE THE FACTOR X**A/(A*BETA(A,B))
1436!-----------------------------------------------------------------------
1437 a0 = min(a, b)
1438 b0 = max(a, b)
1439 IF (a0 >= 1.0e0_dp) THEN
1440 z = a*log(x) - betaln(a, b)
1441 fn_val = exp(z)/a
1442 ELSEIF (b0 >= 8.0e0_dp) THEN
1443 u = gamln1(a0) + algdiv(a0, b0)
1444 z = a*log(x) - u
1445 fn_val = (a0/a)*exp(z)
1446 ELSEIF (b0 > 1.0e0_dp) THEN
1447! PROCEDURE FOR A0 .LT. 1 AND 1 .LT. B0 .LT. 8
1448
1449 u = gamln1(a0)
1450 m = int(b0 - 1.0e0_dp)
1451 IF (m >= 1) THEN
1452 c = 1.0e0_dp
1453 DO i = 1, m
1454 b0 = b0 - 1.0e0_dp
1455 c = c*(b0/(a0 + b0))
1456 END DO
1457 u = log(c) + u
1458 END IF
1459
1460 z = a*log(x) - u
1461 b0 = b0 - 1.0e0_dp
1462 apb = a0 + b0
1463 IF (apb > 1.0e0_dp) THEN
1464 u = a0 + b0 - 1.e0_dp
1465 t = (1.0e0_dp + gam1(u))/apb
1466 ELSE
1467 t = 1.0e0_dp + gam1(apb)
1468 END IF
1469 fn_val = exp(z)*(a0/a)*(1.0e0_dp + gam1(b0))/t
1470 ELSE
1471
1472 ! PROCEDURE FOR A0 .LT. 1 AND B0 .LE. 1
1473
1474 fn_val = x**a
1475 IF (fn_val == 0.0e0_dp) RETURN
1476
1477 apb = a + b
1478 IF (apb > 1.0e0_dp) THEN
1479 u = a + b - 1.e0_dp
1480 z = (1.0e0_dp + gam1(u))/apb
1481 ELSE
1482 z = 1.0e0_dp + gam1(apb)
1483 END IF
1484
1485 c = (1.0e0_dp + gam1(a))*(1.0e0_dp + gam1(b))/z
1486 fn_val = fn_val*c*(b/apb)
1487 END IF
1488
1489! PROCEDURE FOR A0 .LT. 1 AND B0 .GE. 8
1490
1491 IF (fn_val == 0.0e0_dp .OR. a <= 0.1e0_dp*eps) RETURN
1492!-----------------------------------------------------------------------
1493! COMPUTE THE SERIES
1494!-----------------------------------------------------------------------
1495 sum = 0.0e0_dp
1496 n = 0.0e0_dp
1497 c = 1.0e0_dp
1498 tol = eps/a
1499 n = n + 1.0e0_dp
1500 c = c*(0.5e0_dp + (0.5e0_dp - b/n))*x
1501 w = c/(a + n)
1502 sum = sum + w
1503 DO WHILE (abs(w) > tol)
1504 n = n + 1.0e0_dp
1505 c = c*(0.5e0_dp + (0.5e0_dp - b/n))*x
1506 w = c/(a + n)
1507 sum = sum + w
1508 END DO
1509 fn_val = fn_val*(1.0e0_dp + a*sum)
1510
1511 END FUNCTION bpser
1512
1513! **************************************************************************************************
1514!> \brief ...
1515!> \param a ...
1516!> \param b ...
1517!> \param x ...
1518!> \param y ...
1519!> \param n ...
1520!> \param eps ...
1521!> \return ...
1522! **************************************************************************************************
1523 FUNCTION bup(a, b, x, y, n, eps) RESULT(fn_val)
1524!-----------------------------------------------------------------------
1525! EVALUATION OF IX(A,B) - IX(A+N,B) WHERE N IS A POSITIVE INTEGER.
1526! EPS IS THE TOLERANCE USED.
1527!-----------------------------------------------------------------------
1528 REAL(dp), INTENT(IN) :: a, b, x, y
1529 INTEGER, INTENT(IN) :: n
1530 REAL(dp), INTENT(IN) :: eps
1531 REAL(dp) :: fn_val
1532
1533 INTEGER :: i, k, kp1, mu, nm1
1534 REAL(dp) :: ap1, apb, d, l, r, t, w
1535
1536! OBTAIN THE SCALING FACTOR EXP(-MU) AND
1537! EXP(MU)*(X**A*Y**B/BETA(A,B))/A
1538
1539 apb = a + b
1540 ap1 = a + 1.0e0_dp
1541 mu = 0
1542 d = 1.0e0_dp
1543 IF (.NOT. (n == 1 .OR. a < 1.0e0_dp .OR. apb < 1.1e0_dp*ap1)) THEN
1544 mu = int(abs(dxparg(1)))
1545 k = int(dxparg(0))
1546 IF (k < mu) mu = k
1547 t = mu
1548 d = exp(-t)
1549 END IF
1550
1551 fn_val = brcmp1(mu, a, b, x, y)/a
1552 IF (n == 1 .OR. fn_val == 0.0e0_dp) RETURN
1553 nm1 = n - 1
1554 w = d
1555
1556! LET K BE THE INDEX OF THE MAXIMUM TERM
1557
1558 k = 0
1559 IF (b > 1.0e0_dp) THEN
1560 IF (y <= 1.e-4_dp) THEN
1561 k = nm1
1562 DO i = 1, k
1563 l = i - 1
1564 d = ((apb + l)/(ap1 + l))*x*d
1565 w = w + d
1566 END DO
1567 IF (k == nm1) THEN
1568 fn_val = fn_val*w
1569 RETURN
1570 END IF
1571 ELSE
1572 r = (b - 1.0e0_dp)*x/y - a
1573 IF (r >= 1.0e0_dp) THEN
1574 k = nm1
1575 t = nm1
1576 IF (r < t) k = int(r)
1577
1578! ADD THE INCREASING TERMS OF THE SERIES
1579
1580 DO i = 1, k
1581 l = i - 1
1582 d = ((apb + l)/(ap1 + l))*x*d
1583 w = w + d
1584 END DO
1585 IF (k == nm1) THEN
1586 fn_val = fn_val*w
1587 RETURN
1588 END IF
1589 END IF
1590 END IF
1591 END IF
1592
1593! ADD THE REMAINING TERMS OF THE SERIES
1594
1595 kp1 = k + 1
1596 DO i = kp1, nm1
1597 l = i - 1
1598 d = ((apb + l)/(ap1 + l))*x*d
1599 w = w + d
1600 IF (d <= eps*w) EXIT
1601 END DO
1602
1603! TERMINATE THE PROCEDURE
1604
1605 fn_val = fn_val*w
1606
1607 END FUNCTION bup
1608
1609! **************************************************************************************************
1610!> \brief ...
1611!> \param a ...
1612!> \param b ...
1613!> \param x ...
1614!> \param y ...
1615!> \param lambda ...
1616!> \param eps ...
1617!> \return ...
1618! **************************************************************************************************
1619 FUNCTION bfrac(a, b, x, y, lambda, eps) RESULT(fn_val)
1620!-----------------------------------------------------------------------
1621! CONTINUED FRACTION EXPANSION FOR IX(A,B) WHEN A,B .GT. 1.
1622! IT IS ASSUMED THAT LAMBDA = (A + B)*Y - B.
1623!-----------------------------------------------------------------------
1624 REAL(dp), INTENT(IN) :: a, b, x, y, lambda, eps
1625 REAL(dp) :: fn_val
1626
1627 REAL(dp) :: alpha, an, anp1, beta, bn, bnp1, c, c0, &
1628 c1, e, n, p, r, r0, s, t, w, yp1
1629
1630 fn_val = brcomp(a, b, x, y)
1631 IF (fn_val == 0.0e0_dp) RETURN
1632
1633 c = 1.0e0_dp + lambda
1634 c0 = b/a
1635 c1 = 1.0e0_dp + 1.0e0_dp/a
1636 yp1 = y + 1.0e0_dp
1637
1638 n = 0.0e0_dp
1639 p = 1.0e0_dp
1640 s = a + 1.0e0_dp
1641 an = 0.0e0_dp
1642 bn = 1.0e0_dp
1643 anp1 = 1.0e0_dp
1644 bnp1 = c/c1
1645 r = c1/c
1646
1647! CONTINUED FRACTION CALCULATION
1648
1649 DO WHILE (.true.)
1650 n = n + 1.0e0_dp
1651 t = n/a
1652 w = n*(b - n)*x
1653 e = a/s
1654 alpha = (p*(p + c0)*e*e)*(w*x)
1655 IF (alpha <= 0.0e0_dp) THEN
1656! TERMINATION
1657 fn_val = fn_val*r
1658 RETURN
1659 END IF
1660 e = (1.0e0_dp + t)/(c1 + t + t)
1661 beta = n + w/s + e*(c + n*yp1)
1662 p = 1.0e0_dp + t
1663 s = s + 2.0e0_dp
1664
1665! UPDATE AN, BN, ANP1, AND BNP1
1666
1667 t = alpha*an + beta*anp1
1668 an = anp1
1669 anp1 = t
1670 t = alpha*bn + beta*bnp1
1671 bn = bnp1
1672 bnp1 = t
1673 r0 = r
1674 r = anp1/bnp1
1675 IF (abs(r - r0) <= eps*r) THEN
1676! TERMINATION
1677 fn_val = fn_val*r
1678 RETURN
1679 END IF
1680
1681! RESCALE AN, BN, ANP1, AND BNP1
1682
1683 an = an/bnp1
1684 bn = bn/bnp1
1685 anp1 = r
1686 bnp1 = 1.0e0_dp
1687 END DO
1688
1689 END FUNCTION bfrac
1690
1691! **************************************************************************************************
1692!> \brief ...
1693!> \param a ...
1694!> \param b ...
1695!> \param x ...
1696!> \param y ...
1697!> \return ...
1698! **************************************************************************************************
1699 FUNCTION brcomp(a, b, x, y) RESULT(fn_val)
1700!-----------------------------------------------------------------------
1701! EVALUATION OF X**A*Y**B/BETA(A,B)
1702!-----------------------------------------------------------------------
1703 REAL(dp), INTENT(IN) :: a, b, x, y
1704 REAL(dp) :: fn_val
1705
1706 REAL(dp), PARAMETER :: const = 0.398942280401433e0_dp
1707
1708 INTEGER :: i, n
1709 REAL(dp) :: a0, apb, b0, c, e, h, lambda, lnx, lny, &
1710 t, u, v, x0, y0, z
1711
1712!-----------------
1713! CONST = 1/SQRT(2*PI)
1714!-----------------
1715
1716 fn_val = 0.0e0_dp
1717 IF (x == 0.0e0_dp .OR. y == 0.0e0_dp) RETURN
1718 a0 = min(a, b)
1719 IF (a0 >= 8.0e0_dp) THEN
1720!-----------------------------------------------------------------------
1721! PROCEDURE FOR A .GE. 8 AND B .GE. 8
1722!-----------------------------------------------------------------------
1723 IF (a > b) THEN
1724 h = b/a
1725 x0 = 1.0e0_dp/(1.0e0_dp + h)
1726 y0 = h/(1.0e0_dp + h)
1727 lambda = (a + b)*y - b
1728 ELSE
1729 h = a/b
1730 x0 = h/(1.0e0_dp + h)
1731 y0 = 1.0e0_dp/(1.0e0_dp + h)
1732 lambda = a - (a + b)*x
1733 END IF
1734
1735 e = -lambda/a
1736 IF (abs(e) > 0.6e0_dp) THEN
1737 u = e - log(x/x0)
1738 ELSE
1739 u = rlog1(e)
1740 END IF
1741
1742 e = lambda/b
1743 IF (abs(e) > 0.6e0_dp) THEN
1744 v = e - log(y/y0)
1745 ELSE
1746 v = rlog1(e)
1747 END IF
1748
1749 z = exp(-(a*u + b*v))
1750 fn_val = const*sqrt(b*x0)*z*exp(-bcorr(a, b))
1751 RETURN
1752 END IF
1753
1754 IF (x > 0.375e0_dp) THEN
1755 IF (y > 0.375e0_dp) THEN
1756 lnx = log(x)
1757 lny = log(y)
1758 ELSE
1759 lnx = alnrel(-y)
1760 lny = log(y)
1761 END IF
1762 ELSE
1763 lnx = log(x)
1764 lny = alnrel(-x)
1765 END IF
1766
1767 z = a*lnx + b*lny
1768 IF (a0 >= 1.0e0_dp) THEN
1769 z = z - betaln(a, b)
1770 fn_val = exp(z)
1771 RETURN
1772 END IF
1773!-----------------------------------------------------------------------
1774! PROCEDURE FOR A .LT. 1 OR B .LT. 1
1775!-----------------------------------------------------------------------
1776 b0 = max(a, b)
1777 IF (b0 >= 8.0e0_dp) THEN
1778! ALGORITHM FOR B0 .GE. 8
1779 u = gamln1(a0) + algdiv(a0, b0)
1780 fn_val = a0*exp(z - u)
1781 END IF
1782 IF (b0 <= 1.0e0_dp) THEN
1783
1784! ALGORITHM FOR B0 .LE. 1
1785
1786 fn_val = exp(z)
1787 IF (fn_val == 0.0e0_dp) RETURN
1788
1789 apb = a + b
1790 IF (apb > 1.0e0_dp) THEN
1791 u = a + b - 1.e0_dp
1792 z = (1.0e0_dp + gam1(u))/apb
1793 ELSE
1794 z = 1.0e0_dp + gam1(apb)
1795 END IF
1796
1797 c = (1.0e0_dp + gam1(a))*(1.0e0_dp + gam1(b))/z
1798 fn_val = fn_val*(a0*c)/(1.0e0_dp + a0/b0)
1799 RETURN
1800 END IF
1801
1802! ALGORITHM FOR 1 .LT. B0 .LT. 8
1803
1804 u = gamln1(a0)
1805 n = int(b0 - 1.0e0_dp)
1806 IF (n >= 1) THEN
1807 c = 1.0e0_dp
1808 DO i = 1, n
1809 b0 = b0 - 1.0e0_dp
1810 c = c*(b0/(a0 + b0))
1811 END DO
1812 u = log(c) + u
1813 END IF
1814
1815 z = z - u
1816 b0 = b0 - 1.0e0_dp
1817 apb = a0 + b0
1818 IF (apb > 1.0e0_dp) THEN
1819 u = a0 + b0 - 1.e0_dp
1820 t = (1.0e0_dp + gam1(u))/apb
1821 ELSE
1822 t = 1.0e0_dp + gam1(apb)
1823 END IF
1824 fn_val = a0*exp(z)*(1.0e0_dp + gam1(b0))/t
1825
1826 END FUNCTION brcomp
1827
1828! **************************************************************************************************
1829!> \brief ...
1830!> \param mu ...
1831!> \param a ...
1832!> \param b ...
1833!> \param x ...
1834!> \param y ...
1835!> \return ...
1836! **************************************************************************************************
1837 FUNCTION brcmp1(mu, a, b, x, y) RESULT(fn_val)
1838!-----------------------------------------------------------------------
1839! EVALUATION OF EXP(MU) * (X**A*Y**B/BETA(A,B))
1840!-----------------------------------------------------------------------
1841 INTEGER, INTENT(IN) :: mu
1842 REAL(dp), INTENT(IN) :: a, b, x, y
1843 REAL(dp) :: fn_val
1844
1845 REAL(dp), PARAMETER :: const = 0.398942280401433e0_dp
1846
1847 INTEGER :: i, n
1848 REAL(dp) :: a0, apb, b0, c, e, h, lambda, lnx, lny, &
1849 t, u, v, x0, y0, z
1850
1851!-----------------
1852! CONST = 1/SQRT(2*PI)
1853!-----------------
1854
1855 a0 = min(a, b)
1856 IF (a0 >= 8.0e0_dp) THEN
1857!-----------------------------------------------------------------------
1858! PROCEDURE FOR A .GE. 8 AND B .GE. 8
1859!-----------------------------------------------------------------------
1860 IF (a > b) THEN
1861 h = b/a
1862 x0 = 1.0e0_dp/(1.0e0_dp + h)
1863 y0 = h/(1.0e0_dp + h)
1864 lambda = (a + b)*y - b
1865 END IF
1866 h = a/b
1867 x0 = h/(1.0e0_dp + h)
1868 y0 = 1.0e0_dp/(1.0e0_dp + h)
1869 lambda = a - (a + b)*x
1870
1871 e = -lambda/a
1872 IF (abs(e) > 0.6e0_dp) THEN
1873 u = e - log(x/x0)
1874 ELSE
1875 u = rlog1(e)
1876 END IF
1877
1878 e = lambda/b
1879 IF (abs(e) <= 0.6e0_dp) THEN
1880 v = rlog1(e)
1881 ELSE
1882 v = e - log(y/y0)
1883 END IF
1884
1885 z = esum(mu, -(a*u + b*v))
1886 fn_val = const*sqrt(b*x0)*z*exp(-bcorr(a, b))
1887 END IF
1888
1889 IF (x > 0.375e0_dp) THEN
1890 IF (y > 0.375e0_dp) THEN
1891 lnx = log(x)
1892 lny = log(y)
1893 ELSE
1894 lnx = alnrel(-y)
1895 lny = log(y)
1896 END IF
1897 ELSE
1898 lnx = log(x)
1899 lny = alnrel(-x)
1900 END IF
1901 z = a*lnx + b*lny
1902 IF (a0 >= 1.0e0_dp) THEN
1903 z = z - betaln(a, b)
1904 fn_val = esum(mu, z)
1905 RETURN
1906 END IF
1907!-----------------------------------------------------------------------
1908! PROCEDURE FOR A .LT. 1 OR B .LT. 1
1909!-----------------------------------------------------------------------
1910 b0 = max(a, b)
1911 IF (b0 >= 8.0e0_dp) THEN
1912! ALGORITHM FOR B0 .GE. 8
1913
1914 u = gamln1(a0) + algdiv(a0, b0)
1915 fn_val = a0*esum(mu, z - u)
1916 RETURN
1917 END IF
1918 IF (b0 <= 1.0e0_dp) THEN
1919
1920! ALGORITHM FOR B0 .LE. 1
1921
1922 fn_val = esum(mu, z)
1923 IF (fn_val == 0.0e0_dp) RETURN
1924
1925 apb = a + b
1926 IF (apb > 1.0e0_dp) THEN
1927 u = a + b - 1.e0_dp
1928 z = (1.0e0_dp + gam1(u))/apb
1929 ELSE
1930 z = 1.0e0_dp + gam1(apb)
1931 END IF
1932
1933 c = (1.0e0_dp + gam1(a))*(1.0e0_dp + gam1(b))/z
1934 fn_val = fn_val*(a0*c)/(1.0e0_dp + a0/b0)
1935 RETURN
1936 END IF
1937
1938! ALGORITHM FOR 1 .LT. B0 .LT. 8
1939
1940 u = gamln1(a0)
1941 n = int(b0 - 1.0e0_dp)
1942 IF (n >= 1) THEN
1943 c = 1.0e0_dp
1944 DO i = 1, n
1945 b0 = b0 - 1.0e0_dp
1946 c = c*(b0/(a0 + b0))
1947 END DO
1948 u = log(c) + u
1949 END IF
1950
1951 z = z - u
1952 b0 = b0 - 1.0e0_dp
1953 apb = a0 + b0
1954 IF (apb > 1.0e0_dp) THEN
1955 u = a0 + b0 - 1.e0_dp
1956 t = (1.0e0_dp + gam1(u))/apb
1957 ELSE
1958 t = 1.0e0_dp + gam1(apb)
1959 END IF
1960 fn_val = a0*esum(mu, z)*(1.0e0_dp + gam1(b0))/t
1961
1962 END FUNCTION brcmp1
1963
1964! **************************************************************************************************
1965!> \brief ...
1966!> \param a ...
1967!> \param b ...
1968!> \param lambda ...
1969!> \param eps ...
1970!> \return ...
1971! **************************************************************************************************
1972 FUNCTION basym(a, b, lambda, eps) RESULT(fn_val)
1973!-----------------------------------------------------------------------
1974! ASYMPTOTIC EXPANSION FOR IX(A,B) FOR LARGE A AND B.
1975! LAMBDA = (A + B)*Y - B AND EPS IS THE TOLERANCE USED.
1976! IT IS ASSUMED THAT LAMBDA IS NONNEGATIVE AND THAT
1977! A AND B ARE GREATER THAN OR EQUAL TO 15.
1978!-----------------------------------------------------------------------
1979 REAL(dp), INTENT(IN) :: a, b, lambda, eps
1980 REAL(dp) :: fn_val
1981
1982 INTEGER, PARAMETER :: num = 20
1983 REAL(dp), PARAMETER :: e0 = 1.12837916709551e0_dp, &
1984 e1 = .353553390593274e0_dp
1985
1986 INTEGER :: i, im1, imj, j, m, mm1, mmj, n, np1
1987 REAL(dp) :: a0(21), b0(21), bsum, c(21), d(21), &
1988 dsum, f, h, h2, hn, j0, j1, r, r0, r1, &
1989 s, sum, t, t0, t1, u, w, w0, z, z0, &
1990 z2, zn, znm1
1991
1992!------------------------
1993! ****** NUM IS THE MAXIMUM VALUE THAT N CAN TAKE IN THE DO LOOP
1994! ENDING AT STATEMENT 50. IT IS REQUIRED THAT NUM BE EVEN.
1995! THE ARRAYS A0, B0, C, D HAVE DIMENSION NUM + 1.
1996!------------------------
1997! E0 = 2/SQRT(PI)
1998! E1 = 2**(-3/2)
1999!------------------------
2000
2001 fn_val = 0.0e0_dp
2002 IF (a >= b) THEN
2003 h = b/a
2004 r0 = 1.0e0_dp/(1.0e0_dp + h)
2005 r1 = (b - a)/a
2006 w0 = 1.0e0_dp/sqrt(b*(1.0e0_dp + h))
2007 ELSE
2008 h = a/b
2009 r0 = 1.0e0_dp/(1.0e0_dp + h)
2010 r1 = (b - a)/b
2011 w0 = 1.0e0_dp/sqrt(a*(1.0e0_dp + h))
2012 END IF
2013
2014 f = a*rlog1(-lambda/a) + b*rlog1(lambda/b)
2015 t = exp(-f)
2016 IF (t == 0.0e0_dp) RETURN
2017 z0 = sqrt(f)
2018 z = 0.5e0_dp*(z0/e1)
2019 z2 = f + f
2020
2021 a0(1) = (2.0e0_dp/3.0e0_dp)*r1
2022 c(1) = -0.5e0_dp*a0(1)
2023 d(1) = -c(1)
2024 j0 = (0.5e0_dp/e0)*erfc_scaled(z0)
2025 j1 = e1
2026 sum = j0 + d(1)*w0*j1
2027
2028 s = 1.0e0_dp
2029 h2 = h*h
2030 hn = 1.0e0_dp
2031 w = w0
2032 znm1 = z
2033 zn = z2
2034 DO n = 2, num, 2
2035 hn = h2*hn
2036 a0(n) = 2.0e0_dp*r0*(1.0e0_dp + h*hn)/(n + 2.0e0_dp)
2037 np1 = n + 1
2038 s = s + hn
2039 a0(np1) = 2.0e0_dp*r1*s/(n + 3.0e0_dp)
2040
2041 DO i = n, np1
2042 r = -0.5e0_dp*(i + 1.0e0_dp)
2043 b0(1) = r*a0(1)
2044 DO m = 2, i
2045 bsum = 0.0e0_dp
2046 mm1 = m - 1
2047 DO j = 1, mm1
2048 mmj = m - j
2049 bsum = bsum + (j*r - mmj)*a0(j)*b0(mmj)
2050 END DO
2051 b0(m) = r*a0(m) + bsum/m
2052 END DO
2053 c(i) = b0(i)/(i + 1.0e0_dp)
2054
2055 dsum = 0.0e0_dp
2056 im1 = i - 1
2057 DO j = 1, im1
2058 imj = i - j
2059 dsum = dsum + d(imj)*c(j)
2060 END DO
2061 d(i) = -(dsum + c(i))
2062 END DO
2063
2064 j0 = e1*znm1 + (n - 1.0e0_dp)*j0
2065 j1 = e1*zn + n*j1
2066 znm1 = z2*znm1
2067 zn = z2*zn
2068 w = w0*w
2069 t0 = d(n)*w*j0
2070 w = w0*w
2071 t1 = d(np1)*w*j1
2072 sum = sum + (t0 + t1)
2073 IF ((abs(t0) + abs(t1)) <= eps*sum) EXIT
2074 END DO
2075
2076 u = exp(-bcorr(a, b))
2077 fn_val = e0*t*u*sum
2078
2079 END FUNCTION basym
2080
2081END MODULE beta_gamma_psi
beta_gamma_psi
Definition beta_gamma_psi.F:9
beta_gamma_psi::psi
real(dp) function, public psi(xx)
...
Definition beta_gamma_psi.F:718
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34