*
* $Id: besy1.F,v 1.1.1.1 1996/02/15 17:49:08 mclareni Exp $
*
* $Log: besy1.F,v $
* Revision 1.1.1.1  1996/02/15 17:49:08  mclareni
* Kernlib
*
*
#include "kernnum/pilot.h"
      REAL FUNCTION BESY1(RX)
      REAL RX,SX
      CHARACTER*6 ENAME
      LOGICAL MFLAG,RFLAG
#if defined(CERNLIB_NUMHIPRE)
      REAL DBESY1,X,Y,H,ALFA,ZERO,ONE,TWO,EIGHT,B0,B1,B2,P,Q,R,D
      REAL PI1,PI2,PI3,CE,C1(0:14),C2(0:14),C3(0:7),C4(0:10)
#endif
#if defined(CERNLIB_NUMLOPRE)
      DOUBLE PRECISION X,Y,H,ALFA,ZERO,ONE,TWO,EIGHT,B0,B1,B2,P,Q,R,D
      DOUBLE PRECISION PI1,PI2,PI3,CE,C1(0:14),C2(0:14),C3(0:7),C4(0:10)
      DOUBLE PRECISION DBESY1,DX
#endif
 
      DATA ZERO /0.0D0/, ONE /1.0D0/, TWO /2.0D0/, EIGHT /8.0D0/
 
      DATA PI1 /0.79788 45608 0287D0/, PI2 /0.63661 97723 6758D0/
      DATA PI3 /2.35619 44901 923 D0/, CE  /0.57721 56649 0153D0/
 
      DATA C1( 0) /+0.05245 81903 3466D0/
      DATA C1( 1) /+0.04809 64691 5823D0/
      DATA C1( 2) /+0.31327 50823 6157D0/
      DATA C1( 3) /-0.24186 74084 4741D0/
      DATA C1( 4) /+0.07426 67962 1679D0/
      DATA C1( 5) /-0.01296 76273 1174D0/
      DATA C1( 6) /+0.00148 99128 9667D0/
      DATA C1( 7) /-0.00012 22786 8504D0/
      DATA C1( 8) /+0.00000 75626 3023D0/
      DATA C1( 9) /-0.00000 03661 3086D0/
      DATA C1(10) /+0.00000 00142 7732D0/
      DATA C1(11) /-0.00000 00004 5857D0/
      DATA C1(12) /+0.00000 00000 1235D0/
      DATA C1(13) /-0.00000 00000 0028D0/
      DATA C1(14) /+0.00000 00000 0001D0/
 
      DATA C2( 0) /-0.04017 29465 4441D0/
      DATA C2( 1) /-0.44444 71476 3056D0/
      DATA C2( 2) /-0.02271 92444 2842D0/
      DATA C2( 3) /+0.20664 45410 1749D0/
      DATA C2( 4) /-0.08667 16970 5695D0/
      DATA C2( 5) /+0.01763 67030 0316D0/
      DATA C2( 6) /-0.00223 56192 9449D0/
      DATA C2( 7) /+0.00019 70623 0270D0/
      DATA C2( 8) /-0.00001 28858 5330D0/
      DATA C2( 9) /+0.00000 06528 4795D0/
      DATA C2(10) /-0.00000 00264 5074D0/
      DATA C2(11) /+0.00000 00008 7803D0/
      DATA C2(12) /-0.00000 00000 2434D0/
      DATA C2(13) /+0.00000 00000 0057D0/
      DATA C2(14) /-0.00000 00000 0001D0/
 
      DATA C3( 0) /+1.00090 30408 600D0/
      DATA C3( 1) /+0.00089 89898 331D0/
      DATA C3( 2) /-0.00000 39872 843D0/
      DATA C3( 3) /+0.00000 00617 763D0/
      DATA C3( 4) /-0.00000 00018 719D0/
      DATA C3( 5) /+0.00000 00000 882D0/
      DATA C3( 6) /-0.00000 00000 057D0/
      DATA C3( 7) /+0.00000 00000 005D0/
 
      DATA C4( 0) /+0.04677 77870 69525D0/
      DATA C4( 1) /-0.00009 62772 35492D0/
      DATA C4( 2) /+0.00000 09138 61526D0/
      DATA C4( 3) /-0.00000 00209 59781D0/
      DATA C4( 4) /+0.00000 00008 22919D0/
      DATA C4( 5) /-0.00000 00000 46864D0/
      DATA C4( 6) /+0.00000 00000 03515D0/
      DATA C4( 7) /-0.00000 00000 00326D0/
      DATA C4( 8) /+0.00000 00000 00036D0/
      DATA C4( 9) /-0.00000 00000 00005D0/
      DATA C4(10) /+0.00000 00000 00001D0/
#if defined(CERNLIB_NUMHIPRE)
      ROUND(D)  =  D
#endif
#if defined(CERNLIB_NUMLOPRE)
      ROUND(D)  =  SNGL(D+(D-DBLE(SNGL(D))))
#endif
 
      ENAME=' BESY1'
      X=RX
#if defined(CERNLIB_NUMLOPRE)
      GOTO 9
      ENTRY DBESY1(DX)
      ENAME='DBESY1'
      X=DX
#endif
 
    9 IF(X .LE. ZERO) THEN
       CALL KERMTR('C312.1',LGFILE,MFLAG,RFLAG)
       IF(MFLAG) THEN
        SX=X
        IF(LGFILE .EQ. 0) THEN
         WRITE(*,100) ENAME,SX
        ELSE
         WRITE(LGFILE,100) ENAME,SX
        ENDIF
       ENDIF
       IF(.NOT.RFLAG) CALL ABEND
       IF(ENAME .EQ. ' BESY1')  THEN
          BESY1=ZERO
       ELSE
          DBESY1=ZERO
       ENDIF
       RETURN
      ENDIF
      IF(X .LT. EIGHT) THEN
       Y=X/EIGHT
       H=TWO*Y**2-ONE
       ALFA=-TWO*H
       B1=ZERO
       B2=ZERO
       DO 1 I = 14,0,-1
       B0=C1(I)-ALFA*B1-B2
       B2=B1
    1  B1=B0
       P=Y*(B0-B2)
       B1=ZERO
       B2=ZERO
       DO 2 I = 14,0,-1
       B0=C2(I)-ALFA*B1-B2
       B2=B1
    2  B1=B0
       B1=PI2*(CE+LOG(X/TWO))*P-PI2/X+Y*(B0-B2)
      ELSE
       R=ONE/X
       Y=EIGHT*R
       H=TWO*Y**2-ONE
       ALFA=-TWO*H
       B1=ZERO
       B2=ZERO
       DO 3 I = 7,0,-1
       B0=C3(I)-ALFA*B1-B2
       B2=B1
    3  B1=B0
       P=B0-H*B2
       B1=ZERO
       B2=ZERO
       DO 4 I = 10,0,-1
       B0=C4(I)-ALFA*B1-B2
       B2=B1
    4  B1=B0
       Q=Y*(B0-H*B2)
       B0=X-PI3
       B1=PI1*SQRT(R)*(Q*COS(B0)+P*SIN(B0))
      ENDIF
      IF(ENAME .EQ. ' BESY1')  THEN
         BESY1=ROUND(B1)
      ELSE
         DBESY1=B1
      ENDIF
      RETURN
  100 FORMAT(7X,A6,' ... NON-POSITIVE ARGUMENT X = ',E16.6)
      END