pythiaeRHIC/pythia/pygaus.f

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