mchgu.f90

数值计算和数值分析在Fortran下的特殊函数库,是数值计算的必备

MODULE chgu_func
 
! From the book "Computation of Special Functions"
!      by Shanjie Zhang and Jianming Jin
!   Copyright 1996 by John Wiley & Sons, Inc.
! The authors state:
!   "However, we give permission to the reader who purchases this book
!    to incorporate any of these programs into his or her programs
!    provided that the copyright is acknowledged."

! Latest revision - 27 December 2001
! Corrections by Alan Miller (amiller @ bigpond.net.au)
! Variables h0 & r0 were used without values assigned to them

IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(12, 60)
 
CONTAINS


SUBROUTINE chgu(a, b, x, hu, md)

!       =======================================================
!       Purpose: Compute the confluent hypergeometric function
!                U(a,b,x)
!       Input  : a  --- Parameter
!                b  --- Parameter
!                x  --- Argument  ( x > 0 )
!       Output:  HU --- U(a,b,x)
!                MD --- Method code
!       Routines called:
!            (1) CHGUS for small x ( MD=1 )
!            (2) CHGUL for large x ( MD=2 )
!            (3) CHGUBI for integer b ( MD=3 )
!            (4) CHGUIT for numerical integration ( MD=4 )
!       =======================================================

REAL (dp), INTENT(IN OUT)  :: a
REAL (dp), INTENT(IN OUT)  :: b
REAL (dp), INTENT(IN)      :: x
REAL (dp), INTENT(IN OUT)  :: hu
INTEGER, INTENT(OUT)       :: md

LOGICAL    :: il1, il2, il3, bl1, bl2, bl3, bn
REAL (dp)  :: a00, aa, b00, hu1
INTEGER    :: id, id1

aa = a - b + 1.0_dp
il1 = a == INT(a) .AND. a <= 0.0
il2 = aa == INT(aa) .AND. aa <= 0.0
il3 = ABS(a*(a-b+1.0)) / x <= 2.0
bl1 = x <= 5.0 .OR. (x <= 10.0 .AND. a <= 2.0)
bl2 = (x > 5.0 .AND. x <= 12.5) .AND. (a >= 1.0 .AND. b >= a+4.0)
bl3 = x > 12.5 .AND. a >= 5.0 .AND. b >= a + 5.0
bn = b == INT(b) .AND. b /= 0.0
id1 = -100
IF (b /= INT(b)) THEN
  CALL chgus(a, b, x, hu, id1)
  md = 1
  IF (id1 >= 6) RETURN
  hu1 = hu
END IF
IF (il1 .OR. il2 .OR. il3) THEN
  CALL chgul(a, b, x, hu, id)
  md = 2
  IF (id >= 6) RETURN
  IF (id1 > id) THEN
    md = 1
    id = id1
    hu = hu1
  END IF
END IF
IF (a >= 0.0) THEN
  IF (bn .AND. (bl1 .OR. bl2 .OR. bl3)) THEN
    CALL chgubi(a, b, x, hu, id)
    md = 3
  ELSE
    CALL chguit(a, b, x, hu, id)
    md = 4
  END IF
ELSE
  IF (b <= a) THEN
    a00 = a
    b00 = b
    a = a - b + 1.0_dp
    b = 2.0_dp - b
    CALL chguit(a, b, x, hu, id)
    hu = x ** (1.0_dp-b00) * hu
    a = a00
    b = b00
    md = 4
  ELSE IF (bn .AND. (.NOT.il1)) THEN
    CALL chgubi(a, b, x, hu, id)
    md = 3
  END IF
END IF
IF (id < 6) WRITE (*,*) 'No accurate result obtained'
RETURN
END SUBROUTINE chgu



SUBROUTINE chgus(a, b, x, hu, id)

