! ---------------------------------------------------------------------
! Program to integrate a user-defined function f(x) from x1 to x2 by
! the QANC8 subroutine with control of absolute and relative precisions
! ---------------------------------------------------------------------
! SAMPLE RUN:
!
!  (Integrate function cos(x) - 2 sin(x) from x=0 to x=1,
!   with a precision <= 1e-10).
! 
!  Integral Value = -7.792440345582408E-002
!
!  Estimated error =  1.041851518709824E-017
!
!  Error code =  0.000000000000000E+000
!
!  Number of function evaluations =          33
!
! ---------------------------------------------------------------------
! Reference: From Numath Library By Tuan Dang Trong in Fortran 77
!            [BIBLI 18].
!
!                               F90 Release 1.0 By J-P Moreau, Paris
!                                        (www.jpmoreau.fr)
! ---------------------------------------------------------------------
PROGRAM TQANC8

REAL*8  AERROR,CODE,ERROR,RERROR,X1,X2,VAL

  X1=0.D0
  X2=1.D0
  AERROR=1.D-9
  RERROR=1.D-10

  CALL QANC8(X1,X2,AERROR,RERROR,VAL,ERROR,NBRF,CODE)

  print *,' '
  print *,' Integral Value =', VAL
  print *,' '
  print *,' Estimated error =', ERROR
  print *,' '
  print *,' Error Code =', CODE
  print *,' '
  print *,' Number of function evaluations =', NBRF
  print *,' '
     
END

REAL*8 FUNCTION FCT(X)
REAL*8 X
  
  FCT=DCOS(X)-2.D0*DSIN(X)

RETURN
END

! ----------------------------------------------------------------------
      SUBROUTINE QANC8 (A,B,AERR,RERR,RES,ERR,NBF,FLG)
!
!     INTEGRATE A REAL FUNCTION FCT(X) FROM X=A TO X=B, WITH
!     GIVEN ABSOLUTE AND RELATIVE PRECISIONS, AERR, RERR. 

!     INPUTS:
!     FCT     EXTERNAL USER-DEFINED FUNCTION FOR ANY X VALUE
!             IN INTERVAL (A,B)
!     A,B     LIMITS OF INTERVAL
!     AERR,RERR   RESPECTIVELY ABSOLUTE ERROR AND RELATIVE ERROR
!                 REQUIRED BY USER
!
!     OUTPUTS:
!     RES     VALUE OF INTEGRAL
!     ERR     ESTIMATED ERROR
!     NBF     NUMBER OF NECESSARY FCT(X) EVALUATIONS
!     FLG     INDICATOR
!             = 0.0       CORRECT RESULT
!             = NNN.RRR   NO CONVERGENCE DU TO A SINGULARITY.
!             THE SINGULAR POINT ABCISSA IS GIVEN BY FORMULA:
!             XS = B-.RRR*(B-A)

!     REFERENCE :
!     FORSYTHE,G.E. (1977) COMPUTER METHODS FOR MATHEMATICAL
!     COMPUTATIONS. PRENTICE-HALL, INC.
! -----------------------------------------------------------------------
      IMPLICIT REAL *8 (A-H,O-Z)
      DIMENSION QR(31),F(16),X(16),FS(8,30),XS(8,30)
      LMIN = 1
      LMAX = 30
      LOUT = 6
      NMAX = 5000
      NFIN = NMAX-8*(LMAX-LOUT+2**(LOUT+1))
      W0   =   3956.D0/14175.D0
      W1   =  23552.D0/14175.D0
      W2   =  -3712.D0/14175.D0
      W3   =  41984.D0/14175.D0
      W4   = -18160.D0/14175.D0
      FLG  = 0.D0
      RES  = 0.D0
      COR  = 0.D0
      ERR  = 0.D0
      SUM  = 0.D0
      NBF  = 0
      IF (A.EQ.B) RETURN
      L = 0
      NIM = 1
      X0  = A
      X(16) = B
      QP  = 0.D0
      F0   = FCT(X0)
      PAS1  = (B-A)/16.D0
      X(8)  = (X0+X(16))*.5D0
      X(4)  = (X0+X(8))*.5D0
      X(12) = (X(8)+X(16))*.5D0
      X(2)  = (X0+X(4))*.5D0
      X(6)  = (X(4)+X(8))*.5D0
      X(10) = (X(8)+X(12))*.5D0
      X(14) = (X(12)+X(16))*.5D0
      DO 25 J = 2,16,2
      F(J) = FCT(X(J))
   25 CONTINUE
      NBF = 9
   30 X(1)  = (X0+X(2))*.5D0
      F(1) = FCT(X(1))
      DO 35 J = 3,15,2
      X(J)  = (X(J-1)+X(J+1))*.5D0
   35 F(J) = FCT(X(J))
      NBF = NBF+8
      PAS = (X(16)-X0)/16.D0
      QL  = (W0*(F0+F(8))+W1*(F(1)+F(7))+W2*(F(2)+F(6))        &
                                 +W3*(F(3)+F(5))+W4*F(4))*PAS
      QR(L+1) = (W0*(F(8)+F(16))+W1*(F(9)+F(15))               &
             +W2*(F(10)+F(14))+W3*(F(11)+F(13))+W4*F(12))*PAS
      QN = QL + QR(L+1)
      QD = QN - QP
      SUM = SUM + QD
      ERR1 = DABS(QD)/1023.D0
      TOL1 = DMAX1(AERR,RERR*DABS(SUM))*(PAS/PAS1)
      IF (L.LT.LMIN) GO TO 50
      IF (L.GE.LMAX) GO TO 62
      IF (NBF.GT.NFIN) GO TO 60
      IF (ERR1.LE.TOL1) GO TO 70
   50 NIM = 2*NIM
      L = L+1
      DO 52 I = 1,8
      FS(I,L) = F(I+8)
      XS(I,L) = X(I+8)
   52 CONTINUE
      QP = QL
      DO 55 I = 1,8
      F(18-2*I) = F(9-I)
      X(18-2*I) = X(9-I)
   55 CONTINUE
      GO TO 30
   60 NFIN = 2*NFIN
      LMAX = LOUT
      FLG = FLG + (B-X0)/(B-A)
      GO TO 70
   62 FLG = FLG + 1.D0
   70 RES = RES + QN
      ERR = ERR + ERR1
      COR = COR + QD/1023.D0
   72 IF (NIM.EQ.2*(NIM/2)) GO TO 75
      NIM = NIM/2
      L = L-1
      GO TO 72
   75 NIM = NIM+1
      IF (L.LE.0) GO TO 80
      QP = QR(L)
      X0 = X(16)
      F0 = F(16)
      DO 78 I = 1,8
      F(2*I) = FS(I,L)
      X(2*I) = XS(I,L)
   78 CONTINUE
      GO TO 30
   80 RES = RES + COR
      IF (ERR.EQ.0.D0) RETURN
   82 TEMP = DABS(RES) + ERR
      IF (TEMP.NE.DABS(RES)) RETURN
      ERR = 2.D0*ERR
      GO TO 82
      END

! end of file tqanc8.f90