mrswfp.f90

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

MODULE rswfp_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 - 16 December 2002
! Corrections by Alan Miller (amiller @ bigpond.net.au)
! Variables sw & fl were used without values assigned to them.
! There were several calls to routines with the same argument variable
! repeated.   This is only permitted when the called routine does not alter
! either argument - i.e. when they are both INTENT(IN) arguments.   In each
! case, the duplicated argument was nm2 and the second one has been replaced
! with nm3.
! Also, the frequent use of (-1) raised to an integer power has been replaced
! with a much faster calculation of whether that integer is odd or even.

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


SUBROUTINE rswfp(m, n, c, x, cv, kf, r1f, r1d, r2f, r2d)

!    ==============================================================
!    Purpose: Compute prolate spheriodal radial functions of the
!             first and second kinds, and their derivatives
!    Input :  m  --- Mode parameter, m = 0,1,2,...
!             n  --- Mode parameter, n = m,m+1,m+2,...
!             c  --- Spheroidal parameter
!             x  --- Argument of radial function ( x > 1.0 )
!             cv --- Characteristic value
!             KF --- Function code
!                    KF=1 for the first kind
!                    KF=2 for the second kind
!                    KF=3 for both the first and second kinds
!    Output:  R1F --- Radial function of the first kind
!             R1D --- Derivative of the radial function of the first kind
!             R2F --- Radial function of the second kind
!             R2D --- Derivative of the radial function of the second kind
!    Routines called:
!         (1) SDMN for computing expansion coefficients dk
!         (2) RMN1 for computing prolate and oblate radial
!             functions of the first kind
!         (3) RMN2L for computing prolate and oblate radial
!             functions of the second kind for a large argument
!         (4) RMN2SP for computing the prolate radial function
!             of the second kind for a small argument
!    ==============================================================

INTEGER, INTENT(IN)        :: m
INTEGER, INTENT(IN)        :: n
REAL (dp), INTENT(IN OUT)  :: c
REAL (dp), INTENT(IN)      :: x
REAL (dp), INTENT(IN)      :: cv
INTEGER, INTENT(IN)        :: kf
REAL (dp), INTENT(OUT)     :: r1f
REAL (dp), INTENT(OUT)     :: r1d
REAL (dp), INTENT(OUT)     :: r2f
REAL (dp), INTENT(OUT)     :: r2d

REAL (dp)  :: df(200)
INTEGER    :: id, kd

kd = 1
CALL sdmn(m, n, c, cv, kd, df)
IF (kf /= 2) THEN
  CALL rmn1(m, n, c, x, df, kd, r1f, r1d)
END IF
IF (kf > 1) THEN
  CALL rmn2l(m, n, c, x, df, kd, r2f, r2d, id)
  IF (id > -8) THEN
    CALL rmn2sp(m, n, c, x, cv, df, kd, r2f, r2d)
  END IF
END IF
RETURN
END SUBROUTINE rswfp



SUBROUTINE sdmn(m, n, c, cv, kd, df)

!    =====================================================
!    Purpose: Compute the expansion coefficients of the
!             prolate and oblate spheroidal functions, dk
!    Input :  m  --- Mode parameter
!             n  --- Mode parameter
!             c  --- Spheroidal parameter
!             cv --- Characteristic value
!             KD --- Function code
!                    KD=1 for prolate; KD=-1 for oblate
!    Output:  DF(k) --- Expansion coefficients dk;
!                       DF(1), DF(2), ... correspond to
!                       d0, d2, ... for even n-m and d1,
!                       d3, ... for odd n-m
!    =====================================================

INTEGER, INTENT(IN)     :: m
INTEGER, INTENT(IN)     :: n
REAL (dp), INTENT(IN)   :: c
REAL (dp), INTENT(IN)   :: cv
INTEGER, INTENT(IN)     :: kd
REAL (dp), INTENT(OUT)  :: df(200)

REAL (dp)  :: a(200), d(200), g(200)
REAL (dp)  :: cs, d2k, dk0, dk1, dk2, f, f0, f1, f2, fl, fs,  &
              r1, r3, r4, s0, su1, su2, sw
INTEGER    :: i, ip, j, k, k1, kb, nm

nm = 25 + INT(0.5*(n-m) + c)
IF (c < 1.0D-10) THEN
  df(1:nm) = 0D0
  df((n-m)/2+1) = 1.0D0
  RETURN
END IF
cs = c * c * kd
ip = 1
IF (n-m == 2*INT((n-m)/2)) ip = 0
DO  i = 1, nm + 2
  IF (ip == 0) k = 2 * (i-1)
  IF (ip == 1) k = 2 * i - 1
  dk0 = m + k
  dk1 = m + k + 1
  dk2 = 2 * (m+k)
  d2k = 2 * m + k
  a(i) = (d2k+2.0) * (d2k+1.0) / ((dk2+3.0)*(dk2+5.0)) * cs
  d(i) = dk0 * dk1 + (2.0*dk0*dk1-2.0*m*m-1.0) / ((dk2-1.0)*(dk2+ 3.0)) * cs
  g(i) = k * (k-1.0) / ((dk2-3.0)*(dk2-1.0)) * cs