!       ======================================================
!       Purpose: Compute confluent hypergeometric function
!                U(a,b,x) for small argument x
!       Input  : a  --- Parameter
!                b  --- Parameter ( b <> 0,-1,-2,...)
!                x  --- Argument
!       Output:  HU --- U(a,b,x)
!                ID --- Estimated number of significant digits
!       Routine called: GAMMA for computing gamma function
!       ======================================================

REAL (dp), INTENT(IN)      :: a
REAL (dp), INTENT(IN)      :: b
REAL (dp), INTENT(IN)      :: x
REAL (dp), INTENT(OUT)     :: hu
INTEGER, INTENT(OUT)       :: id

REAL (dp), PARAMETER  :: pi = 3.141592653589793_dp
REAL (dp)  :: d1, d2, ga, gab, gb, gb2, h0, hmax, hmin, hu0, hua,  &
              r1, r2, xg1, xg2
INTEGER    :: j

id = -100
CALL gamma(a, ga)
CALL gamma(b, gb)
xg1 = 1.0_dp + a - b
CALL gamma(xg1, gab)
xg2 = 2.0_dp - b
CALL gamma(xg2, gb2)
hu0 = pi / SIN(pi*b)
r1 = hu0 / (gab*gb)
r2 = hu0 * x ** (1.0_dp-b) / (ga*gb2)
hu = r1 - r2
h0 = hu
hmax = 0.0_dp
hmin = 1.0D+300
DO  j = 1, 150
  r1 = r1 * (a+j-1.0_dp) / (j*(b+j-1.0_dp)) * x
  r2 = r2 * (a-b+j) / (j*(1.0_dp-b+j)) * x
  hu = hu + r1 - r2
  hua = ABS(hu)
  IF (hua > hmax) hmax = hua
  IF (hua < hmin) hmin = hua
  IF (ABS(hu-h0) < ABS(hu)*1.0D-15) EXIT
  h0 = hu
END DO

d1 = LOG10(hmax)
IF (hmin /= 0.0) d2 = LOG10(hmin)
id = 15 - ABS(d1-d2)
RETURN
END SUBROUTINE chgus



SUBROUTINE chgul(a, b, x, hu, id)

!       =======================================================
!       Purpose: Compute the confluent hypergeometric function
!                U(a,b,x) for large argument x
!       Input  : a  --- Parameter
!                b  --- Parameter
!                x  --- Argument
!       Output:  HU --- U(a,b,x)
!                ID --- Estimated number of significant digits
!       =======================================================

REAL (dp), INTENT(IN)      :: a
REAL (dp), INTENT(IN)      :: b
REAL (dp), INTENT(IN)      :: x
REAL (dp), INTENT(OUT)     :: hu
INTEGER, INTENT(OUT)       :: id

LOGICAL    :: il1, il2
REAL (dp)  :: aa, r, r0, ra
INTEGER    :: k, nm

id = -100
aa = a - b + 1.0_dp
il1 = a == INT(a) .AND. a <= 0.0
il2 = aa == INT(aa) .AND. aa <= 0.0
IF (il1) nm = ABS(a)
IF (il2) nm = ABS(aa)
IF (il1 .OR. il2) THEN
  hu = 1.0_dp
  r = 1.0_dp
  DO  k = 1, nm
    r = -r * (a+k-1.0_dp) * (a-b+k) / (k*x)
    hu = hu + r
  END DO
  hu = x ** (-a) * hu
  id = 10
ELSE
  hu = 1.0_dp
  r = 1.0_dp
  r0 = r
  DO  k = 1, 25
    r = -r * (a+k-1.0_dp) * (a-b+k) / (k*x)
    ra = ABS(r)
    IF (k > 5 .AND. ra >= r0 .OR. ra < 1.0D-15) EXIT
    r0 = ra
    hu = hu + r
  END DO
  id = ABS(LOG10(ra))
  hu = x ** (-a) * hu
END IF
RETURN
END SUBROUTINE chgul



SUBROUTINE chgubi(a, b, x, hu, id)

!       ======================================================
!       Purpose: Compute confluent hypergeometric function
!                U(a,b,x) with integer b ( b =