pygau2.f

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pygau2.f

 
C*********************************************************************
 
C...PYGAU2
C...Integration by adaptive Gaussian quadrature.
C...Adapted from the CERNLIB DGAUSS routine by K.S. Kolbig.
C...Carbon copy of PYGAUS, but avoids having to use it recursively.
 
      FUNCTION pygau2(F, A, B, EPS)
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      IMPLICIT INTEGER(i-n)
      INTEGER pyk,pychge,pycomp
 
C...Local declarations.
      EXTERNAL f
      DOUBLE PRECISION f,w(12), x(12)
      DATA x( 1) /9.6028985649753623d-1/, w( 1) /1.0122853629037626d-1/
      DATA x( 2) /7.9666647741362674d-1/, w( 2) /2.2238103445337447d-1/
      DATA x( 3) /5.2553240991632899d-1/, w( 3) /3.1370664587788729d-1/
      DATA x( 4) /1.8343464249564980d-1/, w( 4) /3.6268378337836198d-1/
      DATA x( 5) /9.8940093499164993d-1/, w( 5) /2.7152459411754095d-2/
      DATA x( 6) /9.4457502307323258d-1/, w( 6) /6.2253523938647893d-2/
      DATA x( 7) /8.6563120238783174d-1/, w( 7) /9.5158511682492785d-2/
      DATA x( 8) /7.5540440835500303d-1/, w( 8) /1.2462897125553387d-1/
      DATA x( 9) /6.1787624440264375d-1/, w( 9) /1.4959598881657673d-1/
      DATA x(10) /4.5801677765722739d-1/, w(10) /1.6915651939500254d-1/
      DATA x(11) /2.8160355077925891d-1/, w(11) /1.8260341504492359d-1/
      DATA x(12) /9.5012509837637440d-2/, w(12) /1.8945061045506850d-1/
 
C...The Gaussian quadrature algorithm.
      h = 0d0
      IF(b .EQ. a) goto 140
      const = 5d-3 / abs(b-a)
      bb = a
  100 CONTINUE
      aa = bb
      bb = b
  110 CONTINUE
      c1 = 0.5d0*(bb+aa)
      c2 = 0.5d0*(bb-aa)
      s8 = 0d0
      DO 120 i = 1, 4
        u = c2*x(i)
        s8 = s8 + w(i) * (f(c1+u) + f(c1-u))
  120 CONTINUE
      s16 = 0d0
      DO 130 i = 5, 12
        u = c2*x(i)
        s16 = s16 + w(i) * (f(c1+u) + f(c1-u))
  130 CONTINUE
      s16 = c2*s16
      IF(dabs(s16-c2*s8) .LE. eps*(1d0+dabs(s16))) THEN
        h = h + s16
        IF(bb .NE. b) goto 100
      ELSE
        bb = c1
        IF(1d0 + const*abs(c2) .NE. 1d0) goto 110
        h = 0d0
        CALL pyerrm(18,'(PYGAU2:) too high accuracy required')
        goto 140
      ENDIF
  140 CONTINUE
      pygau2 = h
 
      RETURN
      END