	DOUBLE PRECISION FUNCTION BESSJ0(X)
	DOUBLE PRECISION X
	DOUBLE PRECISION Y, Z, AX, XX
	DOUBLE PRECISION P1,P2,P3,P4,P5
	DOUBLE PRECISION Q1,Q2,Q3,Q4,Q5
	DOUBLE PRECISION R1,R2,R3,R4,R5,R6
	DOUBLE PRECISION S1,S2,S3,S4,S5,S6

	DATA P1 /1.D0/
	DATA P2 /-.1098628627D-2/
	DATA P3 /.2734510407D-4/
	DATA P4 /-.2073370639D-5/
	DATA P5 /.2093887211D-6/

	DATA Q1 /-.1562499995D-1/
	DATA Q2 /.1430488765D-3/
	DATA Q3 /-.6911147651D-5/
	DATA Q4 /.7621095161D-6/
	DATA Q5 /-.934945152D-7/

	DATA R1 /57568490574.D0/
	DATA R2 /-13362590354.D0/
	DATA R3 /651619640.7D0/
	DATA R4 /-11214424.18D0/
	DATA R5 /77392.33017D0/
	DATA R6 /-184.9052456D0/

	DATA S1 /57568490411.D0/
	DATA S2 /1029532985.D0/
	DATA S3 /9494680.718D0/
	DATA S4 /59272.64853D0/
	DATA S5 /267.8532712D0/
	DATA S6 /1.D0/

	IF(ABS(X).LT.8.)THEN
		Y = X**2
		BESSJ0 = (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-.785398164
		BESSJ0 = 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)))))
	ENDIF
	RETURN
	END
