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

 
	FUNCTION MBSL1(N,X) 
	DOUBLE PRECISION MBSL1,X 
	DOUBLE PRECISION T,Y,Z,P,Q,S,B0,B1 
	DOUBLE PRECISION A(6),B(6),C(6),D(6),E(5),F(5),G(5),H(5) 
	DATA A/57568490574.0,-13362590354.0,651619640.7, 
     *         -11214424.18,77392.33017,-184.9052456/ 
	DATA B/57568490411.0,1029532985.0,9494680.718, 
     *         59272.64853,267.8532712,1.0/ 
	DATA C/72362614232.0,-7895059235.0,242396853.1, 
     *         -2972611.439,15704.4826,-30.16036606/ 
	DATA D/144725228443.0,2300535178.0,18583304.74, 
     *         99447.43394,376.9991397,1.0/ 
	DATA E/1.0,-0.1098628627D-02,0.2734510407D-04, 
     *         -0.2073370639D-05,0.2093887211D-06/ 
	DATA F/-0.1562499995D-01,0.1430488765D-03,-0.6911147651D-05, 
     *         0.7621095161D-06,-0.934935152D-07/ 
	DATA G/1.0,0.183105D-02,-0.3516396496D-04, 
     *         0.2457520174D-05,-0.240337019D-06/ 
	DATA H/0.4687499995D-01,-0.2002690873D-03,0.8449199096D-05, 
     *         -0.88228987D-06,0.105787412D-06/ 
	T=ABS(X) 
	IF (N.LT.0) N=-N 
	IF (N.NE.1) THEN 
	  IF (T.LT.8.0) THEN 
	    Y=T*T 
	    P=A(6) 
	    Q=B(6) 
	    DO 10 I=5,1,-1 
	      P=P*Y+A(I) 
	      Q=Q*Y+B(I) 
10	    CONTINUE 
	    P=P/Q 
	  ELSE 
	    Z=8.0/T 
	    Y=Z*Z 
	    P=E(5) 
	    Q=F(5) 
	    DO 20 I=4,1,-1 
	      P=P*Y+E(I) 
	      Q=Q*Y+F(I) 
20	    CONTINUE 
	    S=T-0.785398164 
	    P=P*COS(S)-Z*Q*SIN(S) 
	    P=P*SQRT(0.636619772/T) 
	  END IF 
	END IF 
	IF (N.EQ.0) THEN 
	  MBSL1=P 
	  RETURN 
	END IF 
	B0=P 
	IF (T.LT.8.0) THEN 
	  Y=T*T 
	  P=C(6) 
	  Q=D(6) 
	  DO 30 I=5,1,-1 
	    P=P*Y+C(I) 
	    Q=Q*Y+D(I) 
30	  CONTINUE 
	  P=X*P/Q 
	ELSE 
	  Z=8.0/T 
	  Y=Z*Z 
	  P=G(5) 
	  Q=H(5) 
	  DO 40 I=4,1,-1 
	    P=P*Y+G(I) 
	    Q=Q*Y+H(I) 
40	  CONTINUE 
	  S=T-2.356194491 
	  P=P*COS(S)-Z*Q*SIN(S) 
	  P=P*X*SQRT(0.636619772/T)/T 
	END IF 
	IF (N.EQ.1) THEN 
	  MBSL1=P 
	  RETURN 
	END IF 
	B1=P 
	IF (X.EQ.0.0) THEN 
	  MBSL1=0.0 
	  RETURN 
	END IF 
	S=2.0/T 
	IF (T.GT.1.0*N) THEN 
	  IF (X.LT.0.0) B1=-B1 
	  DO 50 I=1,N-1 
	    P=S*I*B1-B0 
	    B0=B1 
	    B1=P 
50	  CONTINUE 
	ELSE 
	  M=SQRT(40.0*N) 
	  M=(N+M)/2 
	  M=M+M 
	  P=0.0 
	  Q=0.0 
	  B0=1.0 
	  B1=0.0 
	  DO 60 I=M-1,0,-1 
	    T=S*(I+1)*B0-B1 
	    B1=B0 
	    B0=T 
	    IF (ABS(B0).GT.1.0D+10) THEN 
	      B0=B0*1.0D-10 
	      B1=B1*1.0D-10 
	      P=P*1.0D-10 
	      Q=Q*1.0D-10 
	    END IF 
	    IF (MOD(I+2,2).EQ.0) Q=Q+B0 
	    IF (i+1.EQ.N) P=B1 
60	  CONTINUE 
	  Q=2.0*Q-B0 
	  P=P/Q 
	END IF 
	IF ((X.LT.0.0).AND.(MOD(N,2).EQ.1)) P=-P 
	MBSL1=P 
	RETURN 
	END