FUNCTION BESSY1(X)
      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,S7
      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 R1,R2,R3,R4,R5,R6/-.4900604943D13,.1275274390D13,-.5153438139
     *D11,
     *    .7349264551D9,-.4237922726D7,.8511937935D4/,
     *    S1,S2,S3,S4,S5,S6,S7/.2499580570D14,.4244419664D12,
     *    .3733650367D10,.2245904002D8,.1020426050D6,.3549632885D3,1.D0/
      IF(X.LT.8.)THEN
        Y=X**2
        BESSY1=X*(R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6)))))/(S1+Y*(S2+Y*
     *      (S3+Y*(S4+Y*(S5+Y*(S6+Y*S7))))))+.636619772
     *      *(BESSJ1(X)*LOG(X)-1./X)
      ELSE
        Z=8./X
        Y=Z**2
        XX=X-2.356194491
        BESSY1=SQRT(.636619772/X)*(SIN(XX)*(P1+Y*(P2+Y*(P3+Y*(P4+Y
     *      *P5))))+Z*COS(XX)*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*Q5)))))
      ENDIF
      RETURN
      END