END DO
fs = 1.0D0
f1 = 0.0D0
f0 = 1.0D-100
kb = 0
df(nm+1) = 0.0D0
DO  k = nm, 1, -1
  f = -((d(k+1)-cv)*f0 + a(k+1)*f1) / g(k+1)
  IF (ABS(f) > ABS(df(k+1))) THEN
    df(k) = f
    fl = df(k+1)
    f1 = f0
    f0 = f
    IF (ABS(f) > 1.0D+100) THEN
      DO  k1 = k, nm
        df(k1) = df(k1) * 1.0D-100
      END DO
      f1 = f1 * 1.0D-100
      f0 = f0 * 1.0D-100
    END IF
  ELSE
    kb = k
    fl = df(k+1)
    f1 = 1.0D-100
    f2 = -(d(1)-cv) / a(1) * f1
    df(1) = f1
    IF (kb == 1) THEN
      fs = f2
    ELSE IF (kb == 2) THEN
      df(2) = f2
      fs = -((d(2)-cv)*f2+g(2)*f1) / a(2)
    ELSE
      df(2) = f2
      DO  j = 3, kb + 1
        f = -((d(j-1)-cv)*f2+g(j-1)*f1) / a(j-1)
        IF (j <= kb) df(j) = f
        IF (ABS(f) > 1.0D+100) THEN
          DO  k1 = 1, j
            df(k1) = df(k1) * 1.0D-100
          END DO
          f = f * 1.0D-100
          f2 = f2 * 1.0D-100
        END IF
        f1 = f2
        f2 = f
      END DO
      fs = f
    END IF
    EXIT
  END IF
END DO

su1 = 0.0D0
r1 = 1.0D0
DO  j = m + ip + 1, 2 * (m+ip)
  r1 = r1 * j
END DO
su1 = df(1) * r1
DO  k = 2, kb
  r1 = -r1 * (k+m+ip-1.5D0) / (k-1.0D0)
  su1 = su1 + r1 * df(k)
END DO
su2 = 0.0D0
sw = su2
DO  k = kb + 1, nm
  IF (k /= 1) r1 = -r1 * (k+m+ip-1.5D0) / (k-1.0D0)
  su2 = su2 + r1 * df(k)
  IF (ABS(sw-su2) < ABS(su2)*1.0D-14) EXIT
  sw = su2
END DO

r3 = 1.0D0
DO  j = 1, (m+n+ip) / 2
  r3 = r3 * (j+0.5D0*(n+m+ip))
END DO
r4 = 1.0D0
DO  j = 1, (n-m-ip) / 2
  r4 = -4.0D0 * r4 * j
END DO
s0 = r3 / (fl*(su1/fs) + su2) / r4
DO  k = 1, kb
  df(k) = fl / fs * s0 * df(k)
END DO
DO  k = kb + 1, nm
  df(k) = s0 * df(k)
END DO
RETURN
END SUBROUTINE sdmn



SUBROUTINE rmn1(m, n, c, x, df, kd, r1f, r1d)

!    =======================================================
!    Purpose: Compute prolate and oblate spheroidal radial
!             functions of the first kind for given m, n, c and x
!    Routines called:
!         (1) SCKB for computing expansion coefficients c2k
!         (2) SPHJ for computing the spherical Bessel
!             functions of the first kind
!    =======================================================

INTEGER, INTENT(IN)        :: m
INTEGER, INTENT(IN)        :: n
REAL (dp), INTENT(IN OUT)  :: c
REAL (dp), INTENT(IN)      :: x
REAL (dp), INTENT(IN)      :: df(200)
INTEGER, INTENT(IN)        :: kd
REAL (dp), INTENT(OUT)     :: r1f
REAL (dp), INTENT(OUT)     :: r1d

REAL (dp)  :: ck(200), sj(0:251), dj(0:251)
REAL (dp)  :: a0, b0, cx, eps, r, r0, r1, r2, r3, reg, sa0, suc, sud, sum, &
              sw, sw1
INTEGER    :: ip, j, k, l, lg, nm, nm1, nm2, nm3, np

eps = 1.0D-14
ip = 1
nm1 = INT((n-m)/2)
IF (n-m == 2*nm1) ip = 0
nm = 25 + nm1 + INT(c)
reg = 1.0D0
IF (m+nm > 80) reg = 1.0D-200
r0 = reg
DO  j = 1, 2 * m + ip
  r0 = r0 * j
END DO
r = r0
suc = r * df(1)
sw = suc
DO  k = 2, nm
  r = r * (m+k-1.0) * (m+k+ip-1.5D0) / (k-1.0D0) / (k+ip-1.5D0)
  suc = suc + r * df(k)
  IF (k > nm1 .AND. ABS(suc-sw) < ABS(suc)*eps) EXIT
  sw = suc
END DO

