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