PROGRAM BESSELC
C======================================================================
C           BESSEL FUNCTION GENERATOR FOR HELMHOLTZ EQUATION
C        THIS PROGRAM COMPUTES BESSEL AND NEUMANN FUNCTIONS OF 
C                         J0, N0, J1, AND N1.
C             CODED BY EIJI FUKUMORI    1996    JAPAN
C======================================================================
      IMPLICIT REAL*8 ( A-H , O-Z )
      COMPLEX*16 H0, H1
      PI = 4.D0 * DATAN( 1.D0 )
      OPEN ( 1, FILE='BESSEL.FUN', STATUS = 'UNKNOWN' )
      WRITE(1,'(9X,1X,1HZ,9X,2HJ0,9X,2HN0,9X,2HJ1,9X,2HN1)')
      DO I = 1 , 150
      X = DFLOAT(I)/10.D0
      CALL BESSEL ( PI, X , H0, H1 )
      WRITE(1,'(5F11.6)') X,H0, H1
      END DO
      CLOSE (1)
      STOP
      END
C
C
      SUBROUTINE BESSEL ( PI, X , H0, H1 )
      IMPLICIT REAL*8 ( A-H , O-Z )
      COMPLEX*16 H0, H1
      PARAMETER ( EULER = 0.5772156649015D0, ALLOWANC=1.D-16 )
C--------- BJ0, BJ1
      X22 = (X/2.D0) **2
      PRODUCT = 1.D0
      SIGNCH = 1.D0
      FACT = 1.D0
      SUM0 = 0.D0
      SUM1 = 0.D0
      M = 1
      CHECK = 1.D0
      DO WHILE ( DABS(CHECK) .GT. ALLOWANC )
      FACT = FACT * DFLOAT(M)
      PRODUCT = PRODUCT * X22
      SIGNCH = -SIGNCH
      DELTA0 = SIGNCH*PRODUCT/(FACT*FACT)
      SUM0 = SUM0 + DELTA0
      IF ( SUM0 .EQ. 0.D0 ) THEN
      CHECK = DELTA0
      ELSE
      CHECK = DELTA0 / SUM0
      END IF
      SUM1 = SUM1 + SIGNCH*PRODUCT/(FACT*FACT*DFLOAT(M+1))
      M = M + 1
      IF ( M .GT. 1000 ) STOP
      END DO
      BJ0 = 1.D0 + SUM0
      BJ1 = (X/2.D0) * ( 1.D0 + SUM1 )
C
C-------- BN0, BN1
      PRODUCT = 1.D0
      SIGNCH = 1.D0
      FACT = 1.D0
      SUM0 = 0.D0
      SUM1 = 0.D0
      M = 1
      CHECK = 1.D0
      DO WHILE ( DABS(CHECK) .GT. ALLOWANC )
      FACT = FACT * DFLOAT(M)
      PRODUCT = PRODUCT * X22
      SIGNCH = -SIGNCH
C..............SUM OF 1/1 + 1/2 + 1/3 +.......+ 1/M
      SUMM = 0.D0
      DO I = 1 , M
      SUMM = SUMM + 1.D0 / DFLOAT(I)
      END DO
      DELTA0 = SIGNCH*PRODUCT/(FACT*FACT) * 2.D0*SUMM
      SUM0 = SUM0 + DELTA0
      IF ( SUM0 .EQ. 0.D0 ) THEN
      CHECK = DELTA0
      ELSE
      CHECK = DELTA0 / SUM0
      END IF
      SUM1 = SUM1 + SIGNCH*PRODUCT/(FACT*FACT*DFLOAT(M+1)) *
     *              (2.D0*SUMM+1.D0/DFLOAT(M+1))
      M = M + 1
      IF ( M .GT. 1000 ) STOP
      END DO
      BN0 = (2.D0/PI)*BJ0*(EULER+DLOG(X/2.D0)) - 1.D0/PI*SUM0
      BN1 = (2.D0/PI)*BJ1*(EULER+DLOG(X/2.D0)) - X/(2.D0*PI) -
     *      1.D0/PI*(X/2.D0)*SUM1 - 2.D0/(PI*X)
C----------- ASSEMBLE THEM INTO COMPLEX NUMBERS
      H0 = DCMPLX(BJ0,BN0)
      H1 = DCMPLX(BJ1,BN1)
      RETURN
      END