PROGRAM J1 C C Compute and print out values of the Bessel function J1(x) C DOUBLE PRECISION X,Y C WRITE (5,10) 10 FORMAT(' BESSEL FUNCTION J1'/) C DO 20 I=0,100 N = I X = FLOAT(I)/10. CALL BESJ1 (X,Y) WRITE (5,15) X,Y 15 FORMAT (' BJ1('F5.1') = 'D19.12) 20 CONTINUE C STOP END C***************************************************************************** SUBROUTINE BESJ1 (X,Y) C C Compute the Bessel function J1(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 370. C C Error: C X < 3: 1.3E-8 C X > 3: 9.0E-8 C DOUBLE PRECISION X,Y,XD3,XD3SQ,F1,THETA C IF (X .LT. 3.D0) THEN XD3SQ = X*X/9D0 Y = X*(.5D0+(-.56249985D0+(.21093573D0+(-.03954289D0 1 +(.00443319D0+(-.00031761D0+(.00001109D0 2 )*XD3SQ)*XD3SQ)*XD3SQ)*XD3SQ)*XD3SQ)*XD3SQ) ELSE XD3 = 3D0/X THETA = X-2.35619449D0+(.12499612D0+(.0000565D0 1 +(-.00637879D0+(.00074348D0+(.00079824D0+(-.00029166D0 2 )*XD3)*XD3)*XD3)*XD3)*XD3)*XD3 F1 = .79788456D0+(.00000156D0+(.01659667D0+(.00017105D0 1 +(-.00249511D0+(.00113653D0+(-.00020033D0 2 )*XD3)*XD3)*XD3)*XD3)*XD3)*XD3 Y = DSQRT(1D0/X)*F1*DCOS(THETA) ENDIF C RETURN END