www.pudn.com > Fortran.zip > MBSL3.FOR, change:1995-04-01,size:1589b

 
	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 
	END IF 
	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 
	MBSL3=P 
	RETURN 
	END