mcsphik.f90

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

MODULE csphik_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 csphik(n, z, nm, csi, cdi, csk, cdk)

!       =======================================================
!       Purpose: Compute modified spherical Bessel functions
!                and their derivatives for a complex argument
!       Input :  z --- Complex argument
!                n --- Order of in(z) & kn(z) ( n = 0,1,2,... )
!       Output:  CSI(n) --- in(z)
!                CDI(n) --- in'(z)
!                CSK(n) --- kn(z)
!                CDK(n) --- kn'(z)
!                NM --- Highest order computed
!       Routines called:
!                MSTA1 and MSTA2 for computing the starting
!                point for backward recurrence
!       =======================================================

INTEGER, INTENT(IN)        :: n
COMPLEX (dp), INTENT(IN)   :: z
INTEGER, INTENT(OUT)       :: nm
COMPLEX (dp), INTENT(OUT)  :: csi(0:n)
COMPLEX (dp), INTENT(OUT)  :: cdi(0:n)
COMPLEX (dp), INTENT(OUT)  :: csk(0:n)
COMPLEX (dp), INTENT(OUT)  :: cdk(0:n)

REAL (dp)     :: a0
COMPLEX (dp)  :: ccosh, cf, cf0, cf1, ci, cs, csi0, csi1, csinh
INTEGER       :: k, m
REAL (dp), PARAMETER  :: pi = 3.141592653589793D0

a0 = ABS(z)
nm = n
IF (a0 < 1.0D-60) THEN
  DO  k = 0, n
    csi(k) = 0.0D0
    cdi(k) = 0.0D0
    csk(k) = 1.0D+300
    cdk(k) = -1.0D+300
  END DO
  csi(0) = 1.0D0
  cdi(1) = 0.3333333333333333D0
  RETURN
END IF
ci = CMPLX(0.0D0, 1.0D0, KIND=dp)
csinh = SIN(ci*z) / ci
ccosh = COS(ci*z)
csi0 = csinh / z
csi1 = (-csinh/z + ccosh) / z
csi(0) = csi0
csi(1) = csi1
IF (n >= 2) THEN
  m = msta1(a0, 200)
  IF (m < n) THEN
    nm = m
  ELSE
    m = msta2(a0, n, 15)
  END IF
  cf0 = 0.0D0
  cf1 = 1.0D0 - 100
  DO  k = m, 0, -1
    cf = (2*k+3) * cf1 / z + cf0
    IF (k <= nm) csi(k) = cf
    cf0 = cf1
    cf1 = cf
  END DO
  IF (ABS(csi0) > ABS(csi1)) cs = csi0 / cf
  IF (ABS(csi0) <= ABS(csi1)) cs = csi1 / cf0
  DO  k = 0, nm
    csi(k) = cs * csi(k)
  END DO
END IF
cdi(0) = csi(1)
DO  k = 1, nm
  cdi(k) = csi(k-1) - (k+1) * csi(k) / z
END DO
csk(0) = 0.5D0 * pi / z * EXP(-z)
csk(1) = csk(0) * (1.0D0+1.0D0/z)
DO  k = 2, nm
  IF (ABS(csi(k-1)) > ABS(csi(k-2))) THEN
    csk(k) = (0.5D0*pi/(z*z) - csi(k)*csk(k-1)) / csi(k-1)
  ELSE
    csk(k) = (csi(k)*csk(k-2) + (k-0.5D0)*pi/z**3) / csi(k-2)
  END IF
END DO
cdk(0) = -csk(1)
DO  k = 1, nm
  cdk(k) = -csk(k-1) - (k+1) * csk(k) / z
END DO
RETURN
END SUBROUTINE csphik



FUNCTION msta1(x, mp) RESULT(fn_val)

!       ===================================================
!       Purpose: Determine the starting point for backward
!                recurrence such that the magnitude of
!                Jn(x) at that point is about 10^(-MP)
!       Input :  x     --- Argument of Jn(x)
!                MP    --- Value of magnitude
!       Output:  MSTA1 --- Starting point
!       ===================================================

REAL (dp), INTENT(IN)      :: x
INTEGER, INTENT(IN)        :: mp
INTEGER                    :: fn_val

REAL (dp)  :: a0, f, f0, f1
INTEGER    :: it, n0, n1, nn

a0 = ABS(x)
n0 = INT(1.1*a0) + 1
f0 = envj(n0,a0) - mp
n1 = n0 + 5
f1 = envj(n1,a0) - mp
DO  it = 1, 20
  nn = n1 - (n1-n0) / (1.0_dp - f0/f1)
  f = envj(nn,a0) - mp
  IF (ABS(nn-n1) < 1) EXIT
  n0 = n1
  f0 = f1
  n1 = nn
  f1 = f
END DO

fn_val = nn
RETURN
END FUNCTION msta1



FUNCTION msta2(x, n, mp) RESULT(fn_val)

!       ===================================================
!       Purpose: Determine the starting point for backward
!                recurrence such that all Jn(x) has MP
!                significant digits
!       Input :  x  --- Argument of Jn(x)
!                n  --- Order of Jn(x)
!                MP --- Significant digit
!       Output:  MSTA2 --- Starting point
!       ===================================================

REAL (dp), INTENT(IN)      :: x
INTEGER, INTENT(IN)        :: n
INTEGER, INTENT(IN)        :: mp
INTEGER                    :: fn_val

REAL (dp)  :: a0, ejn, f, f0, f1, hmp, obj
INTEGER    :: it, n0, n1, nn

a0 = ABS(x)
hmp = 0.5_dp * mp
ejn = envj(n, a0)
IF (ejn <= hmp) THEN
  obj = mp
  n0 = INT(1.1*a0)
ELSE
  obj = hmp + ejn
  n0 = n
END IF
f0 = envj(n0,a0) - obj
n1 = n0 + 5
f1 = envj(n1,a0) - obj
DO  it = 1, 20
  nn = n1 - (n1-n0) / (1.0_dp - f0/f1)
  f = envj(nn, a0) - obj
  IF (ABS(nn-n1) < 1) EXIT
  n0 = n1
  f0 = f1
  n1 = nn
  f1 = f
END DO

fn_val = nn + 10
RETURN
END FUNCTION msta2



FUNCTION envj(n, x) RESULT(fn_val)

INTEGER, INTENT(IN)        :: n
REAL (dp), INTENT(IN)      :: x
REAL (dp)                  :: fn_val

fn_val = 0.5_dp * LOG10(6.28_dp*n) - n * LOG10(1.36_dp*x/n)
RETURN
END FUNCTION envj

END MODULE csphik_func
 
 
 
PROGRAM mcsphik
USE csphik_func
IMPLICIT NONE

! Code converted using TO_F90 by Alan Miller
! Date: 2001-12-25  Time: 11:55:37

!       =============================================================
!       Purpose: This program computes the modified spherical Bessel
!                functions and their derivatives for a complex
!                argument using subroutine CSPHIK
!       Input :  z --- Complex argument
!                n --- Order of in(z) & kn(z) ( 0