Fortran有没有可以调用半整数阶的贝塞尔函数库？

ii08842

用IVF2013调用MKL库的贝塞尔函数，结果看帮助文档中阶数n只能为整型变量。有没有编译了半整数阶的库函数可以直接调用的？

li913

半整数阶贝塞尔是初等函数，利用贝塞尔性质，可以自己写任意阶的代码。

zdhuang

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代码来自这本书，仅供参考。

[Fortran] 纯文本查看 复制代码

查看代码

打印?

001	SUBROUTINE JYV(V,X,VM,BJ,DJ,BY,DY)

002

C

003	C       =======================================================

004	C       Purpose: Compute Bessel functions Jv(x) and Yv(x)

005	C                and their derivatives

006	C       Input :  x --- Argument of Jv(x) and Yv(x)

007	C                v --- Order of Jv(x) and Yv(x)

008	C                      ( v = n+v0, 0 ף v0 < 1, n = 0,1,2,... )

009	C       Output:  BJ(n) --- Jn+v0(x)

010	C                DJ(n) --- Jn+v0'(x)

011	C                BY(n) --- Yn+v0(x)

012	C                DY(n) --- Yn+v0'(x)

013	C                VM --- Highest order computed

014	C       Routines called:

015	C            (1) GAMMA for computing gamma function

016	C            (2) MSTA1 and MSTA2 for computing the starting

017	C                point for backward recurrence

018	C       =======================================================

019

C

020	IMPLICIT DOUBLE PRECISION (A-H,O-Z)

021	DIMENSION BJ(0:*),DJ(0:*),BY(0:*),DY(0:*)

022	EL=.5772156649015329D0

023	PI=3.141592653589793D0

024	RP2=.63661977236758D0

025

X2=X*X

026

N=INT(V)

027

V0=V-N

028	IF (X.LT.1.0D-100) THEN

029	DO 10 K=0,N

030	BJ(K)=0.0D0

031	DJ(K)=0.0D0

032	BY(K)=-1.0D+300

033	10            DY(K)=1.0D+300

034	IF (V0.EQ.0.0) THEN

035	BJ(0)=1.0D0

036	DJ(1)=0.5D0

037

ELSE

038	DJ(0)=1.0D+300

039

ENDIF

040

VM=V

041

RETURN

042

ENDIF

043	IF (X.LE.12.0) THEN

044	DO 25 L=0,1

045

VL=V0+L

046	BJVL=1.0D0

047

R=1.0D0

048	DO 15 K=1,40

049	R=-0.25D0*R*X2/(K*(K+VL))

050	BJVL=BJVL+R

051	IF (DABS(R).LT.DABS(BJVL)*1.0D-15) GO TO 20

052	15            CONTINUE

053	20            VG=1.0D0+VL

054	CALL GAMMA(VG,GA)

055	A=(0.5D0*X)**VL/GA

056	IF (L.EQ.0) BJV0=BJVL*A

057	IF (L.EQ.1) BJV1=BJVL*A

058	25         CONTINUE

059

ELSE

060

K0=11

061	IF (X.GE.35.0) K0=10

062	IF (X.GE.50.0) K0=8

063	DO 40 J=0,1

064	VV=4.0D0*(J+V0)*(J+V0)

065

PX=1.0D0

066

RP=1.0D0

067	DO 30 K=1,K0

068	RP=-0.78125D-2*RP*(VV-(4.0*K-3.0)**2.0)*(VV-

069	&              (4.0*K-1.0)**2.0)/(K*(2.0*K-1.0)*X2)

070	30               PX=PX+RP

071

QX=1.0D0

072

RQ=1.0D0

073	DO 35 K=1,K0

074	RQ=-0.78125D-2*RQ*(VV-(4.0*K-1.0)**2.0)*(VV-

075	&              (4.0*K+1.0)**2.0)/(K*(2.0*K+1.0)*X2)

076	35               QX=QX+RQ

077	QX=0.125D0*(VV-1.0)*QX/X

078	XK=X-(0.5D0*(J+V0)+0.25D0)*PI

079	A0=DSQRT(RP2/X)

080	CK=DCOS(XK)

081	SK=DSIN(XK)

082	IF (J.EQ.0) THEN

083	BJV0=A0*(PX*CK-QX*SK)

084	BYV0=A0*(PX*SK+QX*CK)

085	ELSE IF (J.EQ.1) THEN

086	BJV1=A0*(PX*CK-QX*SK)

087	BYV1=A0*(PX*SK+QX*CK)

