	PROGRAM J0
C
C	Compute and print out values of the Bessel function J0(x)
C
	DOUBLE PRECISION X,Y
C
	WRITE (5,10)
10	FORMAT(' BESSEL FUNCTION J0'/)
C
	DO 20 I=0,100
	  N = I
	  X = FLOAT(I)/10.
	  CALL BESJ0 (X,Y)
	  WRITE (5,15) X,Y
15	  FORMAT (' BJ0('F5.1') = 'D19.12)
20	CONTINUE
C
	STOP
	END
C*****************************************************************************
	SUBROUTINE BESJ0 (X,Y)
C
C	Compute the Bessel function J0(x) in double precision for positive
C	x.  The value is returned in Y.
C
C	The formulae for this subroutine were taken from the
C	   Handbook of Mathematical Functions
C	   U.S. Dept. of Commerce
C	   Page 369-370.
C
C	Error:
C	   X < 3:    5.0E-8
C	   X > 3:    1.6E-8
C
	DOUBLE PRECISION X,Y,XD3,XD3SQ,F1,THETA
C
	IF (X .LT. 3.D0) THEN
	  XD3SQ = X*X/9D0
	  Y = 1D0+(-2.2499997D0+(1.2656208D0+(-.3163866D0
     1	  +(.0444479D0+(-.0039444D0+(.00021D0
     2	  )*XD3SQ)*XD3SQ)*XD3SQ)*XD3SQ)*XD3SQ)*XD3SQ
	ELSE
	  XD3 = 3D0/X
	  THETA = X-.78539816D0+(-.04166397D0+(-.00003954D0
     1	  +(.00262573D0+(-.00054125D0+(-.00029333D0+(.00013558D0
     2	  )*XD3)*XD3)*XD3)*XD3)*XD3)*XD3
	  F1 = .79788456D0+(-.00000077D0+(-.0055274D0+(-.00009512D0
     1	  +(.00137237D0+(-.00072805D0+(.00014476D0
     2	  )*XD3)*XD3)*XD3)*XD3)*XD3)*XD3
	   Y = DSQRT(1D0/X)*F1*DCOS(THETA)
	ENDIF
C
	RETURN
	END
������������������������������������������������������������������������������������������������������������������������������������������������������������