C***********************************************************
C   THIS SUBROUTINE PERFORMS NUMERICAL INTEGRATION OF ANY SUITABLE
C   FUNCTION.  THE NEWTON-COTES METHOD IS USED.
C
C   DEFINITIONS AND EXPLAINATIONS:
C   F(X) = THE FUNCTION, DEFINED ON A CLOSED INTERVAL (A,B).
C          A=LOWER LIMIT OF INTERVAL
C          B=UPPER LIMIT OF INTERVAL
C   W    = B-A = INTERVAL WIDTH.
C   IT IS ASSUME F(X) IS BOUNDED ON THIS INTERVAL AND HAS FINITE 6-TH
C   DERIVATIVE.
C   THE INTERVAL IS PARTITIONED INTO N SUBINTERVALS OF EQUAL LENGTH
C   N    = NUMBER OF SUBINTERVALS. N MUST SATISFY N=3*4*K=12*K,
C          WHERE K IS AN INTEGER. USEFUL VALUES OF K ARE 10,50,100
C          GIVING N=120,600,1200 RESPECTIVELY.
C   F0   = F(A), THE FUNCTION AT X=A.
C   X(I) = A+I*W/N=THE PARTITION POINTS OF (A,B), I=1,...,N.
C   FUNC(I) = F(X(I))=THE FUNCTION VALUES AT THE PARTITION POINTS.
C
C   CALLING ARGUMENTS
C   N      CHOSEN BY THE USER, MUST SATISFY N=12*K.
C   FUNC(N)  COMPUTED BY THE USER, EXTERNAL TO THE ROUTINE.
C   F0     COMPUTED BY THE USER.
C   W
C   XI     THE VALUE OF THE INTEGRAL OF F(X) OVER (A,B)
C   IDIG   THE NUMBER OF DECIMAL DIGITS OF PROBABLE ACCURACY OF XI.
C          IF IDIG IS TOO LOW, INCREASE THE VALUE OF N.
C
C   OTHER NOTATION:
C   XNCI   = I-TH ORDER NEWTON-COTES INTEGRAL, I=1,2,3,4.
C            XNC1 = TRAPEZOIDAL RULE.
C            XNC2 = SIMPSONS RULE.
C            THESE FIRST TWO ARE NOT COMPUTED, BUT THEIR CODE IS IN
C            THE PROGRAM
C   XH     = W/N, THE SIZE OF THE PARTITION.
C   SUMI   = INTERMEDIATE SUMS OF THE FUNC(J) VALUES, I=1,2,3,4.
C***********************************************************
      SUBROUTINE NCNI (N,FUNC,F0,W,XI,IDIG)
C***********************************************************************
C     DECLARATIONS
C***********************************************************************
      DOUBLE PRECISION FUNC, XI, XJ, SUM1, SUM2, SUM3, SUM4, FUNCI, F0
      DOUBLE PRECISION XNC3, XNC4, RR, QQ, R, SF0N, XH
      DIMENSION FUNC(N)
      XI=DBLE(FLOAT(N))/12.0
      I=IDINT(XI)
      XJ=DBLE(FLOAT(I))
      IF (XJ.NE.XI) RETURN
      SUM1=FUNC(1)+FUNC(2)+FUNC(3)
      SUM2=FUNC(2)
      SUM3=0.0
      SUM4=0.0
      N3=N-3
      N4=N-4
      DO 10 I=3,N3,3
   10 SUM3=SUM3+FUNC(I)
      DO 20 I=4,N4,4
      FUNCI=FUNC(I)
      SUM1=SUM1+FUNCI+FUNC(I+1)+FUNC(I+2)+FUNC(I+3)
      SUM2=SUM2+FUNCI+FUNC(I+2)
   20 SUM4=SUM4+FUNCI
      SF0N=F0+FUNC(N)
      XH=DBLE(W)/DBLE(FLOAT(N))
C
C     XNC1=(SF0N+2*SUM1)*XH/2.0
C     XNC2=(SF0N+4*SUM1-2*SUM2)*XH/3.0
C
      XNC3=(SF0N+3.*SUM1-SUM3)*XH*3.0/8.0
      XNC4=(7.*SF0N+32.*SUM1-20.*SUM2+2.*SUM4)*XH*2.0/45.0
      XI=XNC4
      QQ=DABS(XI)
      RR=DABS(XNC3)
      IF (QQ.GT.0.00001) R=XNC3/XI
      IF ((QQ.LT.0.00001).AND.(RR.LT.0.00001)) R=(1.+XNC3)/(1.+XI)
      IF (R.LT.1.0) R=1.0/R
      TEST=R-1.0
      DO 30 IDIG=1,15
      IF (TEST.GT.1.0) GO TO 40
   30 TEST=TEST*10.0
   40 RETURN
      END