ISAJET/isasusy/ssxint.F

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