	function bessj1(x)
c
c       from numerical recipes by w.h. press, b.p. flannery,
c       s.a. teukolsky, and w.t. vetterling, 1986, cambridge press
c
	real*8 y,p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,
     *		s1,s2,s3,s4,s5,s6
	data r1 /72362614232.d0/
        data r2 /-7895059235.d0/
        data r3 /242396853.1d0/
        data r4 /-2972611.439d0/
        data r5 /15704.48260d0/
        data r6 /-30.16036606d0/    
        data s1 /144725228442.d0/
        data s2 /2300535178.d0/
        data s3 /18583304.74d0/
        data s4 /99447.43394d0/
        data s5 /376.9991397d0/
        data s6	/1.d0/
        data p1 /1.d0/
        data p2 /.183105d-2/
        data p3 /-.3516396496d-4/
        data p4 /.2457520174d-5/
        data p5 /-.240337019d-6/
        data q1 /.04687499995d0/
        data q2 /-.2002690873d-3/
        data q3 /.8449199096d-5/
        data q4 /-.88228987d-6/
        data q5 /.105787412d-6/
	if(abs(x).lt.8.)then
		y=x**2
		bessj1=x*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))
     *			/(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6)))))
	else
		ax=abs(x)
		z=8./ax
		y=z**2
		xx=ax-2.356194491
		bessj1=sqrt(.636619772/ax)*(cos(xx)*(p1+y*(p2+y*
     *             (p3+y*(p4+y*p5))))-z*sin(xx)*(q1+y*(q2+y*
     *             (q3+y*(q4+y*q5)))))*sign(1.,x)
	endif
	return
	end
