mcchg.f90

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

MODULE cchg_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."
 
IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(12, 60)
 
CONTAINS
 

SUBROUTINE cchg(a, b, z, chg)

!    ===================================================
!    Purpose: Compute confluent hypergeometric function
!             M(a,b,z) with real parameters a, b and a
!             complex argument z
!    Input :  a --- Parameter
!             b --- Parameter
!             z --- Complex argument
!    Output:  CHG --- M(a,b,z)
!    Routine called: GAMMA for computing gamma function
!    ===================================================

REAL (dp), INTENT(IN OUT)     :: a
REAL (dp), INTENT(IN)         :: b
COMPLEX (dp), INTENT(IN OUT)  :: z
COMPLEX (dp), INTENT(OUT)     :: chg

COMPLEX (dp)  :: cfac, chg1, chg2, chw, ci, cr, cr1, cr2, crg, cs1, cs2,  &
                 cy0, cy1, z0
REAL (dp)     :: a0, a1, ba, g1, g2, g3, phi, x, x0, y
INTEGER       :: i, j, k, la, m, n, nl, ns
REAL (dp), PARAMETER  :: pi = 3.141592653589793_dp

ci = (0.0_dp, 1.0_dp)
a0 = a
a1 = a
z0 = z
IF (b == 0.0 .OR. b == -INT(ABS(b))) THEN
  chg = (1.0D+300, 0.0_dp)
ELSE IF (a == 0.0_dp .OR. z == 0.0_dp) THEN
  chg = (1.0_dp, 0.0_dp)
ELSE IF (a == -1.0_dp) THEN
  chg = 1.0_dp - z / b
ELSE IF (a == b) THEN
  chg = EXP(z)
ELSE IF (a-b == 1.0_dp) THEN
  chg = (1.0_dp + z/b) * EXP(z)
ELSE IF (a == 1.0_dp .AND. b == 2.0_dp) THEN
  chg = (EXP(z) - 1.0_dp) / z
ELSE IF (a == INT(a) .AND. a < 0.0_dp) THEN
  m = INT(-a)
  cr = (1.0_dp,0.0_dp)
  chg = (1.0_dp,0.0_dp)
  DO  k = 1, m
    cr = cr * (a+k-1) / k / (b+k-1) * z
    chg = chg + cr
  END DO
ELSE
  x0 = REAL(z)
  IF (x0 < 0.0_dp) THEN
    a = b - a
    a0 = a
    z = -z
  END IF
  IF (a < 2.0_dp) nl = 0
  IF (a >= 2.0_dp) THEN
    nl = 1
    la = INT(a)
    a = a - la - 1.0_dp
  END IF
  DO  n = 0, nl
    IF (a0 >= 2.0_dp) a = a + 1.0_dp
    IF (ABS(z) < 20.0_dp + ABS(b) .OR. a < 0.0_dp) THEN
      chg = (1.0_dp,0.0_dp)
      crg = (1.0_dp,0.0_dp)
      DO  j = 1, 500
        crg = crg * (a+j-1) / (j*(b+j-1)) * z
        chg = chg + crg
        IF (ABS((chg-chw)/chg) < 1.d-15) GO TO 40
        chw = chg
      END DO
    ELSE
      CALL gamma(a, g1)
      CALL gamma(b, g2)
      ba = b - a
      CALL gamma(ba, g3)
      cs1 = (1.0_dp,0.0_dp)
      cs2 = (1.0_dp,0.0_dp)
      cr1 = (1.0_dp,0.0_dp)
      cr2 = (1.0_dp,0.0_dp)
      DO  i = 1, 8
        cr1 = -cr1 * (a+i-1) * (a-b+i) / (z*i)
        cr2 = cr2 * (b-a+i-1) * (i-a) / (z*i)
        cs1 = cs1 + cr1
        cs2 = cs2 + cr2
      END DO
      x = REAL(z)
      y = AIMAG(z)
      IF (x == 0.0 .AND. y >= 0.0) THEN
        phi = 0.5_dp * pi
      ELSE IF (x == 0.0 .AND. y <= 0.0) THEN
        phi = -0.5_dp * pi
      ELSE
        phi = ATAN(y/x)
      END IF
      IF (phi > -0.5*pi .AND. phi < 1.5*pi) ns = 1
      IF (phi > -1.5*pi .AND. phi <= -0.5*pi) ns = -1
      cfac = EXP(ns*ci*pi*a)
      IF (y == 0.0_dp) cfac = COS(pi*a)
      chg1 = g2 / g3 * z ** (-a) * cfac * cs1
      chg2 = g2 / g1 * EXP(z) * z ** (a-b) * cs2
      chg = chg1 + chg2
    END IF
    40 IF (n == 0) cy0 = chg
    IF (n == 1) cy1 = chg
  END DO
  IF (a0 >= 2.0_dp) THEN
    DO  i = 1, la - 1
      chg = ((2.0_dp*a-b+z)*cy1 + (b-a)*cy0) / a
      cy0 = cy1
      cy1 = chg
      a = a + 1.0_dp
    END DO
  END IF
  IF (x0 < 0.0_dp) chg = chg * EXP(-z)
END IF
a = a1
z = z0
RETURN
END SUBROUTINE cchg


SUBROUTINE gamma(x, ga)

!    ==================================================
!    Purpose: Compute gamma function