pythiaeRHIC/pythia/pygau2.f

0001
0002 C*********************************************************************
0003
0004 C...PYGAU2
0005 C...Integration by adaptive Gaussian quadrature.
0006 C...Adapted from the CERNLIB DGAUSS routine by K.S. Kolbig.
0007 C...Carbon copy of PYGAUS, but avoids having to use it recursively.
0008
0009       FUNCTION PYGAU2(F, A, B, EPS)
0010
0011 C...Double precision and integer declarations.
0012       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0013       IMPLICIT INTEGER(I-N)
0014       INTEGER PYK,PYCHGE,PYCOMP
0015
0016 C...Local declarations.
0017       EXTERNAL F
0018       DOUBLE PRECISION F,W(12), X(12)
0019       DATA X( 1) /9.6028985649753623D-1/, W( 1) /1.0122853629037626D-1/
0020       DATA X( 2) /7.9666647741362674D-1/, W( 2) /2.2238103445337447D-1/
0021       DATA X( 3) /5.2553240991632899D-1/, W( 3) /3.1370664587788729D-1/
0022       DATA X( 4) /1.8343464249564980D-1/, W( 4) /3.6268378337836198D-1/
0023       DATA X( 5) /9.8940093499164993D-1/, W( 5) /2.7152459411754095D-2/
0024       DATA X( 6) /9.4457502307323258D-1/, W( 6) /6.2253523938647893D-2/
0025       DATA X( 7) /8.6563120238783174D-1/, W( 7) /9.5158511682492785D-2/
0026       DATA X( 8) /7.5540440835500303D-1/, W( 8) /1.2462897125553387D-1/
0027       DATA X( 9) /6.1787624440264375D-1/, W( 9) /1.4959598881657673D-1/
0028       DATA X(10) /4.5801677765722739D-1/, W(10) /1.6915651939500254D-1/
0029       DATA X(11) /2.8160355077925891D-1/, W(11) /1.8260341504492359D-1/
0030       DATA X(12) /9.5012509837637440D-2/, W(12) /1.8945061045506850D-1/
0031
0032 C...The Gaussian quadrature algorithm.
0033       H = 0D0
0034       IF(B .EQ. A) GOTO 140
0035       CONST = 5D-3 / ABS(B-A)
0036       BB = A
0037   100 CONTINUE
0038       AA = BB
0039       BB = B
0040   110 CONTINUE
0041       C1 = 0.5D0*(BB+AA)
0042       C2 = 0.5D0*(BB-AA)
0043       S8 = 0D0
0044       DO 120 I = 1, 4
0045         U = C2*X(I)
0046         S8 = S8 + W(I) * (F(C1+U) + F(C1-U))
0047   120 CONTINUE
0048       S16 = 0D0
0049       DO 130 I = 5, 12
0050         U = C2*X(I)
0051         S16 = S16 + W(I) * (F(C1+U) + F(C1-U))
0052   130 CONTINUE
0053       S16 = C2*S16
0054       IF(DABS(S16-C2*S8) .LE. EPS*(1D0+DABS(S16))) THEN
0055         H = H + S16
0056         IF(BB .NE. B) GOTO 100
0057       ELSE
0058         BB = C1
0059         IF(1D0 + CONST*ABS(C2) .NE. 1D0) GOTO 110
0060         H = 0D0
0061         CALL PYERRM(18,'(PYGAU2:) too high accuracy required')
0062         GOTO 140
0063       ENDIF
0064   140 CONTINUE
0065       PYGAU2 = H
0066
0067       RETURN
0068       END