(git:b279b6b)
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 ! **************************************************************************************************
9 MODULE 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 
22 CONTAINS
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 
1218 20 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 
1227 30 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
1248 50 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 
1254 80 w = fpser(a0, b0, x0, eps)
1255  w1 = 0.5e0_dp + (0.5e0_dp - w)
1256  GO TO 220
1257 
1258 90 w1 = apser(a0, b0, x0, eps)
1259  w = 0.5e0_dp + (0.5e0_dp - w1)
1260  GO TO 220
1261 
1262 100 w = bpser(a0, b0, x0, eps)
1263  w1 = 0.5e0_dp + (0.5e0_dp - w)
1264  GO TO 220
1265 
1266 110 w1 = bpser(b0, a0, y0, eps)
1267  w = 0.5e0_dp + (0.5e0_dp - w1)
1268  GO TO 220
1269 
1270 120 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 
1274 130 w1 = bup(b0, a0, y0, x0, n, eps)
1275  b0 = b0 + n
1276 131 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 
1281 140 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
1286 141 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 
1292 150 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
1296 151 CALL bgrat(a0, b0, x0, y0, w, eps, ierr1)
1297  w1 = 0.5e0_dp + (0.5e0_dp - w)
1298  GO TO 220
1299 
1300 180 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 
1306 220 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 
2081 END 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