*
* $Id: dgauss.F,v 1.1.1.1 1996/02/15 17:48:17 mclareni Exp $
*
* $Log: dgauss.F,v $
* Revision 1.1.1.1  1996/02/15 17:48:17  mclareni
* Kernlib
*
*
#include "kernnum/pilot.h"
      DOUBLE PRECISION FUNCTION DGAUSS(F,A,B,EPS)
      DOUBLE PRECISION F,A,B,EPS
      DOUBLE PRECISION W(12),X(12),AA,BB,C1,C2,U,S8,S16,CONST
      LOGICAL MFLAG,RFLAG
      EXTERNAL F
C
C     ******************************************************************
C
C     ADAPTIVE DOUBLE PRECISION GAUSSIAN QUADRATURE.
C
C     DGAUSS IS SET EQUAL TO THE APPROXIMATE VALUE OF THE INTEGRAL OF
C     THE FUNCTION F OVER THE INTERVAL (A,B), WITH ACCURACY PARAMETER
C     EPS.
C
C     ******************************************************************
C
#if defined(CERNLIB_NUMHIPRE)
C
      DATA W / 0.10122 85362 90376 25915 25313 543D0,
     1         0.22238 10344 53374 47054 43559 944D0,
     2         0.31370 66458 77887 28733 79622 020D0,
     3         0.36268 37833 78361 98296 51504 493D0,
     4         0.27152 45941 17540 94851 78057 246D-1,
     5         0.62253 52393 86478 92862 84383 699D-1,
     6         0.95158 51168 24927 84809 92510 760D-1,
     7         0.12462 89712 55533 87205 24762 822D0,
     8         0.14959 59888 16576 73208 15017 305D0,
     9         0.16915 65193 95002 53818 93120 790D0,
     A         0.18260 34150 44923 58886 67636 680D0,
     B         0.18945 06104 55068 49628 53967 232D0/
C
      DATA X / 0.96028 98564 97536 23168 35608 686D0,
     1         0.79666 64774 13626 73959 15539 365D0,
     2         0.52553 24099 16328 98581 77390 492D0,
     3         0.18343 46424 95649 80493 94761 424D0,
     4         0.98940 09349 91649 93259 61541 735D0,
     5         0.94457 50230 73232 57607 79884 155D0,
     6         0.86563 12023 87831 74388 04678 977D0,
     7         0.75540 44083 55003 03389 51011 948D0,
     8         0.61787 62444 02643 74844 66717 640D0,
     9         0.45801 67776 57227 38634 24194 430D0,
     A         0.28160 35507 79258 91323 04605 015D0,
     B         0.95012 50983 76374 40185 31933 543D-1/
#endif
#if defined(CERNLIB_NUMLOPRE)
      DATA W / 0.10122 85362 90376 259D0,
     1         0.22238 10344 53374 471D0,
     2         0.31370 66458 77887 287D0,
     3         0.36268 37833 78361 983D0,
     4         0.27152 45941 17540 949D-1,
     5         0.62253 52393 86478 929D-1,
     6         0.95158 51168 24927 848D-1,
     7         0.12462 89712 55533 872D0,
     8         0.14959 59888 16576 732D0,
     9         0.16915 65193 95002 538D0,
     A         0.18260 34150 44923 589D0,
     B         0.18945 06104 55068 496D0/
 
      DATA X / 0.96028 98564 97536 232D0,
     1         0.79666 64774 13626 740D0,
     2         0.52553 24099 16328 986D0,
     3         0.18343 46424 95649 805D0,
     4         0.98940 09349 91649 933D0,
     5         0.94457 50230 73232 576D0,
     6         0.86563 12023 87831 744D0,
     7         0.75540 44083 55003 034D0,
     8         0.61787 62444 02643 748D0,
     9         0.45801 67776 57227 386D0,
     A         0.28160 35507 79258 913D0,
     B         0.95012 50983 76374 402D-1/
#endif
C
C     ******************************************************************
C
C  START.
      DGAUSS=0.0D0
      IF(B.EQ.A) RETURN
      CONST=0.005D0/(B-A)
      BB=A
C
C  COMPUTATIONAL LOOP.
    1 AA=BB
      BB=B
    2    C1=0.5D0*(BB+AA)
         C2=0.5D0*(BB-AA)
         S8=0.0D0
         DO 3 I=1,4
            U=C2*X(I)
            S8=S8+W(I)*(F(C1+U)+F(C1-U))
    3    CONTINUE
         S8=C2*S8
         S16=0.0D0
         DO 4 I=5,12
            U=C2*X(I)
            S16=S16+W(I)*(F(C1+U)+F(C1-U))
    4    CONTINUE
         S16=C2*S16
         IF( ABS(S16-S8) .LE. EPS*(1.+ABS(S16)) ) GO TO 5
         BB=C1
         IF( 1.D0+ABS(CONST*C2) .NE. 1.D0) GO TO 2
      DGAUSS=0.0D0
      CALL KERMTR('D103.1',LGFILE,MFLAG,RFLAG)
      IF(MFLAG) THEN
         IF(LGFILE.EQ.0) THEN
            WRITE(*,6)
         ELSE
            WRITE(LGFILE,6)
         ENDIF
      ENDIF
      IF(.NOT. RFLAG) CALL ABEND
      RETURN
    5 DGAUSS=DGAUSS+S16
      IF(BB.NE.B) GO TO 1
      RETURN
C
    6 FORMAT( 4X, 'FUNCTION DGAUSS ... TOO HIGH ACCURACY REQUIRED')
      END