!  "bessj0" function that returns the Bessel function J0(x) for any real x
 !
 !
 !   REFERENCE: "Numerical recipes in FORTRAN,
 !               The Art of Scientific Computing, Second Edition"
 !               Press, Teukolsky, Vetterling, Flannery
 !               Cambridge University Press
 !               ISBN 0-521-43064-X
 !       PAGES:  225 ("FUNCTION bessj0")
 !
 !       INPUT:  "x" the input argument  (REAL)
 !      OUTPUT:  "bessj0" is the value of J0(x) (REAL)

 REAL bessj0 x ;
 REAL ax xx z ;
 DOUBLE p1 p2 p3 p4 p5 := 1.D0 -.1098628627D-2 .2734510407D-4
                         -.2073370639D-5 .2093887211D-6 ;
 DOUBLE q1 q2 q3 q4 q5 := -.1562499995D-1 .1430488765D-3 -.6911147651D-5
                           .7621095161D-6 -.934945152D-7 ;
 DOUBLE r1 r2 r3 r4 r5 r6 := 57568490574.D0 -13362590354.D0
       651619640.7D0 -11214424.18D0 77392.33017D0 -184.9052456D0 ;
 DOUBLE s1 s2 s3 s4 s5 s6 := 57568490411.D0 1029532985.D0
       9494680.718D0 59272.64853D0 267.8532712D0 1.D0 ;
 DOUBLE y ;

 :: >>x<< ;
 IF x ABS 8. < THEN
    EVALUATE y := x x * R_TO_D ;
    EVALUATE bessj0 := r6 y * r5 + y * r4 + y * r3 + y * r2 + y * r1 +
                       s6 y * s5 + y * s4 + y * s3 + y * s2 + y * s1 + /
                       D_TO_R ;
 ELSE
    EVALUATE ax := x ABS ;
    EVALUATE z  xx := 8. ax / ax .785398164 - ;
    EVALUATE y := z z * R_TO_D ;
    EVALUATE bessj0 := p5 y * p4 + y * p3 + y * p2 + y * p1 +
                       xx COS R_TO_D *
                       q5 y * q4 + y * q3 + y * q2 + y * q1 +
                       xx SIN z * R_TO_D * -
                       .636619772 ax / SQRT R_TO_D *
                       D_TO_R ;
 ENDIF ;
 :: <<bessj0>> ;
 QUIT " Function *bessj0* XREF " .