PROGRAM J0ROOT C C Compute the roots of the Bessel function J0(x) and display C the results. C DOUBLE PRECISION X,Y C WRITE (5,10) 10 FORMAT(' ROOTS OF BESSEL FUNCTION J0'/) C DO 20 I=1,150 CALL BESJ0R (I,X) CALL BESJ0 (X,Y) WRITE (5,15) I,X,Y 15 FORMAT (' ROOT 'I3' = 'F23.18', J0(X) = 'F15.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 C***************************************************************************** SUBROUTINE BESJ0R (N,ROOT) C C Compute the Nth root of the Bessel function J0(x) in double C precision. C DOUBLE PRECISION SD,PID,BETAD,T1D,T2D,T3D,T4D,ROOT DOUBLE PRECISION ZEROS(150),START(20) C PARAMETER ( PID = 3.1415 92653 58979 324 D+0 ) PARAMETER ( T1D = 1.2500 00000 00000 000 D-1 ) PARAMETER ( T2D = 1.6145 83333 33333 333 D-1 ) PARAMETER ( T3D = 2.4602 86458 33333 333 D-1 ) PARAMETER ( T4D = 1.8244 38767 20610 118 D+0 ) C DATA START / 2.4048255572D0, 5.5200781103D0, 8.6537279129D0, 1 11.7915344391D0, 14.9309177086D0, 18.0710639379D0, 2 21.2116366299D0, 24.3524715308D0, 27.4934791320D0, 3 30.6346064684D0, 33.7758202136D0, 36.9170983537D0, 4 40.0584257646D0, 43.1997917132D0, 46.3411883717D0, 5 49.4826098974D0, 52.6240518411D0, 55.7655107550D0, 6 58.9069839261D0, 62.0484691902D0 / C IF (N .LE. 20) THEN ROOT = START(N) ELSE SD = N BETAD = (SD-.25D0)*PID ROOT = BETAD + T1D/BETAD + T2D/BETAD**3 + T3D/BETAD**5 1 + T4D/BETAD**7 ENDIF C RETURN END