REAL*8 FUNCTION BESSJ0(X)
C From Numerical Recipes, p. 172
C Modified 26 Feb 1991 (J.T. Armstrong):  added ten terms of series 
C representation for J0 for small arguments.  The series representation
C has double-precision accuracy (15 digits) for arguments less than 
C 1.8, and single-precision accuracy (8 digits) for arguments less
C than 3.4.  The Numerical Recipes algorithms have single-precision
C accuracy.
      REAL*8 X,Y,P1,P2,P3,P4,P5,Q1,Q2,Q3,Q4,Q5,R1,R2,R3,R4,R5,R6,
     *    S1,S2,S3,S4,S5,S6,T1,T2,T3,T4,T5,T6,T7,T8,T9,T10
      DATA P1,P2,P3,P4,P5/1.D0,-.1098628627D-2,.2734510407D-4,
     *    -.2073370639D-5,.2093887211D-6/, 
     *    Q1,Q2,Q3,Q4,Q5/-.1562499995D-1,.1430488765D-3,-.6911147651D-5,
     *    .7621095161D-6,-.934945152D-7/
      DATA R1,R2,R3,R4,R5,R6/57568490574.D0,-13362590354.D0,
     *    651619640.7D0,-11214424.18D0,77392.33017D0,-184.9052456D0/,
     *    S1,S2,S3,S4,S5,S6/57568490411.D0,1029532985.D0,
     *    9494680.718D0,59272.64853D0,267.8532712D0,1.D0/
      DATA T1,T2,T3,T4,T5,T6,T7,T8,T9,T10
     *    /1D0,1D0,.25D0,.111111111111111D0,.0625D0,
     *    .04D0,.277777777777778D-1,.204081632653061D-1,
     *    .015625D0,.123456790123457D-1/
      IF(ABS(X).LT.3.4) THEN
        Y=(X**2)/4D0
        BESSJ0=T1-Y*T2*(T1-Y*T3*(T1-Y*T4*(T1-Y*T5*(T1-Y*T6*(T1-Y*T7*
     1         (T1-Y*T8*(T1-Y*T9*(T1-Y*T10))))))))
      ELSE IF(ABS(X).LT.8.)THEN
        Y=X**2
        BESSJ0=(R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6)))))
     1        /(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*
     1      (P4+Y*P5))))-Z*SIN(XX)*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*Q5)))))
      ENDIF
      RETURN
      END
C=======================================================================
      REAL*8 FUNCTION BESSJ1(X)
C From Numerical Recipes, p. 173
C Modified 26 Feb 1991 (J.T. Armstrong):  added thirteen terms of series 
C representation for J1 for small arguments.  The series representation
C has double-precision accuracy (15 digits) for arguments less than 
C about 2.5, and single-precision accuracy (8 digits) for arguments less
C than 5.5.  The Numerical Recipes algorithms have single-precision
C accuracy.
      REAL*8 X,Y,P1,P2,P3,P4,P5,Q1,Q2,Q3,Q4,Q5,R1,R2,R3,R4,R5,R6,
     *    S1,S2,S3,S4,S5,S6,T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11,T12,T13
      DATA R1,R2,R3,R4,R5,R6/72362614232.D0,-7895059235.D0,
     *    242396853.1D0,-2972611.439D0,15704.48260D0,-30.16036606D0/,
     *    S1,S2,S3,S4,S5,S6/144725228442.D0,2300535178.D0,
     *    18583304.74D0,99447.43394D0,376.9991397D0,1.D0/
      DATA P1,P2,P3,P4,P5/1.D0,.183105D-2,-.3516396496D-4,
     *    .2457520174D-5,-.240337019D-6/,
     *     Q1,Q2,Q3,Q4,Q5/.04687499995D0,-.2002690873D-3,
     *    .8449199096D-5,-.88228987D-6,.105787412D-6/
      DATA T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11,T12,T13
     *    /1D0,.5D0,.166666666666667D0,.833333333333333D-1,
     *    .05D0,.333333333333333D-1,.238095238095238D-1,
     *    .178571428571429D-1,.138888888888889D-1,.111111111111111D-1,
     *    .909090909090909D-2,.757575757575758D-2,.641025641025641D-2/
      IF(ABS(X).LT.5.5) THEN
        Y=(X**2)/4D0
        BESSJ1=(X/2D0)*(T1-T2*Y*(T1-T3*Y*(T1-T4*Y*(T1-T5*Y*(T1-T6*Y*
     1         (T1-T7*Y*(T1-T8*Y*(T1-T9*Y*(T1-T10*Y*(T1-T11*Y*
     2         (T1-T12*Y*(T1-T13*Y))))))))))))
      ELSE IF(ABS(X).LT.8.)THEN
        Y=X**2
        BESSJ1=X*(R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6)))))
     1          /(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*
     1      (P4+Y*P5))))-Z*SIN(XX)*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*Q5)))))
     2      *DSIGN(1.0D0,X)
      ENDIF
      RETURN
      END