ISAJET/isasusy/ssdint.F

   1 #include "isajet/pilot.h"
   2       DOUBLE PRECISION FUNCTION SSDINT(XL,F,XR)
   3 C-----------------------------------------------------------------------
   4 C          Integrate F over (XL,XR) using adaptive Gaussian quadrature.
   5 C          TOLABS = 1e-35: absolute error for convergence.
   6 C          TOLREL = 5e-5:  relative error for convergence.
   7 C          TOLBIN = 1e-3:  relative bin size limit. Set contribution to
   8 C                          zero if bin falls below this.
   9 C     Note quadrature constants R and W have been converted to explicit
  10 C     double precision (.xxxxxDxx) form.
  11 C
  12 C     Bisset's XINTH
  13 C-----------------------------------------------------------------------
  14 #if defined(CERNLIB_IMPNONE)
  15       IMPLICIT NONE
  16 #endif
  17 #include "isajet/sslun.inc"
  18       EXTERNAL F
  19       INTEGER NMAX
  20       DOUBLE PRECISION TOLABS,TOLREL,TOLBIN,XMIN,XLIMS(100)
  21       DOUBLE PRECISION R(93),W(93)
  22       INTEGER PTR(4),NORD(4)
  23       INTEGER ICOUNT,IBAD
  24       DOUBLE PRECISION XL,XR,F
  25       DOUBLE PRECISION AA,BB,TVAL,VAL,TOL
  26       INTEGER NLIMS,I,J,KKK
  27 C
  28       DATA PTR,NORD/4,10,22,46,  6,12,24,48/
  29       DATA (R(KKK),KKK=1,48)/
  30      $ .2386191860D0,.6612093865D0,.9324695142D0,.1252334085D0,
  31      $ .3678314990D0,.5873179543D0,.7699026742D0,.9041172563D0,
  32      $ .9815606342D0,.0640568929D0,.1911188675D0,.3150426797D0,
  33      $ .4337935076D0,.5454214714D0,.6480936519D0,.7401241916D0,
  34      $ .8200019860D0,.8864155270D0,.9382745520D0,.9747285560D0,
  35      $ .9951872200D0,.0323801710D0,.0970046992D0,.1612223561D0,
  36      $ .2247637903D0,.2873624873D0,.3487558863D0,.4086864820D0,
  37      $ .4669029048D0,.5231609747D0,.5772247261D0,.6288673968D0,
  38      $ .6778723796D0,.7240341309D0,.7671590325D0,.8070662040D0,
  39      $ .8435882616D0,.8765720203D0,.9058791367D0,.9313866907D0,
  40      $ .9529877032D0,.9705915925D0,.9841245837D0,.9935301723D0,
  41      $ .9987710073D0,.0162767488D0,.0488129851D0,.0812974955D0/
  42       DATA (R(KKK),KKK=49,93)/
  43      $ .1136958501D0,.1459737146D0,.1780968824D0,.2100313105D0,
  44      $ .2417431561D0,.2731988126D0,.3043649444D0,.3352085229D0,
  45      $ .3656968614D0,.3957976498D0,.4254789884D0,.4547094222D0,
  46      $ .4834579739D0,.5116941772D0,.5393881083D0,.5665104186D0,
  47      $ .5930323648D0,.6189258401D0,.6441634037D0,.6687183100D0,
  48      $ .6925645366D0,.7156768123D0,.7380306437D0,.7596023411D0,
  49      $ .7803690438D0,.8003087441D0,.8194003107D0,.8376235112D0,
  50      $ .8549590334D0,.8713885059D0,.8868945174D0,.9014606353D0,
  51      $ .9150714231D0,.9277124567D0,.9393703398D0,.9500327178D0,
  52      $ .9596882914D0,.9683268285D0,.9759391746D0,.9825172636D0,
  53      $ .9880541263D0,.9925439003D0,.9959818430D0,.9983643759D0,
  54      $ .9996895039/
  55         DATA (W(KKK),KKK=1,48)/
  56      $ .4679139346D0,.3607615730D0,.1713244924D0,.2491470458D0,
  57      $ .2334925365D0,.2031674267D0,.1600783285D0,.1069393260D0,
  58      $ .0471753364D0,.1279381953D0,.1258374563D0,.1216704729D0,
  59      $ .1155056681D0,.1074442701D0,.0976186521D0,.0861901615D0,
  60      $ .0733464814D0,.0592985849D0,.0442774388D0,.0285313886D0,
  61      $ .0123412298D0,.0647376968D0,.0644661644D0,.0639242386D0,
  62      $ .0631141923D0,.0620394232D0,.0607044392D0,.0591148397D0,
  63      $ .0572772921D0,.0551995037D0,.0528901894D0,.0503590356D0,
  64      $ .0476166585D0,.0446745609D0,.0415450829D0,.0382413511D0,
  65      $ .0347772226D0,.0311672278D0,.0274265097D0,.0235707608D0,
  66      $ .0196161605D0,.0155793157D0,.0114772346D0,.0073275539D0,
  67      $ .0031533461D0,.0325506145D0,.0325161187D0,.0324471637D0/
  68       DATA (W(KKK),KKK=49,93)/
  69      $ .0323438226D0,.0322062048D0,.0320344562D0,.0318287589D0,
  70      $ .0315893308D0,.0313164256D0,.0310103326D0,.0306713761D0,
  71      $ .0302999154D0,.0298963441D0,.0294610900D0,.0289946142D0,
  72      $ .0284974111D0,.0279700076D0,.0274129627D0,.0268268667D0,
  73      $ .0262123407D0,.0255700360D0,.0249006332D0,.0242048418D0,
  74      $ .0234833991D0,.0227370697D0,.0219666444D0,.0211729399D0,
  75      $ .0203567972D0,.0195190811D0,.0186606796D0,.0177825023D0,
  76      $ .0168854799D0,.0159705629D0,.0150387210D0,.0140909418D0,
  77      $ .0131282296D0,.0121516047D0,.0111621020D0,.0101607705D0,
  78      $ .0091486712D0,.0081268769D0,.0070964708D0,.0060585455D0,
  79      $ .0050142027D0,.0039645543D0,.0029107318D0,.0018539608D0,
  80      $ .0007967921/
  81 C
  82       DATA TOLABS,TOLREL,TOLBIN,NMAX/1.D-35,5.D-5,1D-5,100/
  83 C
  84       SSDINT=0
  85       NLIMS=2
  86       XLIMS(1)=XL
  87       XLIMS(2)=XR
  88       ICOUNT=0
  89       IBAD=0
  90       XMIN=TOLBIN*ABS(XR-XL)
  91 C
  92 10    AA=(XLIMS(NLIMS)-XLIMS(NLIMS-1))/2
  93       BB=(XLIMS(NLIMS)+XLIMS(NLIMS-1))/2
  94       TVAL=0
  95       DO 15 I=1,3
  96 15    TVAL=TVAL+W(I)*(F(BB+AA*R(I))+F(BB-AA*R(I)))
  97       TVAL=TVAL*AA
  98       DO 25 J=1,4
  99         VAL=0
 100         DO 20 I=PTR(J),PTR(J)-1+NORD(J)
 101           ICOUNT=ICOUNT+1
 102           IF(ICOUNT.GT.1E5) THEN
 103             WRITE(LOUT,*) 'ERROR IN SSDINT: SET SSDINT = 0'
 104             SSDINT=0.
 105             RETURN
 106           ENDIF
 107 20      VAL=VAL+W(I)*(F(BB+AA*R(I))+F(BB-AA*R(I)))
 108         VAL=VAL*AA
 109         TOL=MAX(TOLABS,TOLREL*ABS(VAL))
 110         IF(ABS(TVAL-VAL).LT.TOL) THEN
 111           SSDINT=SSDINT+VAL
 112           NLIMS=NLIMS-2
 113           IF (NLIMS.NE.0) GO TO 10
 114           GO TO 999
 115         ELSEIF(ABS(AA).LT.XMIN.AND.J.EQ.4) THEN
 116 C           Bin is too small -- set VAL = 0. -- FEP
 117           IBAD=IBAD+1
 118           NLIMS=NLIMS-2
 119           IF (NLIMS.NE.0) GO TO 10
 120           GO TO 999
 121         ENDIF
 122 25    TVAL=VAL
 123       IF(NMAX.EQ.2) THEN
 124         SSDINT=VAL
 125         GO TO 999
 126       ENDIF
 127       IF(NLIMS.GT.(NMAX-2)) THEN
 128         WRITE(LOUT,50) SSDINT,NMAX,BB-AA,BB+AA
 129 50      FORMAT (' ERROR IN SSDINT, SSDINT,NMAX,XL,XR=',G15.7,I5,2G15.7)
 130         RETURN
 131       ENDIF
 132       XLIMS(NLIMS+1)=BB
 133       XLIMS(NLIMS+2)=BB+AA
 134       XLIMS(NLIMS)=BB
 135       NLIMS=NLIMS+2
 136       GO TO 10
 137 C
 138 999   IF(IBAD.GT.0) THEN
 139         WRITE(LOUT,*) 'WARNING IN SSXINT: BAD CONVERGENCE FOR ',IBAD,
 140      $  ' INTERVALS.'
 141       ENDIF
 142       RETURN
 143       END