Fortran有没有可以调用半整数阶的贝塞尔函数库？
用IVF2013调用MKL库的贝塞尔函数，结果看帮助文档中阶数n只能为整型变量。有没有编译了半整数阶的库函数可以直接调用的？半整数阶贝塞尔是初等函数，利用贝塞尔性质，可以自己写任意阶的代码。 JYV -- Evaluate a sequence of Bessel functions of the first and
second kinds and their derivatives with arbitrary real orders
and real arguments.
书名《COMPUTATION OF SPECIAL FUNCTIONS》
作者 Shanjie Zhang and Jianming Jin
以下Fortran77代码来自这本书，仅供参考。
SUBROUTINE JYV(V,X,VM,BJ,DJ,BY,DY)
C
C =======================================================
C Purpose: Compute Bessel functions Jv(x) and Yv(x)
C and their derivatives
C Input :x --- Argument of Jv(x) and Yv(x)
C v --- Order of Jv(x) and Yv(x)
C ( v = n+v0, 0 ף v0 < 1, n = 0,1,2,... )
C Output:BJ(n) --- Jn+v0(x)
C DJ(n) --- Jn+v0'(x)
C BY(n) --- Yn+v0(x)
C DY(n) --- Yn+v0'(x)
C VM --- Highest order computed
C Routines called:
C (1) GAMMA for computing gamma function
C (2) MSTA1 and MSTA2 for computing the starting
C point for backward recurrence
C =======================================================
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION BJ(0:*),DJ(0:*),BY(0:*),DY(0:*)
EL=.5772156649015329D0
PI=3.141592653589793D0
RP2=.63661977236758D0
X2=X*X
N=INT(V)
V0=V-N
IF (X.LT.1.0D-100) THEN
DO 10 K=0,N
BJ(K)=0.0D0
DJ(K)=0.0D0
BY(K)=-1.0D+300
10 DY(K)=1.0D+300
IF (V0.EQ.0.0) THEN
BJ(0)=1.0D0
DJ(1)=0.5D0
ELSE
DJ(0)=1.0D+300
ENDIF
VM=V
RETURN
ENDIF
IF (X.LE.12.0) THEN
DO 25 L=0,1
VL=V0+L
BJVL=1.0D0
R=1.0D0
DO 15 K=1,40
R=-0.25D0*R*X2/(K*(K+VL))
BJVL=BJVL+R
IF (DABS(R).LT.DABS(BJVL)*1.0D-15) GO TO 20
15 CONTINUE
20 VG=1.0D0+VL
CALL GAMMA(VG,GA)
A=(0.5D0*X)**VL/GA
IF (L.EQ.0) BJV0=BJVL*A
IF (L.EQ.1) BJV1=BJVL*A
25 CONTINUE
ELSE
K0=11
IF (X.GE.35.0) K0=10
IF (X.GE.50.0) K0=8
DO 40 J=0,1
VV=4.0D0*(J+V0)*(J+V0)
PX=1.0D0
RP=1.0D0
DO 30 K=1,K0
RP=-0.78125D-2*RP*(VV-(4.0*K-3.0)**2.0)*(VV-
& (4.0*K-1.0)**2.0)/(K*(2.0*K-1.0)*X2)
30 PX=PX+RP
QX=1.0D0
RQ=1.0D0
DO 35 K=1,K0
RQ=-0.78125D-2*RQ*(VV-(4.0*K-1.0)**2.0)*(VV-
& (4.0*K+1.0)**2.0)/(K*(2.0*K+1.0)*X2)
35 QX=QX+RQ
QX=0.125D0*(VV-1.0)*QX/X
XK=X-(0.5D0*(J+V0)+0.25D0)*PI
A0=DSQRT(RP2/X)
CK=DCOS(XK)
SK=DSIN(XK)
IF (J.EQ.0) THEN
BJV0=A0*(PX*CK-QX*SK)
BYV0=A0*(PX*SK+QX*CK)
ELSE IF (J.EQ.1) THEN
BJV1=A0*(PX*CK-QX*SK)
BYV1=A0*(PX*SK+QX*CK)
ENDIF
40 CONTINUE
ENDIF
BJ(0)=BJV0
BJ(1)=BJV1
DJ(0)=V0/X*BJ(0)-BJ(1)
DJ(1)=-(1.0D0+V0)/X*BJ(1)+BJ(0)
IF (N.GE.2.AND.N.LE.INT(0.9*X)) THEN
F0=BJV0
F1=BJV1
DO 45 K=2,N
F=2.0D0*(K+V0-1.0D0)/X*F1-F0
BJ(K)=F
F0=F1
45 F1=F
ELSE IF (N.GE.2) THEN
M=MSTA1(X,200)
IF (M.LT.N) THEN
N=M
ELSE
M=MSTA2(X,N,15)
ENDIF
F2=0.0D0
F1=1.0D-100
DO 50 K=M,0,-1
F=2.0D0*(V0+K+1.0D0)/X*F1-F2
IF (K.LE.N) BJ(K)=F
F2=F1
50 F1=F
IF (DABS(BJV0).GT.DABS(BJV1)) THEN
CS=BJV0/F
ELSE
CS=BJV1/F2
ENDIF
DO 55 K=0,N
55 BJ(K)=CS*BJ(K)
ENDIF
DO 60 K=2,N
60 DJ(K)=-(K+V0)/X*BJ(K)+BJ(K-1)
IF (X.LE.12.0D0) THEN
IF (V0.NE.0.0) THEN
DO 75 L=0,1
VL=V0+L
BJVL=1.0D0
R=1.0D0
DO 65 K=1,40
R=-0.25D0*R*X2/(K*(K-VL))
BJVL=BJVL+R
IF (DABS(R).LT.DABS(BJVL)*1.0D-15) GO TO 70
65 CONTINUE
70 VG=1.0D0-VL
CALL GAMMA(VG,GB)
B=(2.0D0/X)**VL/GB
IF (L.EQ.0) BJU0=BJVL*B
IF (L.EQ.1) BJU1=BJVL*B
75 CONTINUE
PV0=PI*V0
PV1=PI*(1.0D0+V0)
BYV0=(BJV0*DCOS(PV0)-BJU0)/DSIN(PV0)
BYV1=(BJV1*DCOS(PV1)-BJU1)/DSIN(PV1)
ELSE
EC=DLOG(X/2.0D0)+EL
CS0=0.0D0
W0=0.0D0
R0=1.0D0
DO 80 K=1,30
W0=W0+1.0D0/K
R0=-0.25D0*R0/(K*K)*X2
80 CS0=CS0+R0*W0
BYV0=RP2*(EC*BJV0-CS0)
CS1=1.0D0
W1=0.0D0
R1=1.0D0
DO 85 K=1,30
W1=W1+1.0D0/K
R1=-0.25D0*R1/(K*(K+1))*X2
85 CS1=CS1+R1*(2.0D0*W1+1.0D0/(K+1.0D0))
BYV1=RP2*(EC*BJV1-1.0D0/X-0.25D0*X*CS1)
ENDIF
ENDIF
BY(0)=BYV0
BY(1)=BYV1
DO 90 K=2,N
BYVK=2.0D0*(V0+K-1.0D0)/X*BYV1-BYV0
BY(K)=BYVK
BYV0=BYV1
90 BYV1=BYVK
DY(0)=V0/X*BY(0)-BY(1)
DO 95 K=1,N
95 DY(K)=-(K+V0)/X*BY(K)+BY(K-1)
VM=N+V0
RETURN
END
SUBROUTINE GAMMA(X,GA)
C
C ==================================================
C Purpose: Compute gamma function ג(x)
C Input :x--- Argument of ג(x)
C ( x is not equal to 0,-1,-2,תתת)
C Output:GA --- ג(x)
C ==================================================
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION G(26)
PI=3.141592653589793D0
IF (X.EQ.INT(X)) THEN
IF (X.GT.0.0D0) THEN
GA=1.0D0
M1=X-1
DO 10 K=2,M1
10 GA=GA*K
ELSE
GA=1.0D+300
ENDIF
ELSE
IF (DABS(X).GT.1.0D0) THEN
Z=DABS(X)
M=INT(Z)
R=1.0D0
DO 15 K=1,M
15 R=R*(Z-K)
Z=Z-M
ELSE
Z=X
ENDIF
DATA G/1.0D0,0.5772156649015329D0,
& -0.6558780715202538D0, -0.420026350340952D-1,
& 0.1665386113822915D0,-.421977345555443D-1,
& -.96219715278770D-2, .72189432466630D-2,
& -.11651675918591D-2, -.2152416741149D-3,
& .1280502823882D-3, -.201348547807D-4,
& -.12504934821D-5, .11330272320D-5,
& -.2056338417D-6, .61160950D-8,
& .50020075D-8, -.11812746D-8,
& .1043427D-9, .77823D-11,
& -.36968D-11, .51D-12,
& -.206D-13, -.54D-14, .14D-14, .1D-15/
GR=G(26)
DO 20 K=25,1,-1
20 GR=GR*Z+G(K)
GA=1.0D0/(GR*Z)
IF (DABS(X).GT.1.0D0) THEN
GA=GA*R
IF (X.LT.0.0D0) GA=-PI/(X*GA*DSIN(PI*X))
ENDIF
ENDIF
RETURN
END
INTEGER FUNCTION MSTA1(X,MP)
C
C ===================================================
C Purpose: Determine the starting point for backward
C recurrence such that the magnitude of
C Jn(x) at that point is about 10^(-MP)
C Input :x --- Argument of Jn(x)
C MP --- Value of magnitude
C Output:MSTA1 --- Starting point
C ===================================================
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
A0=DABS(X)
N0=INT(1.1*A0)+1
F0=ENVJ(N0,A0)-MP
N1=N0+5
F1=ENVJ(N1,A0)-MP
DO 10 IT=1,20
NN=N1-(N1-N0)/(1.0D0-F0/F1)
F=ENVJ(NN,A0)-MP
IF(ABS(NN-N1).LT.1) GO TO 20
N0=N1
F0=F1
N1=NN
10 F1=F
20 MSTA1=NN
RETURN
END
INTEGER FUNCTION MSTA2(X,N,MP)
C
C ===================================================
C Purpose: Determine the starting point for backward
C recurrence such that all Jn(x) has MP
C significant digits
C Input :x--- Argument of Jn(x)
C n--- Order of Jn(x)
C MP --- Significant digit
C Output:MSTA2 --- Starting point
C ===================================================
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
A0=DABS(X)
HMP=0.5D0*MP
EJN=ENVJ(N,A0)
IF (EJN.LE.HMP) THEN
OBJ=MP
N0=INT(1.1*A0)
ELSE
OBJ=HMP+EJN
N0=N
ENDIF
F0=ENVJ(N0,A0)-OBJ
N1=N0+5
F1=ENVJ(N1,A0)-OBJ
DO 10 IT=1,20
NN=N1-(N1-N0)/(1.0D0-F0/F1)
F=ENVJ(NN,A0)-OBJ
IF (ABS(NN-N1).LT.1) GO TO 20
N0=N1
F0=F1
N1=NN
10 F1=F
20 MSTA2=NN+10
RETURN
END
REAL*8 FUNCTION ENVJ(N,X)
DOUBLE PRECISION X
ENVJ=0.5D0*DLOG10(6.28D0*N)-N*DLOG10(1.36D0*X/N)
RETURN
END
页:
[1]