求助-复变量的贝塞尔第一类和第二类函数的Fortran子程序
求助群里大神们，谁那里有可以求解复变量的贝塞尔第一类和第二类函数的Fortran子程序？肯求提供，感激涕零！去网上搜一搜，应该有很多的。 风平老涡 发表于 2021-1-26 00:34
去网上搜一搜，应该有很多的。
翻墙搜了很多，但都不能运行，编译不能通过。我也查了很多书带的子程序，但是可能代码太老也是没法编译 本帖最后由 风平老涡 于 2021-2-9 00:25 编辑
jlg1206 发表于 2021-1-31 18:57
翻墙搜了很多，但都不能运行，编译不能通过。我也查了很多书带的子程序，但是可能代码太老也是没法编译 ...
由于字数的限置，我把整个程序分成几个回帖。 jlg1206 发表于 2021-1-31 18:57
翻墙搜了很多，但都不能运行，编译不能通过。我也查了很多书带的子程序，但是可能代码太老也是没法编译 ...
MODULE Complex_Bessel
! SUBROUTINE CBSSLJ evaluates Bessel functions of complex ORDER
! SUBROUTINE CBSSLJ(x, nu, w)
! x = argument, nu = order, w = function value
! Latest revision - 20 July 2001
! amiller @ bigpond.net.au
! WARNING (from the NSWC Math. Library Manual)
! Precision:
! The real and imaginary parts are normally accurate to 11-12 significant
! digits when the function value is not near zero. The exception is when
! Real(z) is small, Imag(nu) is small, and 14 < |z| < 17.5 + 0.5|nu|^2
! In this region, the real part of the result is still accurate, but all
! accuracy may be lost in the imaginary part.
USE constants_NSWC
IMPLICIT NONE
CONTAINS
FUNCTION cpabs(x, y) RESULT(fn_val)
! --------------------------------------
! EVALUATION OF SQRT(X*X + Y*Y)
! --------------------------------------
REAL (dp), INTENT(IN):: x, y
REAL (dp) :: fn_val
! Local variable
REAL (dp):: a
IF (ABS(x) > ABS(y)) THEN
a = y / x
fn_val = ABS(x) * SQRT(1.0_dp + a*a)
RETURN
END IF
IF (y /= 0.0_dp) THEN
a = x / y
fn_val = ABS(y) * SQRT(1.0_dp + a*a)
RETURN
END IF
fn_val = 0.0_dp
RETURN
END FUNCTION cpabs
SUBROUTINE dcrec(x, y, u, v)
!-----------------------------------------------------------------------
! COMPLEX RECIPROCAL U + I*V = 1/(X + I*Y)
!-----------------------------------------------------------------------
REAL (dp), INTENT(IN) :: x, y
REAL (dp), INTENT(OUT):: u, v
! Local variables
REAL (dp):: t, d
IF (ABS(x) <= ABS(y)) THEN
t = x / y
d = y + t * x
u = t / d
v = -1.0_dp / d
RETURN
END IF
t = y / x
d = x + t * y
u = 1.0_dp / d
v = -t / d
RETURN
END SUBROUTINE dcrec
FUNCTION cdiv(a,b) RESULT(fn_val)
!-----------------------------------------------------------------------
! COMPLEX DIVISION A/B WHERE B IS NONZERO
!-----------------------------------------------------------------------
COMPLEX (dp), INTENT(IN):: a, b
COMPLEX (dp) :: fn_val
REAL (dp):: ai, ar, bi, br, d, t
REAL (dp):: u, v
ar = REAL(a, KIND=dp)
ai = AIMAG(a)
br = REAL(b, KIND=dp)
bi = AIMAG(b)
IF (ABS(br) >= ABS(bi)) THEN
t = bi / br
d = br + t * bi
u = (ar+ai*t) / d
v = (ai-ar*t) / d
fn_val = CMPLX(u,v, KIND=dp)
RETURN
END IF
t = br / bi
d = bi + t * br
u = (ar*t+ai) / d
v = (ai*t-ar) / d
fn_val = CMPLX(u,v, KIND=dp)
RETURN
END FUNCTION cdiv
FUNCTION dsin1(x) RESULT(fn_val)
!-----------------------------------------------------------------------
! REAL (dp) EVALUATION OF SIN(PI*X)
! --------------
! THE EXPANSION FOR SIN(PI*A) (ABS(A) <= PI/4) USING A1,...,A13
! IS ACCURATE TO WITHIN 2 UNITS OF THE 40-TH SIGNIFICANT DIGIT, AND
! THE EXPANSION FOR COS(PI*A) (ABS(A) <= PI/4) USING B1,...,B13
! IS ACCURATE TO WITHIN 4 UNITS OF THE 40-TH SIGNIFICANT DIGIT.
!-----------------------------------------------------------------------
REAL (dp), INTENT(IN):: x
REAL (dp) :: fn_val
! Local variables
REAL (dp) :: a, t, w
REAL (dp), PARAMETER :: pi = 3.141592653589793238462643383279502884197D+00
REAL (dp), PARAMETER :: a1 = -.1028083791780141522795259479153765743002D+00, &
a2= .3170868848763100170457042079710451905600D-02, &
a3= -.4657026956105571623449026167864697920000D-04,&
a4= .3989844942879455643410226655783424000000D-06, &
a5= -.2237397227721999776371894030796800000000D-08,&
a6= .8847045483056962709715066675200000000000D-11, &
a7= -.2598715447506450292885585920000000000000D-13,&
a8= .5893449774331011070033920000000000000000D-16 ,&
a9= -.1062975472045522550784000000000000000000D-18, &
a10 = .1561182648301780992000000000000000000000D-21, &
a11 = -.1903193516670976000000000000000000000000D-24, &
a12 = .1956617650176000000000000000000000000000D-27, &
a13 = -.1711276032000000000000000000000000000000D-30
REAL (dp), PARAMETER :: b1 = -.3084251375340424568385778437461297229882D+00, &
b2= .1585434424381550085228521039855226435920D-01, &
b3= -.3259918869273900136414318317506279360000D-03,&
b4= .3590860448591510079069203991239232000000D-05, &
b5= -.2461136950494199754009084061808640000000D-07,&
b6= .1150115912797405152263195572224000000000D-09, &
b7= -.3898073171259675439899172864000000000000D-12,&
b8= .1001886461636271969091584000000000000000D-14, &
b9= -.2019653396886572027084800000000000000000D-17,&
b10 = .3278483561466560512000000000000000000000D-20, &
b11 = -.4377345082051788800000000000000000000000D-23,&
b12 = .4891532381388800000000000000000000000000D-26, &
b13 = -.4617089843200000000000000000000000000000D-29
INTEGER:: max, n
!------------------------
! ****** MAX IS A MACHINE DEPENDENT CONSTANT. MAX IS THE
! LARGEST POSITIVE INTEGER THAT MAY BE USED.
! MAX = IPMPAR(3)
max = HUGE(3)
!------------------------
a = ABS(x)
t = MAX
IF (a >= t) THEN
fn_val = 0.0_dp
RETURN
END IF
n = a
t = n
a = a - t
IF (a <= 0.75_dp) THEN
IF (a < 0.25_dp) GO TO 10
! 0.25 <= A <= 0.75
a = 0.25_dp + (0.25_dp-a)
t = 16._dp * a * a
fn_val = (((((((((((((b13*t + b12)*t + b11)*t + b10)*t + b9)*t + b8)*t&
+ b7)*t + b6)*t + b5)*t + b4)*t + b3)*t + b2)*t + b1)*t +&
0.5_dp) + 0.5_dp
GO TO 20
END IF
! A < 0.25ORA > 0.75
a = 0.25_dp + (0.75_dp-a)
10 t = 16._dp * a * a
w = (((((((((((((a13*t + a12)*t + a11)*t + a10)*t + a9)*t + a8)*t + a7)*t&
+ a6)*t + a5)*t + a4)*t + a3)*t + a2)*t + a1)*t + 0.5_dp) + 0.5_dp
fn_val = pi * a * w
! TERMINATION
20 IF (x < 0.0) fn_val = -fn_val
IF (MOD(n,2) /= 0) fn_val = -fn_val
RETURN
END FUNCTION dsin1
FUNCTION dcos1 (x) RESULT(fn_val)
!-----------------------------------------------------------------------
! REAL (dp) EVALUATION OF COS(PI*X)
! --------------
! THE EXPANSION FOR SIN(PI*A) (ABS(A) .LE. PI/4) USING A1,...,A13
! IS ACCURATE TO WITHIN 2 UNITS OF THE 40-TH SIGNIFICANT DIGIT, AND
! THE EXPANSION FOR COS(PI*A) (ABS(A) .LE. PI/4) USING B1,...,B13
! IS ACCURATE TO WITHIN 4 UNITS OF THE 40-TH SIGNIFICANT DIGIT.
!-----------------------------------------------------------------------
REAL (dp), INTENT(IN):: x
REAL (dp) :: fn_val
REAL (dp):: a, t, w
INTEGER :: MAX, n
!------------------------
REAL (dp), PARAMETER:: pi = 3.141592653589793238462643383279502884197_dp
!------------------------
REAL (dp), PARAMETER:: &
a1= -.1028083791780141522795259479153765743002D+00,&
a2=.3170868848763100170457042079710451905600D-02,&
a3= -.4657026956105571623449026167864697920000D-04,&
a4=.3989844942879455643410226655783424000000D-06,&
a5= -.2237397227721999776371894030796800000000D-08,&
a6=.8847045483056962709715066675200000000000D-11,&
a7= -.2598715447506450292885585920000000000000D-13,&
a8=.5893449774331011070033920000000000000000D-16,&
a9= -.1062975472045522550784000000000000000000D-18,&
a10 =.1561182648301780992000000000000000000000D-21,&
a11 = -.1903193516670976000000000000000000000000D-24,&
a12 =.1956617650176000000000000000000000000000D-27,&
a13 = -.1711276032000000000000000000000000000000D-30
!------------------------
REAL (dp), PARAMETER:: &
b1= -.3084251375340424568385778437461297229882D+00,&
b2=.1585434424381550085228521039855226435920D-01,&
b3= -.3259918869273900136414318317506279360000D-03,&
b4=.3590860448591510079069203991239232000000D-05,&
b5= -.2461136950494199754009084061808640000000D-07,&
b6=.1150115912797405152263195572224000000000D-09,&
b7= -.3898073171259675439899172864000000000000D-12,&
b8=.1001886461636271969091584000000000000000D-14,&
b9= -.2019653396886572027084800000000000000000D-17,&
b10 =.3278483561466560512000000000000000000000D-20,&
b11 = -.4377345082051788800000000000000000000000D-23,&
b12 =.4891532381388800000000000000000000000000D-26,&
b13 = -.4617089843200000000000000000000000000000D-29
!------------------------
! ****** MAX IS A MACHINE DEPENDENT CONSTANT. MAX IS THE
! LARGEST POSITIVE INTEGER THAT MAY BE USED.
MAX = HUGE(0)
!------------------------
a = ABS(x)
t = MAX
IF (a < t) GO TO 10
fn_val = 1.d0
RETURN
10 n = a
t = n
a = a - t
IF (a > 0.75D0) GO TO 20
IF (a < 0.25D0) GO TO 21
! 0.25 .LE. A .LE. 0.75
a = 0.25D0 + (0.25D0 - a)
t = 16.d0*a*a
w = (((((((((((((a13*t + a12)*t + a11)*t + a10)*t + a9)*t +&
a8)*t + a7)*t + a6)*t + a5)*t + a4)*t + a3)*t + a2)*t + a1)*t + 0.5D0) + 0.5D0
fn_val = pi*a*w
GO TO 30
! A .LT. 0.25ORA .GT. 0.75
20 a = 0.25D0 + (0.75D0 - a)
n = n - 1
21 t = 16.d0*a*a
fn_val = (((((((((((((b13*t + b12)*t + b11)*t + b10)*t + b9)*t + b8)*t + &
b7)*t + b6)*t + b5)*t + b4)*t + b3)*t + b2)*t + b1)*t + 0.5D0) + 0.5D0
! TERMINATION
30 IF (MOD(n,2) /= 0) fn_val = -fn_val
RETURN
END FUNCTION dcos1
FUNCTION drexp(x) RESULT(fn_val)
!-----------------------------------------------------------------------
! EVALUATION OF THE FUNCTION EXP(X) - 1
!-----------------------------------------------------------------------
REAL (dp), INTENT(IN):: x
REAL (dp) :: fn_val
! Local variables
REAL (dp) :: e, w, z
REAL (dp) :: a0 = .248015873015873015873016D-04, &
a1 = -.344452080605731005808147D-05, a2 = .206664230430046597475413D-06,&
a3 = -.447300111094328162971036D-08, a4 = .114734027080634968083920D-11,&
b1 = -.249994190011341852652396D+00, b2 = .249987228833107957725728D-01,&
b3 = -.119037506846942249362528D-02, b4 = .228908693387350391768682D-04
REAL (dp) :: c1 = .1666666666666666666666666666666667D+00, &
c2 = .4166666666666666666666666666666667D-01, &
c3 = .8333333333333333333333333333333333D-02, &
c4 = .1388888888888888888888888888888889D-02, &
c5 = .1984126984126984126984126984126984D-03
!---------------------------
IF (ABS(x) <= 0.15_dp) THEN
! Z IS A MINIMAX APPROXIMATION OF THE SERIES
! C6 + C7*X + C8*X**2 + ....
! THIS APPROXIMATION IS ACCURATE TO WITHIN
! 1 UNIT OF THE 23-RD SIGNIFICANT DIGIT.
! THE RESULTING VALUE FOR W IS ACCURATE TO
! WITHIN 1 UNIT OF THE 33-RD SIGNIFICANT DIGIT.
z = ((((a4*x + a3)*x + a2)*x + a1)*x + a0) /&
((((b4*x + b3)*x + b2)*x + b1)*x + 1._dp)
w = ((((((z*x + c5)*x + c4)*x + c3)*x + c2)*x + c1)*x + 0.5_dp)*x + 1._dp
fn_val = x * w
RETURN
END IF
IF (x >= 0.0_dp) THEN
e = EXP(x)
fn_val = e * (0.5_dp + (0.5_dp - 1.0_dp/e))
RETURN
END IF
IF (x >= -77._dp) THEN
fn_val = (EXP(x)-0.5_dp) - 0.5_dp
RETURN
END IF
fn_val = -1._dp
RETURN
END FUNCTION drexp
SUBROUTINE dcgama(mo, z, w)
!-----------------------------------------------------------------------
! EVALUATION OF THE COMPLEX GAMMA AND LOGGAMMA FUNCTIONS
! ---------------
! MO IS AN INTEGER.Z AND W ARE INTERPRETED AS REAL (dp)
! COMPLEX NUMBERS.IT IS ASSUMED THAT Z(1) AND Z(2) ARE THE REAL
! AND IMAGINARY PARTS OF THE COMPLEX NUMBER Z, AND THAT W(1) AND
! W(2) ARE THE REAL AND IMAGINARY PARTS OF W.
! W = GAMMA(Z) IF MO = 0
! W = LN(GAMMA(Z)) OTHERWISE
!-----------------------------------------------------------------------
INTEGER, INTENT(IN) :: mo
COMPLEX (dp), INTENT(IN) :: z
COMPLEX (dp), INTENT(OUT):: w
! Local variables
REAL (dp), PARAMETER :: c0(30)&
= (/ .8333333333333333333333333333333333333333D-01,&
-.2777777777777777777777777777777777777778D-02,&
.7936507936507936507936507936507936507937D-03,&
-.5952380952380952380952380952380952380952D-03,&
.8417508417508417508417508417508417508418D-03,&
-.1917526917526917526917526917526917526918D-02,&
.6410256410256410256410256410256410256410D-02,&
-.2955065359477124183006535947712418300654D-01,&
.1796443723688305731649384900158893966944D+00,&
-.1392432216905901116427432216905901116427D+01,&
.1340286404416839199447895100069013112491D+02,&
-.1568482846260020173063651324520889738281D+03,&
.2193103333333333333333333333333333333333D+04,&
-.3610877125372498935717326521924223073648D+05,&
.6914722688513130671083952507756734675533D+06,&
-.1523822153940741619228336495888678051866D+08,&
.3829007513914141414141414141414141414141D+09,&
-.1088226603578439108901514916552510537473D+11,&
.3473202837650022522522522522522522522523D+12,&
-.1236960214226927445425171034927132488108D+14,&
.4887880647930793350758151625180229021085D+15,&
-.2132033396091937389697505898213683855747D+17,&
.1021775296525700077565287628053585500394D+19,&
-.5357547217330020361082770919196920448485D+20,&
.3061578263704883415043151051329622758194D+22,&
-.1899991742639920405029371429306942902947D+24,&
.1276337403382883414923495137769782597654D+26,&
-.9252847176120416307230242348347622779519D+27,&
.7218822595185610297836050187301637922490D+29,&
-.6045183405995856967743148238754547286066D+31 /),&
dlpi = 1.144729885849400174143427351353058711647_dp,&
hl2p =.9189385332046727417803297364056176398614_dp,&
pi = 3.141592653589793238462643383279502884197_dp,&
pi2 = 6.283185307179586476925286766559005768394_dp
REAL (dp) :: a, a1, a2, c, cn, cut, d, eps, et, e2t, h1, h2, q1, q2, s, sn,&
s1, s2, t, t1, t2, u, u1, u2, v1, v2, w1, w2, x, y, y2
INTEGER :: j, k, max, n, nm1
!---------------------------
! DLPI = LOG(PI)
! HL2P = 0.5 * LOG(2*PI)
!---------------------------
! ****** MAX AND EPS ARE MACHINE DEPENDENT CONSTANTS.
! MAX IS THE LARGEST POSITIVE INTEGER THAT MAY
! BE USED, AND EPS IS THE SMALLEST NUMBER SUCH
! THAT1._dp + EPS > 1._dp.
! MAX = IPMPAR(3)
max = HUGE(3)
eps = EPSILON(1.0_dp)
!---------------------------
x = REAL(z, KIND=dp)
y = AIMAG(z)
IF (x < 0._dp) THEN
!-----------------------------------------------------------------------
! CASE WHEN THE REAL PART OF Z IS NEGATIVE
!-----------------------------------------------------------------------
y = ABS(y)
t = -pi * y
et = EXP(t)
e2t = et * et
! SETA1 = (1 + E2T)/2ANDA2 = (1 - E2T)/2
a1 = 0.5_dp * (1._dp+e2t)
t2 = t + t
IF (t2 >= -0.15_dp) THEN
a2 = -0.5_dp * drexp(t2)
ELSE
a2 = 0.5_dp * (0.5_dp+(0.5_dp-e2t))
END IF
! COMPUTE SIN(PI*X) AND COS(PI*X)
u = MAX
IF (ABS(x) >= MIN(u,1._dp/eps)) GO TO 80
k = ABS(x)
u = x + k
k = MOD(k,2)
IF (u <= -0.5_dp) THEN
u = 0.5_dp + (0.5_dp+u)
k = k + 1
END IF
u = pi * u
sn = SIN(u)
cn = COS(u)
IF (k == 1) THEN
sn = -sn
cn = -cn
END IF
! SETH1 + H2*ITOPI/SIN(PI*Z)ORLOG(PI/SIN(PI*Z))
a1 = sn * a1
a2 = cn * a2
a = a1 * a1 + a2 * a2
IF (a == 0._dp) GO TO 80
IF (mo == 0) THEN
h1 = a1 / a
h2 = -a2 / a
c = pi * et
h1 = c * h1
h2 = c * h2
ELSE
h1 = (dlpi+t) - 0.5_dp * LOG(a)
h2 = -ATAN2(a2,a1)
END IF
IF (AIMAG(z) >= 0._dp) THEN
x = 1.0 - x
y = -y
ELSE
h2 = -h2
x = 1.0 - x
END IF
END IF
!-----------------------------------------------------------------------
! CASE WHEN THE REAL PART OF Z IS NONNEGATIVE
!-----------------------------------------------------------------------
w1 = 0._dp
w2 = 0._dp
n = 0
t = x
y2 = y * y
a = t * t + y2
cut = 225._dp
IF (eps > 1.d-30) cut = 144._dp
IF (eps > 1.d-20) cut = 64._dp
IF (a < cut) THEN
IF (a == 0._dp) GO TO 80
10 n = n + 1
t = t + 1._dp
a = t * t + y2
IF (a < cut) GO TO 10
! LET S1 + S2*I BE THE PRODUCT OF THE TERMS (Z+J)/(Z+N)
u1 = (x*t+y2) / a
u2 = y / a
s1 = u1
s2 = n * u2
IF (n >= 2) THEN
u = t / a
nm1 = n - 1
DO j = 1, nm1
v1 = u1 + j * u
v2 = (n-j) * u2
c = s1 * v1 - s2 * v2
d = s1 * v2 + s2 * v1
s1 = c
s2 = d
END DO
END IF
! SETW1 + W2*I = LOG(S1 + S2*I)WHEN MO IS NONZERO
s = s1 * s1 + s2 * s2
IF (mo /= 0) THEN
w1 = 0.5_dp * LOG(s)
w2 = ATAN2(s2,s1)
END IF
END IF
! SETV1 + V2*I = (Z - 0.5) * LOG(Z + N) - Z
t1 = 0.5_dp * LOG(a) - 1._dp
t2 = ATAN2(y,t)
u = x - 0.5_dp
v1 = (u*t1-0.5_dp) - y * t2
v2 = u * t2 + y * t1
! LET A1 + A2*I BE THE ASYMPTOTIC SUM
u1 = t / a
u2 = -y / a
q1 = u1 * u1 - u2 * u2
q2 = 2._dp * u1 * u2
a1 = 0._dp
a2 = 0._dp
DO j = 1, 30
t1 = a1
t2 = a2
a1 = a1 + c0(j) * u1
a2 = a2 + c0(j) * u2
IF (a1 == t1) THEN
IF (a2 == t2) GO TO 40
END IF
t1 = u1 * q1 - u2 * q2
t2 = u1 * q2 + u2 * q1
u1 = t1
u2 = t2
END DO
!-----------------------------------------------------------------------
! GATHERING TOGETHER THE RESULTS
!-----------------------------------------------------------------------
40 w1 = (((a1+hl2p)-w1)+v1) - n
w2 = (a2-w2) + v2
IF (REAL(z, KIND=dp) < 0.0_dp) GO TO 60
IF (mo == 0) THEN
! CASE WHEN THE REAL PART OF Z IS NONNEGATIVE AND MO = 0
a = EXP(w1)
w1 = a * COS(w2)
w2 = a * SIN(w2)
IF (n == 0) GO TO 70
c = (s1*w1+s2*w2) / s
d = (s1*w2-s2*w1) / s
w1 = c
w2 = d
GO TO 70
END IF
! CASE WHEN THE REAL PART OF Z IS NONNEGATIVE AND MO IS NONZERO.
! THE ANGLE W2 IS REDUCED TO THE INTERVAL -PI < W2 <= PI.
50 IF (w2 <= pi) THEN
k = 0.5_dp - w2 / pi2
w2 = w2 + pi2 * k
GO TO 70
END IF
k = w2 / pi2 - 0.5_dp
u = k + 1
w2 = w2 - pi2 * u
IF (w2 <= -pi) w2 = pi
GO TO 70
! CASE WHEN THE REAL PART OF Z IS NEGATIVE AND MO IS NONZERO
60 IF (mo /= 0) THEN
w1 = h1 - w1
w2 = h2 - w2
GO TO 50
END IF
! CASE WHEN THE REAL PART OF Z IS NEGATIVE AND MO = 0
a = EXP(-w1)
t1 = a * COS(-w2)
t2 = a * SIN(-w2)
w1 = h1 * t1 - h2 * t2
w2 = h1 * t2 + h2 * t1
IF (n /= 0) THEN
c = w1 * s1 - w2 * s2
d = w1 * s2 + w2 * s1
w1 = c
w2 = d
END IF
! TERMINATION
70 w = CMPLX(w1, w2, KIND=dp)
RETURN
!-----------------------------------------------------------------------
! THE REQUESTED VALUE CANNOT BE COMPUTED
!-----------------------------------------------------------------------
80 w = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
RETURN
END SUBROUTINE dcgama
FUNCTION cgam0(z) RESULT(fn_val)
!-----------------------------------------------------------------------
! EVALUATION OF 1/GAMMA(1 + Z)FOR ABS(Z) < 1.0
!-----------------------------------------------------------------------
COMPLEX (dp), INTENT(IN):: z
COMPLEX (dp) :: fn_val
COMPLEX (dp):: w
INTEGER :: i, k, n
!-----------------------
REAL (dp):: x, y
REAL, PARAMETER :: a(25) = (/ .577215664901533_dp, -.655878071520254_dp,&
-.420026350340952D-01, .166538611382291_dp, -.421977345555443D-01,&
-.962197152787697D-02, .721894324666310D-02, -.116516759185907D-02,&
-.215241674114951D-03, .128050282388116D-03, -.201348547807882D-04,&
-.125049348214267D-0, .113302723198170D-05, -.205633841697761D-0,&
.611609510448142D-08, .500200764446922D-08, -.118127457048702D-08,&
.104342671169110D-09, .778226343990507D-11, -.369680561864221D-11,&
.510037028745448D-12, -.205832605356651D-13, -.534812253942302D-14,&
.122677862823826D-14, -.118125930169746D-15 /)
!-----------------------
n = 25
x = REAL(z, KIND=dp)
y = AIMAG(z)
IF (x*x+y*y <= 0.04D0) n = 14
k = n
w = a(n)
DOi = 2, n
k = k - 1
w = a(k) + z * w
END DO
fn_val = 1.0D0 + z * w
RETURN
END FUNCTION cgam0
FUNCTION dgamma(a) RESULT(fn_val)
!-----------------------------------------------------------------------
! EVALUATION OF THE GAMMA FUNCTION FOR
! REAL (dp) ARGUMENTS
! -----------
! DGAMMA(A) IS ASSIGNED THE VALUE 0 WHEN THE GAMMA FUNCTION CANNOT
! BE COMPUTED.
!-----------------------------------------------------------------------
! WRITTEN BY ALFRED H. MORRIS, JR.
! NAVAL SURFACE WEAPONS CENTER
! DAHLGREN, VIRGINIA
!-----------------------------------------------------------------------
REAL (dp), INTENT(IN) :: a
REAL (dp) :: fn_val
! Local variables
REAL (dp), PARAMETER :: d = 0.41893853320467274178032973640562_dp,&
pi = 3.14159265358979323846264338327950_dp
REAL (dp) :: s, t, x, w
INTEGER :: j, n
!-----------------------------------------------------------------------
! D = 0.5*(LN(2*PI) - 1)
!-----------------------------------------------------------------------
fn_val = 0.0_dp
x = a
IF (ABS(a) <= 20._dp) THEN
!-----------------------------------------------------------------------
! EVALUATION OF DGAMMA(A) FOR ABS(A) <= 20
!-----------------------------------------------------------------------
t = 1.0_dp
n = x
n = n - 1
! LET T BE THE PRODUCT OF A-J WHEN A >= 2
IF (n < 0) THEN
GO TO 40
ELSE IF (n == 0) THEN
GO TO 30
END IF
DO j = 1, n
x = x - 1._dp
t = x * t
END DO
30 x = x - 1._dp
GO TO 60
! LET T BE THE PRODUCT OF A+J WHEN A < 1
40 t = a
IF (a <= 0._dp) THEN
n = -n - 1
IF (n /= 0) THEN
DO j = 1, n
x = x + 1._dp
t = x * t
END DO
END IF
x = (x+0.5_dp) + 0.5_dp
t = x * t
IF (t == 0._dp) RETURN
END IF
! THE FOLLOWING CODE CHECKS IF 1/T CAN OVERFLOW. THIS
! CODE MAY BE OMITTED IF DESIRED.
IF (ABS(t) < 1.d-33) THEN
IF (ABS(t)*HUGE(1.0_dp) <= 1.000000001_dp) RETURN
fn_val = 1._dp / t
RETURN
END IF
! COMPUTE DGAMMA(1 + X) FOR 0 <= X < 1
60 fn_val = 1._dp / (1._dp + dgam1(x))
! TERMINATION
IF (a >= 1._dp) THEN
fn_val = fn_val * t
RETURN
END IF
fn_val = fn_val / t
RETURN
END IF
!-----------------------------------------------------------------------
! EVALUATION OF DGAMMA(A) FOR ABS(A) > 20
!-----------------------------------------------------------------------
IF (ABS(a) >= 1.d3) RETURN
IF (a <= 0._dp) THEN
s = dsin1(a) / pi
IF (s == 0._dp) RETURN
x = -a
END IF
! COMPUTE THE MODIFIED ASYMPTOTIC SUM
w = dpdel(x)
! FINAL ASSEMBLY
w = (d+w) + (x-0.5_dp) * (LOG(x)-1._dp)
IF (w > dxparg(0)) RETURN
fn_val = EXP(w)
IF (a < 0._dp) fn_val = (1._dp/(fn_val*s)) / x
RETURN
END FUNCTION dgamma
FUNCTION glog(x) RESULT(fn_val)
! -------------------
! EVALUATION OF LN(X) FOR X >= 15
! -------------------
REAL (dp), INTENT(IN) :: x
REAL (dp) :: fn_val
! Local variables
REAL (dp), PARAMETER :: c1 = .286228750476730_dp, c2 = .399999628131494+dp,&
c3 = .666666666752663_dp
REAL (dp), PARAMETER :: w(163) = (/ .270805020110221007D+01, .277258872223978124D+01,&
.283321334405621608D+01, .289037175789616469D+01, .294443897916644046D+01, .299573227355399099D+01, .304452243772342300D+01, &
.309104245335831585D+01, .313549421592914969D+01, .317805383034794562D+01, .321887582486820075D+01, .325809653802148205D+01, &
.329583686600432907D+01, .333220451017520392D+01, .336729582998647403D+01, .340119738166215538D+01, .343398720448514625D+01, &
.346573590279972655D+01, .349650756146648024D+01, .352636052461616139D+01, .355534806148941368D+01, .358351893845611000D+01, &
.361091791264422444D+01, .363758615972638577D+01, .366356164612964643D+01, .368887945411393630D+01, .371357206670430780D+01, &
.373766961828336831D+01, .376120011569356242D+01, .378418963391826116D+01, .380666248977031976D+01, .382864139648909500D+01, &
.385014760171005859D+01, .387120101090789093D+01, .389182029811062661D+01, .391202300542814606D+01, .393182563272432577D+01, &
.395124371858142735D+01, .397029191355212183D+01, .398898404656427438D+01, .400733318523247092D+01, .402535169073514923D+01, &
.404305126783455015D+01, .406044301054641934D+01, .407753744390571945D+01, .409434456222210068D+01, .411087386417331125D+01, &
.412713438504509156D+01, .414313472639153269D+01, .415888308335967186D+01, .417438726989563711D+01, .418965474202642554D+01, &
.420469261939096606D+01, .421950770517610670D+01, .423410650459725938D+01, .424849524204935899D+01, .426267987704131542D+01, &
.427666611901605531D+01, .429045944114839113D+01, .430406509320416975D+01, .431748811353631044D+01, .433073334028633108D+01, &
.434380542185368385D+01, .435670882668959174D+01, .436944785246702149D+01, .438202663467388161D+01, .439444915467243877D+01, &
.440671924726425311D+01, .441884060779659792D+01, .443081679884331362D+01, .444265125649031645D+01, .445434729625350773D+01, &
.446590811865458372D+01, .447733681447820647D+01, .448863636973213984D+01, .449980967033026507D+01, .451085950651685004D+01, &
.452178857704904031D+01, .453259949315325594D+01, .454329478227000390D+01, .455387689160054083D+01, .456434819146783624D+01, &
.457471097850338282D+01, .458496747867057192D+01, .459511985013458993D+01, .460517018598809137D+01, .461512051684125945D+01, &
.462497281328427108D+01, .463472898822963577D+01, .464439089914137266D+01, .465396035015752337D+01, .466343909411206714D+01, &
.467282883446190617D+01, .468213122712421969D+01, .469134788222914370D+01, .470048036579241623D+01, .470953020131233414D+01, &
.471849887129509454D+01, .472738781871234057D+01, .473619844839449546D+01, .474493212836325007D+01, .475359019110636465D+01, &
.476217393479775612D+01, .477068462446566476D+01, .477912349311152939D+01, .478749174278204599D+01, .479579054559674109D+01, &
.480402104473325656D+01, .481218435537241750D+01, .482028156560503686D+01, .482831373730230112D+01, .483628190695147800D+01, &
.484418708645859127D+01, .485203026391961717D+01, .485981240436167211D+01, .486753445045558242D+01, .487519732320115154D+01, &
.488280192258637085D+01, .489034912822175377D+01, .489783979995091137D+01, .490527477843842945D+01, .491265488573605201D+01, &
.491998092582812492D+01, .492725368515720469D+01, .493447393313069176D+01, .494164242260930430D+01, .494875989037816828D+01, &
.495582705760126073D+01, .496284463025990728D+01, .496981329957600062D+01, .497673374242057440D+01, .498360662170833644D+01, &
.499043258677873630D+01, .499721227376411506D+01, .500394630594545914D+01, .501063529409625575D+01, .501727983681492433D+01, &
.502388052084627639D+01, .503043792139243546D+01, .503695260241362916D+01, .504342511691924662D+01, .504985600724953705D+01, &
.505624580534830806D+01, .506259503302696680D+01, .506890420222023153D+01, .507517381523382692D+01, .508140436498446300D+01, &
.508759633523238407D+01, .509375020080676233D+01, .509986642782419842D+01, .510594547390058061D+01, .511198778835654323D+01, &
.511799381241675511D+01, .512396397940325892D+01, .512989871492307347D+01, .513579843705026176D+01, .514166355650265984D+01, &
.514749447681345304D+01, .515329159449777895D+01, .515905529921452903D+01, .516478597392351405D+01, .517048399503815178D+01, &
.517614973257382914D+01 /)
REAL (dp) :: t, t2, z
INTEGER :: n
! -------------------
IF (x < 178.0_dp) THEN
n = x
t = (x-n) / (x+n)
t2 = t * t
z = (((c1*t2 + c2)*t2 + c3)*t2 + 2.0) * t
fn_val = w(n-14) + z
RETURN
END IF
fn_val = LOG(x)
RETURN
END FUNCTION glog
SUBROUTINE cbsslj(z,cnu,w)
!-----------------------------------------------------------------------
! EVALUATION OF THE COMPLEX BESSEL FUNCTION J (Z)
! CNU
!-----------------------------------------------------------------------
! WRITTEN BY
! ANDREW H. VAN TUYL AND ALFRED H. MORRIS, JR.
! NAVAL SURFACE WARFARE CENTER
! OCTOBER, 1991
! A MODIFICATION OF THE PROCEDURE DEVELOPED BY ALLEN V. HERSHEY
! (NAVAL SURFACE WARFARE CENTER) IN 1978 FOR HANDLING THE DEBYE
! APPROXIMATION IS EMPLOYED.
!-----------------------------------------------------------------------
COMPLEX (dp), INTENT(IN) :: z
COMPLEX (dp), INTENT(IN) :: cnu
COMPLEX (dp), INTENT(OUT):: w
COMPLEX (dp):: c, nu, s, sm1, sm2, t, tsc, w0, w1, zn, zz
!-----------------------
REAL (dp) :: a, cn1, cn2, e, fn
REAL (dp) :: pn, qm, qn, qnp1
REAL (dp) :: r, rn2, r2, sn, t1, t2
REAL (dp) :: u, v, x, y
INTEGER :: i, k, m, n
REAL (dp), PARAMETER:: pi = 3.141592653589793238462643383279502884197_dp
!-----------------------
x = REAL(z, KIND=dp)
y = AIMAG(z)
r = cpabs(x,y)
cn1 = REAL(cnu, KIND=dp)
cn2 = AIMAG(cnu)
rn2 = cn1 * cn1 + cn2 * cn2
pn = INT(cn1)
fn = cn1 - pn
sn = 1.0_dp
! CALCULATION WHEN ORDER IS AN INTEGER
IF (fn == 0.0_dp .AND. cn2 == 0.0_dp) THEN
n = pn
pn = ABS(pn)
cn1 = pn
IF (n < 0 .AND. n /= (n/2)*2) sn = -1.0_dp
END IF
! SELECTION OF METHOD
IF (r > 17.5D0) THEN
IF (r > 17.5D0 + 0.5D0*rn2) GO TO 10
GO TO 20
END IF
! USE MACLAURIN EXPANSION AND RECURSION
IF (cn1 < 0.0D0) THEN
qn = -1.25D0 * (r + 0.5D0*ABS(cn2) - ABS(y-0.5D0*cn2))
IF (cn1 < qn) THEN
qn = 1.25D0 * (r - MAX(1.2D0*r,ABS(y-cn2)))
IF (cn1 < qn) THEN
qn = MIN(pn, REAL(-INT(1.25D0*(r-ABS(cn2))), KIND=dp))
GO TO 60
END IF
END IF
END IF
r2 = r * r
qm = 0.0625D0 * r2 * r2 - cn2 * cn2
qn = MAX(pn, REAL(INT(SQRT(MAX(0.0D0,qm))), KIND=dp))
GO TO 60
! USE ASYMPTOTIC EXPANSION
10 CALL cbja(z,cnu,w)
RETURN
! CALCULATION FOR 17.5 < ABS(Z) <= 17.5 + 0.5*ABS(CNU)**2
20 n = 0
IF (ABS(cn2) < 0.8D0*ABS(y)) THEN
qm = -1.25D0 * (r + 0.5D0*ABS(cn2) - ABS(y-0.5D0*cn2))
IF (cn1 < qm) THEN
qm = 1.25D0 * (r - MAX(1.2D0*r, ABS(y-cn2)))
IF (cn1 < qm) n = 1
END IF
END IF
qn = pn
a = 4.d-3 * r * r
zz = z
IF (x < 0.0D0) zz = -z
! CALCULATION OF ZONE OF EXCLUSION OF DEBYE APPROXIMATION
30 nu = CMPLX(qn+fn,cn2, KIND=dp)
zn = nu / z
t2 = AIMAG(zn) * AIMAG(zn)
u = 1.0D0 - REAL(zn, KIND=dp)
t1 = u * u + t2
u = 1.0D0 + DBLE(zn)
t2 = u * u + t2
u = t1 * t2
v = a * u / (t1*t1 + t2*t2)
IF (u*v*v <= 1.0D0) THEN
! THE ARGUMENT LIES INSIDE THE ZONE OF EXCLUSION
qn = qn + 1.0D0
IF (n == 0) GO TO 30
! USE MACLAURIN EXPANSION WITH FORWARD RECURRENCE
qn = MIN(pn, REAL(-INT(1.25D0*(r-ABS(cn2))), KIND=dp))
ELSE
! USE BACKWARD RECURRENCE STARTING FROM THE ASYMPTOTIC EXPANSION
qnp1 = qn + 1.0D0
IF (ABS(qn) < ABS(pn)) THEN
IF (r >= 17.5D0 + 0.5D0*(qnp1*qnp1 + cn2*cn2)) THEN
nu = CMPLX(qn+fn,cn2, KIND=dp)
CALL cbja(zz,nu,sm1)
nu = CMPLX(qnp1+fn,cn2, KIND=dp)
CALL cbja(zz,nu,sm2)
GO TO 40
END IF
END IF
! USE BACKWARD RECURRENCE STARTING FROM THE DEBYE APPROXIMATION
nu = CMPLX(qn+fn,cn2, KIND=dp)
CALL cbdb(zz,nu,fn,sm1)
IF (qn == pn) GO TO 50
nu = CMPLX(qnp1+fn,cn2, KIND=dp)
CALL cbdb(zz,nu,fn,sm2)
40 nu = CMPLX(qn+fn,cn2, KIND=dp)
tsc = 2.0D0 * nu * sm1 / zz - sm2
sm2 = sm1
sm1 = tsc
qn = qn - 1.0D0
IF (qn /= pn) GO TO 40
50 w = sm1
IF (sn < 0.0D0) w = -w
IF (x >= 0.0D0) RETURN
nu = pi * CMPLX(-cn2,cn1, KIND=dp)
IF (y < 0.0D0) nu = -nu
w = EXP(nu) * w
RETURN
END IF
! USE MACLAURIN EXPANSION WITH FORWARD OR BACKWARD RECURRENCE.
60 m = qn - pn
IF (ABS(m) <= 1) THEN
nu = CMPLX(cn1,cn2, KIND=dp)
CALL cbjm(z,nu,w)
ELSE
nu = CMPLX(qn+fn,cn2, KIND=dp)
CALL cbjm(z,nu,w1)
w0 = 0.25D0 * z * z
IF (m <= 0) THEN
! FORWARD RECURRENCE
m = ABS(m)
nu = nu + 1.0D0
CALL cbjm(z,nu,w)
DOi = 2, m
c = nu * (nu+1.0D0)
t = (c/w0) * (w-w1)
w1 = w
w = t
nu = nu + 1.0D0
END DO
ELSE
! BACKWARD RECURRENCE
nu = nu - 1.0D0
CALL cbjm(z,nu,w)
DOi = 2, m
c = nu * (nu+1.0D0)
t = (w0/c) * w1
w1 = w
w = w - t
nu = nu - 1.0D0
END DO
END IF
END IF
! FINAL ASSEMBLY
IF (fn == 0.0D0 .AND. cn2 == 0.0D0) THEN
k = pn
IF (k == 0) RETURN
e = sn / dgamma(pn+1.0D0)
w = e * w * (0.5D0*z) ** k
RETURN
END IF
s = cnu * LOG(0.5D0*z)
w = EXP(s) * w
IF (rn2 <= 0.81D0) THEN
w = w * cgam0(cnu)
RETURN
END IF
CALL dcgama(0,cnu,t)
w = cdiv(w,cnu*t)
RETURN
END SUBROUTINE cbsslj
SUBROUTINE cbjm(z,cnu,w)
!-----------------------------------------------------------------------
! COMPUTATION OF(Z/2)**(-CNU) * GAMMA(CNU + 1) * J(CNU,Z)
! -----------------
! THE MACLAURIN EXPANSION IS USED.IT IS ASSUMED THAT CNU IS NOT
! A NEGATIVE INTEGER.
!-----------------------------------------------------------------------
COMPLEX (dp), INTENT(IN) :: z
COMPLEX (dp), INTENT(IN) :: cnu
COMPLEX (dp), INTENT(OUT):: w
COMPLEX (dp):: nu, nup1, p, s, sn, t, ti
!--------------------------
REAL (dp):: a, a0, eps, inu, m, rnu
INTEGER :: i, imin, k, km1, km2
!--------------------------
! ****** EPS IS A MACHINE DEPENDENT CONSTANT. EPS IS THE
! SMALLEST NUMBER SUCH THAT 1.0 + EPS .GT. 1.0 .
eps = EPSILON(0.0_dp)
!--------------------------
s = -0.25D0 * (z*z)
nu = cnu
rnu = REAL(nu, KIND=dp)
inu = AIMAG(nu)
a = 0.5D0 + (0.5D0+rnu)
nup1 = CMPLX(a,inu, KIND=dp)
IF (a > 0.0D0) THEN
m = 1.0D0
t = s / nup1
w = 1.0D0 + t
ELSE
! ADD 1.0 AND THE FIRST K-1 TERMS
k = INT(-a) + 2
km1 = k - 1
w = (1.0D0,0.0D0)
t = w
DOi = 1, km1
m = i
t = t * (s/(m*(nu+m)))
w = w + t
IF (anorm(t) <= eps*anorm(w)) GO TO 20
END DO
GO TO 50
! CHECK IF THE (K-1)-ST AND K-TH TERMS CAN BE IGNORED.
! IF SO THEN THE SUMMATION IS COMPLETE.
20 IF (i /= km1) THEN
imin = i + 1
IF (imin < k-5) THEN
ti = t
m = km1
t = s / (nu+m)
a0 = anorm(t) / m
t = t * (s/(nu+(m+1.0D0)))
a = anorm(t) / (m*(m+1.0D0))
a = MAX(a,a0)
t = (1.0D0,0.0D0)
km2 = k - 2
DOi = imin, km2
m = i
t = t * (s/(m*(nu+m)))
IF (a*anorm(t) < 0.5D0) RETURN
END DO
t = t * ti
imin = km2
END IF
! ADD THE (K-1)-ST TERM
a = 1.0D0
p = (1.0D0,0.0D0)
sn = p
DOi = imin, km1
m = i
a = a * m
p = p * (nu+m)
sn = s * sn
END DO
t = t * (cdiv(sn,p)/a)
w = w + t
END IF
END IF
! ADD THE REMAINING TERMS
50 m = m + 1.0D0
t = t * (s/(m*(nu+m)))
w = w + t
IF (anorm(t) > eps*anorm(w)) GO TO 50
RETURN
END SUBROUTINE cbjm
页:
[1]
2