第二类修正贝塞尔函数(Fortran代码)
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
调试日期：2011年9月13日星期二
程序说明：计算第二类修正贝塞尔函数的Fortran代码，参看徐士良先生的《Fortran常用程序算法集》
PROGRAM BSL_XSL
DOUBLE PRECISION MBSL4,X
OPEN(1,FILE='BSL.DAT',ACTION='WRITE')
DO X=0.05,3,0.05
WRITE(1,*),X,MBSL4(0,X-0.01),MBSL4(1,X)
ENDDO
CLOSE(1)
ENDPROGRAM
FUNCTION MBSL3(N,X)
DOUBLE PRECISION MBSL3,X
DOUBLE PRECISION T,Y,P,B0,B1,Q,A(7),B(7),C(9),D(9)
DATA A/1.0,3.5156229,3.0899424,1.2067492,
* 0.2659732,0.0360768,0.0045813/
DATA B/0.5,0.87890594,0.51498869,0.15084934,
* 0.02658773,0.00301532,0.00032411/
DATA C/0.39894228,0.01328592,0.00225319,
* -0.00157565,0.00916281,-0.02057706,
* 0.02635537,-0.01647663,0.00392377/
DATA D/0.39894228,-0.03988024,-0.00362018,
* 0.00163801,-0.01031555,0.02282967,
* -0.02895312,0.01787654,-0.00420059/
IF (N.LT.0) N=-N
T=ABS(X)
IF (N.NE.1) THEN
IF (T.LT.3.75) THEN
Y=(X/3.75)*(X/3.75)
P=A(7)
DO 10 I=6,1,-1
10 P=P*Y+A(I)
ELSE
Y=3.75/T
P=C(9)
DO 20 I=8,1,-1
20 P=P*Y+C(I)
P=P*EXP(T)/SQRT(T)
END IF
END IF
IF (N.EQ.0) THEN
MBSL3=P
RETURN
Q=P
IF (T.LT.3.75) THEN
Y=(X/3.75)*(X/3.75)
P=B(7)
DO 30 I=6,1,-1
30 P=P*Y+B(I)
P=P*T
ELSE
Y=3.75/T
P=D(9)
DO 40 I=8,1,-1
40 P=P*Y+D(I)
P=P*EXP(T)/SQRT(T)
END IF
IF (X.LT.0.0) P=-P
IF (N.EQ.1) THEN
MBSL3=P
RETURN
END IF
IF (X+1.0.EQ.1.0) THEN
MBSL3=0.0
RETURN
END IF
Y=2.0/T
T=0.0
B1=1.0
B0=0.0
M=SQRT(40.0*N)
M=N+M
M=M+M
DO 50 I=M,1,-1
P=B0+I*Y*B1
B0=B1
B1=P
IF (ABS(B1).GT.1.0D+10) THEN
T=T*1.0D-10
B0=B0*1.0D-10
B1=B1*1.0D-10
END IF
IF (I.EQ.N) T=B0
50 CONTINUE
P=T*Q/B1
IF ((X.LT.0.0).AND.(MOD(N,2).EQ.1)) P=-P
RETURN
ENDFUNCTION
FUNCTION MBSL4(N,X)
DOUBLE PRECISION MBSL4,X,MBSL3
DOUBLE PRECISION Y,P,B0,B1,A(7),B(7),C(7),D(7)
DATA A/-0.57721566,0.4227842,0.23069756,0.0348859,
* 0.00262698,0.0001075,0.0000074/
DATA B/1.0,0.15443144,-0.67278579,-0.18156897,
* -0.01919402,-0.00110404,-0.00004686/
DATA C/1.25331414,-0.07832358,0.02189568,-0.01062446, * 0.00587872,-0.0025154,0.00053208/
DATA D/1.25331414,0.23498619,-0.0365562,0.01504268,
* -0.00780353,0.00325614,-0.00068245/
IF (N.LT.0) N=-N
IF (X.LT.0.0) X=-X
IF (X+1.0.EQ.1.0) THEN
MBSL4=1.0D+35
RETURN
END IF
IF (N.NE.1) THEN
IF (X.LE.2.0) THEN
Y=X*X/4.0
P=A(7)
DO 10 I=6,1,-1
10 P=P*Y+A(I)
P=P-MBSL3(0,X)*LOG(X/2.0)
ELSE
Y=2.0/X
P=C(7)
DO 20 I=6,1,-1
20 P=P*Y+C(I)
P=P*EXP(-X)/SQRT(X)
END IF
END IF
IF (N.EQ.0) THEN
MBSL4=P
RETURN
END IF
B0=P
IF (X.LE.2.0) THEN