pygaus.f

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

 
C*********************************************************************
 
C...PYGAUS
C...Integration by adaptive Gaussian quadrature.
C...Adapted from the CERNLIB DGAUSS routine by K.S. Kolbig.
 
      FUNCTION pygaus(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,'(PYGAUS:) too high accuracy required')
        goto 140
      ENDIF
  140 CONTINUE
      pygaus = h
 
      RETURN
      END