C     *****************************************************************
C     INTEGRATION BY NEWTON-COTES RULE (N=6)
C
C     A, B = ARE THE INTEGRATION LIMITS
C     SUM  = THE RESULT OF THE INTEGRATION
C     THE FUNCTION TO BE INTEGRATED MUST BE DEFINED IN
C     A SUBPROGRAM FUNCTION FUN(ARGUMENTS)
C
C     SPECIFICALLY, THIS PROGRAM EVALUATES BESSEL FUNCTION
C     JN(X), OF ORDER N AND ARGUMENT X
C     *****************************************************************

      DATA PIE/3.1415927/

      A=0.0
      B=PIE
      N=0
      DO 20 I=1,20
         X=0.1*FLOAT(I)
         CALL INTEGRATION(A,B,N,X,SUM)
         WRITE(6,10) X,SUM
 10      FORMAT(2X,'X = ',F6.3,2X,'J0 = ',E20.12,/)
 20   CONTINUE
      STOP
      END

C     *****************************************************************
C     SUBROUTINE FOR INTEGRATION USING NEWTON-COTES RULE (N=6)
C     A IS THE LOWER LIMIT OF INTEGRATION
C     B IS THE UPPER LIMIT OF INTEGRATION
C     H IS THE INTERVAL SIZE
C     *****************************************************************
      SUBROUTINE INTEGRATION(AA,BB,NO,X,SUM)
      IMPLICIT REAL(A-H,O-Z)
      DIMENSION C(0:6)
      DATA NC /840/
      DATA ( C(I),I=0,6 )/41.,216.,27.,272.,27.,216.,41./

C     START COMPUTATION

      N=6
      M=240    ! M MUST BE A MULTIPLE OF N
      H=(BB-AA)/FLOAT(M)
      SUM=0.0
      NN=M/N
      T=AA
      DO 20 I=1,NN
         DO 10 J=0,N
            A=C(J)*FUN(T,NO,X)
            SUM = SUM + A
            T = T + H
 10      CONTINUE
         T = T - H
 20   CONTINUE
      SUM = SUM*H*FLOAT(N)/FLOAT(NC)
      RETURN
      END

C     *****************************************************************
C     THE FUNCTION TO BE INTEGRATED
C     *****************************************************************
      FUNCTION FUN(THETA,N,X)
      DATA PIE/3.1415927/

      FUN = COS ( X*SIN(THETA) - FLOAT(N)*THETA )/PIE
      RETURN
      END