088

ENDIF

089	40         CONTINUE

090

ENDIF

091	BJ(0)=BJV0

092	BJ(1)=BJV1

093	DJ(0)=V0/X*BJ(0)-BJ(1)

094	DJ(1)=-(1.0D0+V0)/X*BJ(1)+BJ(0)

095	IF (N.GE.2.AND.N.LE.INT(0.9*X)) THEN

096

F0=BJV0

097

F1=BJV1

098	DO 45 K=2,N

099	F=2.0D0*(K+V0-1.0D0)/X*F1-F0

100

BJ(K)=F

101

F0=F1

102	45            F1=F

103	ELSE IF (N.GE.2) THEN

104	M=MSTA1(X,200)

105	IF (M.LT.N) THEN

106

N=M

107

ELSE

108	M=MSTA2(X,N,15)

109

ENDIF

110

F2=0.0D0

111	F1=1.0D-100

112	DO 50 K=M,0,-1

113	F=2.0D0*(V0+K+1.0D0)/X*F1-F2

114	IF (K.LE.N) BJ(K)=F

115

F2=F1

116	50            F1=F

117	IF (DABS(BJV0).GT.DABS(BJV1)) THEN

118

CS=BJV0/F

119

ELSE

120	CS=BJV1/F2

121

ENDIF

122	DO 55 K=0,N

123	55            BJ(K)=CS*BJ(K)

124

ENDIF

125	DO 60 K=2,N

126	60         DJ(K)=-(K+V0)/X*BJ(K)+BJ(K-1)

127	IF (X.LE.12.0D0) THEN

128	IF (V0.NE.0.0) THEN

129	DO 75 L=0,1

130

VL=V0+L

131	BJVL=1.0D0

132

R=1.0D0

133	DO 65 K=1,40

134	R=-0.25D0*R*X2/(K*(K-VL))

135	BJVL=BJVL+R

136	IF (DABS(R).LT.DABS(BJVL)*1.0D-15) GO TO 70

137	65               CONTINUE

138	70               VG=1.0D0-VL

139	CALL GAMMA(VG,GB)

140	B=(2.0D0/X)**VL/GB

141	IF (L.EQ.0) BJU0=BJVL*B

142	IF (L.EQ.1) BJU1=BJVL*B

143	75            CONTINUE

144

PV0=PI*V0

145	PV1=PI*(1.0D0+V0)

146	BYV0=(BJV0*DCOS(PV0)-BJU0)/DSIN(PV0)

147	BYV1=(BJV1*DCOS(PV1)-BJU1)/DSIN(PV1)

148

ELSE

149	EC=DLOG(X/2.0D0)+EL

150

CS0=0.0D0

151

W0=0.0D0

152

R0=1.0D0

153	DO 80 K=1,30

154	W0=W0+1.0D0/K

155	R0=-0.25D0*R0/(K*K)*X2

156	80               CS0=CS0+R0*W0

157	BYV0=RP2*(EC*BJV0-CS0)

158

CS1=1.0D0

159

W1=0.0D0

160

R1=1.0D0

161	DO 85 K=1,30

162	W1=W1+1.0D0/K

163	R1=-0.25D0*R1/(K*(K+1))*X2

164	85               CS1=CS1+R1*(2.0D0*W1+1.0D0/(K+1.0D0))

165	BYV1=RP2*(EC*BJV1-1.0D0/X-0.25D0*X*CS1)

166

ENDIF

167

ENDIF

168	BY(0)=BYV0

169	BY(1)=BYV1

170	DO 90 K=2,N

171	BYVK=2.0D0*(V0+K-1.0D0)/X*BYV1-BYV0

172	BY(K)=BYVK

173

BYV0=BYV1

174	90         BYV1=BYVK

175	DY(0)=V0/X*BY(0)-BY(1)

176	DO 95 K=1,N

177	95         DY(K)=-(K+V0)/X*BY(K)+BY(K-1)

178

VM=N+V0

179

RETURN

180

END

181

182

183	SUBROUTINE GAMMA(X,GA)

184

C

185	C       ==================================================

186	C       Purpose: Compute gamma function ג(x)

187	C       Input :  x  --- Argument of ג(x)

188	C                       ( x is not equal to 0,-1,-2,תתת)

189	C       Output:  GA --- ג(x)

190	C       ==================================================

191

C

192	IMPLICIT DOUBLE PRECISION (A-H,O-Z)

193	DIMENSION G(26)

194	PI=3.141592653589793D0