IF (x == 0.0) THEN
  CALL sckb(m, n, c, df, ck)
  sum = 0.0D0
  DO  j = 1, nm
    sum = sum + ck(j)
    IF (ABS(sum-sw1) < ABS(sum)*eps) EXIT
    sw1 = sum
  END DO
  r1 = 1.0D0
  DO  j = 1, (n+m+ip) / 2
    r1 = r1 * (j+0.5D0*(n+m+ip))
  END DO
  r2 = 1.0D0
  DO  j = 1, m
    r2 = 2.0D0 * c * r2 * j
  END DO
  r3 = 1.0D0
  DO  j = 1, (n-m-ip) / 2
    r3 = r3 * j
  END DO
  sa0 = (2*(m+ip)+1) * r1 / (2.0**n*c**ip*r2*r3)
  IF (ip == 0) THEN
    r1f = sum / (sa0*suc) * df(1) * reg
    r1d = 0.0D0
  ELSE IF (ip == 1) THEN
    r1f = 0.0D0
    r1d = sum / (sa0*suc) * df(1) * reg
  END IF
  RETURN
END IF
cx = c * x
nm2 = 2 * nm + m
CALL sphj(nm2, cx, nm3, sj, dj)
a0 = (1.0D0-kd/(x*x)) ** (0.5D0*m) / suc
r1f = 0.0D0
DO  k = 1, nm
  l = 2 * k + m - n - 2 + ip
  IF (l == 4*INT(l/4)) lg = 1
  IF (l /= 4*INT(l/4)) lg = -1
  IF (k == 1) THEN
    r = r0
  ELSE
    r = r * (m+k-1.0) * (m+k+ip-1.5D0) / (k-1.0D0) / (k+ip-1.5D0)
  END IF
  np = m + 2 * k - 2 + ip
  r1f = r1f + lg * r * df(k) * sj(np)
  IF (k > nm1 .AND. ABS(r1f-sw) < ABS(r1f)*eps) EXIT
  sw = r1f
END DO

r1f = r1f * a0
b0 = kd * m / x ** 3.0D0 / (1.0 - kd/(x*x)) * r1f
sud = 0.0D0
DO  k = 1, nm
  l = 2 * k + m - n - 2 + ip
  IF (l == 4*INT(l/4)) lg = 1
  IF (l /= 4*INT(l/4)) lg = -1
  IF (k == 1) THEN
    r = r0
  ELSE
    r = r * (m+k-1.0) * (m+k+ip-1.5D0) / (k-1.0D0) / (k+ip-1.5D0)
  END IF
  np = m + 2 * k - 2 + ip
  sud = sud + lg * r * df(k) * dj(np)
  IF (k > nm1 .AND. ABS(sud-sw) < ABS(sud)*eps) EXIT
  sw = sud
END DO

r1d = b0 + a0 * c * sud
RETURN
END SUBROUTINE rmn1



SUBROUTINE sckb(m, n, c, df, ck)

!       ======================================================
!       Purpose: Compute the expansion coefficients of the
!                prolate and oblate spheroidal functions, c2k
!       Input :  m  --- Mode parameter
!                n  --- Mode parameter
!                c  --- Spheroidal parameter
!                DF(k) --- Expansion coefficients dk
!       Output:  CK(k) --- Expansion coefficients ck;
!                          CK(1), CK(2), ... correspond to c0, c2, ...
!       ======================================================


INTEGER, INTENT(IN)        :: m
INTEGER, INTENT(IN)        :: n
REAL (dp), INTENT(IN OUT)  :: c
REAL (dp), INTENT(IN)      :: df(200)
REAL (dp), INTENT(OUT)     :: ck(200)

REAL (dp)  :: d1, d2, d3, fac, r, r1, reg, sum, sw
INTEGER    :: i, i1, i2, ip, k, nm

IF (c <= 1.0D-10) c = 1.0D-10
nm = 25 + INT(0.5*(n-m)+c)
ip = 1
IF (n-m == 2*INT((n-m)/2)) ip = 0
reg = 1.0D0
IF (m+nm > 80) reg = 1.0D-200
fac = -0.5D0 ** m
DO  k = 0, nm - 1
  fac = -fac
  i1 = 2 * k + ip + 1
  r = reg
  DO  i = i1, i1 + 2 * m - 1
    r = r * i
  END DO
  i2 = k + m + ip
  DO  i = i2, i2 + k - 1
    r = r * (i+0.5D0)
  END DO
  sum = r * df(k+1)
  DO  i = k + 1, nm
    d1 = 2.0D0 * i + ip
    d2 = 2.0D0 * m + d1
    d3 = i + m + ip - 0.5D0
    r = r * d2 * (d2-1.0D0) * i * (d3+k) / (d1*(d1-1.0D0)*(i-k)*d3 )
    sum = sum + r * df(i+1)
    IF (ABS(sw-sum) < ABS(sum)*1.0D-14) EXIT
    sw = sum
  END DO

  r1 = reg
  DO  i = 2, m + k
    r1 = r1 * i
  END DO
  ck(k+1) = fac * sum / r1
END DO
RETURN
END SUBROUTINE sckb



SUBROUTINE sphj(n, x, nm, sj, dj)

!       =======================================================
!       Purpose: Compute spherical Bessel functions jn(x) and their derivatives
!       Input :  x --- Argument of jn(x)
!                n --- Order of jn(x)  ( n = 0,1,