]>
Commit | Line | Data |
---|---|---|
0795afa3 | 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 |