195	IF (X.EQ.INT(X)) THEN

196	IF (X.GT.0.0D0) THEN

197

GA=1.0D0

198

M1=X-1

199	DO 10 K=2,M1

200	10               GA=GA*K

201

ELSE

202	GA=1.0D+300

203

ENDIF

204

ELSE

205	IF (DABS(X).GT.1.0D0) THEN

206

Z=DABS(X)

207

M=INT(Z)

208

R=1.0D0

209	DO 15 K=1,M

210	15               R=R*(Z-K)

211

Z=Z-M

212

ELSE

213

Z=X

214

ENDIF

215	DATA G/1.0D0,0.5772156649015329D0,

216	&          -0.6558780715202538D0, -0.420026350340952D-1,

217	&          0.1665386113822915D0,-.421977345555443D-1,

218	&          -.96219715278770D-2, .72189432466630D-2,

219	&          -.11651675918591D-2, -.2152416741149D-3,

220	&          .1280502823882D-3, -.201348547807D-4,

221	&          -.12504934821D-5, .11330272320D-5,

222	&          -.2056338417D-6, .61160950D-8,

223	&          .50020075D-8, -.11812746D-8,

224	&          .1043427D-9, .77823D-11,

225	&          -.36968D-11, .51D-12,

226	&          -.206D-13, -.54D-14, .14D-14, .1D-15/

227

GR=G(26)

228	DO 20 K=25,1,-1

229	20            GR=GR*Z+G(K)

230	GA=1.0D0/(GR*Z)

231	IF (DABS(X).GT.1.0D0) THEN

232

GA=GA*R

233	IF (X.LT.0.0D0) GA=-PI/(X*GA*DSIN(PI*X))

234

ENDIF

235

ENDIF

236

RETURN

237

END

238

239

240	INTEGER FUNCTION MSTA1(X,MP)

241

C

242	C       ===================================================

243	C       Purpose: Determine the starting point for backward

244	C                recurrence such that the magnitude of

245	C                Jn(x) at that point is about 10^(-MP)

246	C       Input :  x     --- Argument of Jn(x)

247	C                MP    --- Value of magnitude

248	C       Output:  MSTA1 --- Starting point

249	C       ===================================================

250

C

251	IMPLICIT DOUBLE PRECISION (A-H,O-Z)

252	A0=DABS(X)

253	N0=INT(1.1*A0)+1

254	F0=ENVJ(N0,A0)-MP

255

N1=N0+5

256	F1=ENVJ(N1,A0)-MP

257	DO 10 IT=1,20

258	NN=N1-(N1-N0)/(1.0D0-F0/F1)

259	F=ENVJ(NN,A0)-MP

260	IF(ABS(NN-N1).LT.1) GO TO 20

261

N0=N1

262

F0=F1

263

N1=NN

264	10        F1=F

265	20     MSTA1=NN

266

RETURN

267

END

268

269

270	INTEGER FUNCTION MSTA2(X,N,MP)

271

C

272	C       ===================================================

273	C       Purpose: Determine the starting point for backward

274	C                recurrence such that all Jn(x) has MP

275	C                significant digits

276	C       Input :  x  --- Argument of Jn(x)

277	C                n  --- Order of Jn(x)

278	C                MP --- Significant digit

279	C       Output:  MSTA2 --- Starting point

280	C       ===================================================

281

C

282	IMPLICIT DOUBLE PRECISION (A-H,O-Z)

283	A0=DABS(X)

284	HMP=0.5D0*MP

285	EJN=ENVJ(N,A0)

286	IF (EJN.LE.HMP) THEN

287

OBJ=MP

288	N0=INT(1.1*A0)

289

ELSE

290	OBJ=HMP+EJN

291

N0=N

292

ENDIF

293	F0=ENVJ(N0,A0)-OBJ

294

N1=N0+5

295	F1=ENVJ(N1,A0)-OBJ

296	DO 10 IT=1,20

297	NN=N1-(N1-N0)/(1.0D0-F0/F1)

298	F=ENVJ(NN,A0)-OBJ

299	IF (ABS(NN-N1).LT.1) GO TO 20

300

N0=N1

301

F0=F1

302

N1=NN

303	10         F1=F

304	20      MSTA2=NN+10

305

RETURN

306

END

307

308	REAL*8 FUNCTION ENVJ(N,X)

309	DOUBLE PRECISION X

310	ENVJ=0.5D0*DLOG10(6.28D0*N)-N*DLOG10(1.36D0*X/N)

311

RETURN

312

END