第二类修正贝塞尔函数(Fortran代码).doc

调试日期：2011年9月13日星期二程序说明：计算第二类修正贝塞尔函数的Fortran代码，参看徐士良先生的Fortran常用程序算法集PROGRAM BSL_XSLDOUBLE PRECISION MBSL4,XOPEN(1,FILE=BSL.DAT,ACTION=WRITE)DO X=0.05,3,0.05 WRITE(1,*),X,MBSL4(0,X-0.01),MBSL4(1,X)ENDDOCLOSE(1)ENDPROGRAMFUNCTION MBSL3(N,X)DOUBLE PRECISION MBSL3,XDOUBLE 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=-NT=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,-110 P=P*Y+A(I) ELSE Y=3.75/T P=C(9) DO 20 I=8,1,-120 P=P*Y+C(I) P=P*EXP(T)/SQRT(T) END IFEND IFIF (N.EQ.0) THEN MBSL3=P RETURNEND IFQ=PIF (T.LT.3.75) THEN Y=(X/3.75)*(X/3.75) P=B(7) DO 30 I=6,1,-130 P=P*Y+B(I) P=P*TELSE Y=3.75/T P=D(9) DO 40 I=8,1,-140 P=P*Y+D(I) P=P*EXP(T)/SQRT(T)END IFIF (X.LT.0.0) P=-PIF (N.EQ.1) THEN MBSL3=P RETURNEND IFIF (X+1.0.EQ.1.0) THEN MBSL3=0.0 RETURNEND IFY=2.0/TT=0.0B1=1.0B0=0.0M=SQRT(40.0*N)M=N+MM=M+MDO 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=B050CONTINUEP=T*Q/B1IF (X.LT.0.0).AND.(MOD(N,2).EQ.1) P=-PMBSL3=PRETURNENDFUNCTIONFUNCTION MBSL4(N,X)DOUBLE PRECISION MBSL4,X,MBSL3DOUBLE 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=-NIF (X.LT.0.0) X=-XIF (X+1.0.EQ.1.0) THEN MBSL4=1.0D+35 RETURNEND IFIF (N.NE.1) THEN IF (X.LE.2.0) THEN Y=X*X/4.0 P=A(7) DO 10 I=6,1,-110 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,-120 P=P*Y+C(I) P=P*EXP(-X)/SQRT(X) END IFEND IFIF (N.EQ.0) THEN MBSL4=P RETURNEND IFB0=PIF (X.LE.2.0) THEN Y=X*X/4.0 P=B(7) DO 30 I=6,1,-130 P=P*Y+B(I) P=P/X+MBSL3(1,X)*LOG(X/2.0)ELSE Y=2.0/X P=D(7) DO 40 I=6,1,-140 P=P*Y+D(I) P=P*EXP(-X)/SQRT(X)