C++ Interface to Tauola
gaus_integr.f
1* THIS IS ITERATIVE INTEGRATION PROCEDURE
2* ORIGINATES PROBABLY FROM CERN LIBRARY
3* IT SUBDIVIDES INEGRATION RANGE UNTIL REQUIRED PRECISION IS REACHED
4* PRECISION IS A DIFFERENCE FROM 8 AND 16 POINT GAUSS ITEGR. RESULT
5* EEPS POSITIVE TREATED AS ABSOLUTE PRECISION
6* EEPS NEGATIVE TREATED AS RELATIVE PRECISION
7 DOUBLE PRECISION FUNCTION gaus(F,A,B,EEPS)
8* *************************
9 IMPLICIT REAL*8(a-h,o-z)
10 EXTERNAL f
11 dimension w(12),x(12)
12 DATA const /1.0d-19/
13 DATA w
14 1/0.10122 85362 90376, 0.22238 10344 53374, 0.31370 66458 77887,
15 2 0.36268 37833 78362, 0.02715 24594 11754, 0.06225 35239 38648,
16 3 0.09515 85116 82493, 0.12462 89712 55534, 0.14959 59888 16577,
17 4 0.16915 65193 95003, 0.18260 34150 44924, 0.18945 06104 55069/
18 DATA x
19 1/0.96028 98564 97536, 0.79666 64774 13627, 0.52553 24099 16329,
20 2 0.18343 46424 95650, 0.98940 09349 91650, 0.94457 50230 73233,
21 3 0.86563 12023 87832, 0.75540 44083 55003, 0.61787 62444 02644,
22 4 0.45801 67776 57227, 0.28160 35507 79259, 0.09501 25098 37637/
23
24C-----------------------------------------------------------------------------
25C NEW INTEGRATION METHOD - constant step size (provided by fourth parameter)
26C comment out following 5 lines to return to old solutions of advantages too
27C-----------------------------------------------------------------------------
28 result=0
29 eps=0
30c CALL STEPGAUSS(F,A,B,6,RESULT,EPS)
31 CALL changegauss(f,a,b,2,result,eps)
32 gaus=result
33 RETURN
34C-----------------------------------------------------------------------------
35C PREVIOUS INTEGRATION METHOD - step size depends on EPSSQ
36C-----------------------------------------------------------------------------
37 eps=dabs(eeps)
38 delta=const*dabs(a-b)
39 gaus=0.d0
40 aa=a
41 5 y=b-aa
42 IF(dabs(y) .LE. delta) RETURN
43 2 bb=aa+y
44 c1=0.5d0*(aa+bb)
45 c2=c1-aa
46 s8=0.d0
47 s16=0.d0
48 DO 1 i=1,4
49 u=x(i)*c2
50 1 s8=s8+w(i)*(f(c1+u)+f(c1-u))
51 DO 3 i=5,12
52 u=x(i)*c2
53 3 s16=s16+w(i)*(f(c1+u)+f(c1-u))
54 s8=s8*c2
55 s16=s16*c2
56 IF(eeps.LT.0d0) THEN
57 IF(dabs(s16-s8) .GT. eps*dabs(s16)) GO TO 4
58 ELSE
59 IF(dabs(s16-s8) .GT. eps) GO TO 4
60 ENDIF
61 gaus=gaus+s16
62 aa=bb
63 GO TO 5
64 4 y=0.5d0*y
65 IF(dabs(y) .GT. delta) GOTO 2
66 print 7
67 gaus=0.d0
68 RETURN
69 7 FORMAT(1x,36hgaus ... too high accuracy required)
70 END
71
72
73 DOUBLE PRECISION FUNCTION gaus2(F,A,B,EEPS)
74* *************************
75 IMPLICIT REAL*8(a-h,o-z)
76 EXTERNAL f
77 dimension w(12),x(12)
78 DATA const /1.0d-19/
79 DATA w
80 1/0.10122 85362 90376, 0.22238 10344 53374, 0.31370 66458 77887,
81 2 0.36268 37833 78362, 0.02715 24594 11754, 0.06225 35239 38648,
82 3 0.09515 85116 82493, 0.12462 89712 55534, 0.14959 59888 16577,
83 4 0.16915 65193 95003, 0.18260 34150 44924, 0.18945 06104 55069/
84 DATA x
85 1/0.96028 98564 97536, 0.79666 64774 13627, 0.52553 24099 16329,
86 2 0.18343 46424 95650, 0.98940 09349 91650, 0.94457 50230 73233,
87 3 0.86563 12023 87832, 0.75540 44083 55003, 0.61787 62444 02644,
88 4 0.45801 67776 57227, 0.28160 35507 79259, 0.09501 25098 37637/
89
90C-----------------------------------------------------------------------------
91C NEW INTEGRATION METHOD - constant step size (provided by fourth parameter)
92C comment out following 5 lines to return to old solutions of advantages too
93C-----------------------------------------------------------------------------
94 result=0
95 eps=0
96c CALL STEPGAUSS2(F,A,B,6,RESULT,EPS)
97 CALL changegauss2(f,a,b,2,result,eps)
98 gaus2=result
99 RETURN
100C-----------------------------------------------------------------------------
101C PREVIOUS INTEGRATION METHOD - step size depends on EPSSQ
102C-----------------------------------------------------------------------------
103 eps=dabs(eeps)
104 delta=const*dabs(a-b)
105 gaus2=0.d0
106 aa=a
107 5 y=b-aa
108 IF(dabs(y) .LE. delta) RETURN
109 2 bb=aa+y
110 c1=0.5d0*(aa+bb)
111 c2=c1-aa
112 s8=0.d0
113 s16=0.d0
114 DO 1 i=1,4
115 u=x(i)*c2
116 1 s8=s8+w(i)*(f(c1+u)+f(c1-u))
117 DO 3 i=5,12
118 u=x(i)*c2
119 3 s16=s16+w(i)*(f(c1+u)+f(c1-u))
120 s8=s8*c2
121 s16=s16*c2
122 IF(eeps.LT.0d0) THEN
123 IF(dabs(s16-s8) .GT. eps*dabs(s16)) GO TO 4
124 ELSE
125 IF(dabs(s16-s8) .GT. eps) GO TO 4
126 ENDIF
127 gaus2=gaus2+s16
128 aa=bb
129 GO TO 5
130 4 y=0.5d0*y
131 IF(dabs(y) .GT. delta) GOTO 2
132 print 7
133 gaus2=0.d0
134 RETURN
135 7 FORMAT(1x,36hgaus2 ... too high accuracy required)
136 END
137
138 DOUBLE PRECISION FUNCTION gaus3(F,A,B,EEPS)
139* *************************
140 IMPLICIT REAL*8(a-h,o-z)
141 EXTERNAL f
142 dimension w(12),x(12)
143 DATA const /1.0d-19/
144 DATA w
145 1/0.10122 85362 90376, 0.22238 10344 53374, 0.31370 66458 77887,
146 2 0.36268 37833 78362, 0.02715 24594 11754, 0.06225 35239 38648,
147 3 0.09515 85116 82493, 0.12462 89712 55534, 0.14959 59888 16577,
148 4 0.16915 65193 95003, 0.18260 34150 44924, 0.18945 06104 55069/
149 DATA x
150 1/0.96028 98564 97536, 0.79666 64774 13627, 0.52553 24099 16329,
151 2 0.18343 46424 95650, 0.98940 09349 91650, 0.94457 50230 73233,
152 3 0.86563 12023 87832, 0.75540 44083 55003, 0.61787 62444 02644,
153 4 0.45801 67776 57227, 0.28160 35507 79259, 0.09501 25098 37637/
154
155C-----------------------------------------------------------------------------
156C NEW INTEGRATION METHOD - constant step size (provided by fourth parameter)
157C comment out following 5 lines to return to old solutions of advantages too
158C-----------------------------------------------------------------------------
159 result=0
160 eps=0
161c CALL STEPGAUSS3(F,A,B,6,RESULT,EPS)
162 CALL changegauss3(f,a,b,2,result,eps)
163 gaus3=result
164 RETURN
165C-----------------------------------------------------------------------------
166C PREVIOUS INTEGRATION METHOD - step size depends on EPSSQ
167C-----------------------------------------------------------------------------
168 eps=dabs(eeps)
169 delta=const*dabs(a-b)
170 gaus3=0.d0
171 aa=a
172 5 y=b-aa
173 IF(dabs(y) .LE. delta) RETURN
174 2 bb=aa+y
175 c1=0.5d0*(aa+bb)
176 c2=c1-aa
177 s8=0.d0
178 s16=0.d0
179 DO 1 i=1,4
180 u=x(i)*c2
181 1 s8=s8+w(i)*(f(c1+u)+f(c1-u))
182 DO 3 i=5,12
183 u=x(i)*c2
184 3 s16=s16+w(i)*(f(c1+u)+f(c1-u))
185 s8=s8*c2
186 s16=s16*c2
187 IF(eeps.LT.0d0) THEN
188 IF(dabs(s16-s8) .GT. eps*dabs(s16)) GO TO 4
189 ELSE
190 IF(dabs(s16-s8) .GT. eps) GO TO 4
191 ENDIF
192 gaus3=gaus3+s16
193 aa=bb
194 GO TO 5
195 4 y=0.5d0*y
196 IF(dabs(y) .GT. delta) GOTO 2
197 print 7
198 gaus3=0.d0
199 RETURN
200 7 FORMAT(1x,36hgaus3 ... too high accuracy required)
201 END
202
203C--------------------------------------------------------------------------------------------------
204C--------------------------------------------------------------------------------------------------
205C--------------------------------------------------------------------------------------------------
206C OBSOLETE 1.13 901108 3.40 CERN PROGRAM LIBRARY OBSOLETE PAM
207C CERN OBSOLETE PROGRAMS AND SUBROUTINES PAM-FILE.
208C MANY OF THE ROUTINES ARE FOR ONE FORTRAN DIALECT ONLY.
209C THE FLAG PATCH F77 SHOULD BE USED FOR THE FORTRAN 77 VERSIONS
210C SHOULD THEY EXIST - IF NOT FORTRAN 4 WILL BE OBTAINED.
211C THE CERN LIBRARIES GENLIB AND KERNLIB ARE USUALLY REQUIRED.
212C OBSOLETE ROUTINES OF THE PROGRAM LIBRARY ARE OFTEN REQUIRED.
213C
214 SUBROUTINE stepgauss(F,A,B,LAMBDA,RESULT,EPS)
215 implicit real*8 (a-h,o-z)
216 dimension w(12),x(12)
217 DATA w
218 1/0.10122 85362 90376, 0.22238 10344 53374, 0.31370 66458 77887,
219 2 0.36268 37833 78362, 0.02715 24594 11754, 0.06225 35239 38648,
220 3 0.09515 85116 82493, 0.12462 89712 55534, 0.14959 59888 16577,
221 4 0.16915 65193 95003, 0.18260 34150 44924, 0.18945 06104 55069/
222 DATA x
223 1/0.96028 98564 97536, 0.79666 64774 13627, 0.52553 24099 16329,
224 2 0.18343 46424 95650, 0.98940 09349 91650, 0.94457 50230 73233,
225 3 0.86563 12023 87832, 0.75540 44083 55003, 0.61787 62444 02644,
226 4 0.45801 67776 57227, 0.28160 35507 79259, 0.09501 25098 37637/
227 result=0.
228 eps=0.
229 au=a
230 c=(b-a)/float(lambda)
231 DO 1 i = 1,lambda
232 fi=i
233 ao=a+fi*c
234 c1=0.5*(au+ao)
235 c2=c1-au
236 s8=0.
237 s16=0.
238 DO 2 j = 1,4
239 u=x(j)*c2
240 2 s8=s8+w(j)*(f(c1+u)+f(c1-u))
241 DO 3 j = 5,12
242 u=x(j)*c2
243 3 s16=s16+w(j)*(f(c1+u)+f(c1-u))
244 result=result+c2*s16
245 eps=eps+abs(c2*(s16-s8))
246 1 au=ao
247 RETURN
248 END
249
250 SUBROUTINE stepgauss2(F,A,B,LAMBDA,RESULT,EPS)
251 implicit real*8 (a-h,o-z)
252 dimension w(12),x(12)
253 DATA w
254 1/0.10122 85362 90376, 0.22238 10344 53374, 0.31370 66458 77887,
255 2 0.36268 37833 78362, 0.02715 24594 11754, 0.06225 35239 38648,
256 3 0.09515 85116 82493, 0.12462 89712 55534, 0.14959 59888 16577,
257 4 0.16915 65193 95003, 0.18260 34150 44924, 0.18945 06104 55069/
258 DATA x
259 1/0.96028 98564 97536, 0.79666 64774 13627, 0.52553 24099 16329,
260 2 0.18343 46424 95650, 0.98940 09349 91650, 0.94457 50230 73233,
261 3 0.86563 12023 87832, 0.75540 44083 55003, 0.61787 62444 02644,
262 4 0.45801 67776 57227, 0.28160 35507 79259, 0.09501 25098 37637/
263 result=0.
264 eps=0.
265 au=a
266 c=(b-a)/float(lambda)
267 DO 1 i = 1,lambda
268 fi=i
269 ao=a+fi*c
270 c1=0.5*(au+ao)
271 c2=c1-au
272 s8=0.
273 s16=0.
274 DO 2 j = 1,4
275 u=x(j)*c2
276 2 s8=s8+w(j)*(f(c1+u)+f(c1-u))
277 DO 3 j = 5,12
278 u=x(j)*c2
279 3 s16=s16+w(j)*(f(c1+u)+f(c1-u))
280 result=result+c2*s16
281 eps=eps+abs(c2*(s16-s8))
282 1 au=ao
283 RETURN
284 END
285
286 SUBROUTINE stepgauss3(F,A,B,LAMBDA,RESULT,EPS)
287 implicit real*8 (a-h,o-z)
288 dimension w(12),x(12)
289 DATA w
290 1/0.10122 85362 90376, 0.22238 10344 53374, 0.31370 66458 77887,
291 2 0.36268 37833 78362, 0.02715 24594 11754, 0.06225 35239 38648,
292 3 0.09515 85116 82493, 0.12462 89712 55534, 0.14959 59888 16577,
293 4 0.16915 65193 95003, 0.18260 34150 44924, 0.18945 06104 55069/
294 DATA x
295 1/0.96028 98564 97536, 0.79666 64774 13627, 0.52553 24099 16329,
296 2 0.18343 46424 95650, 0.98940 09349 91650, 0.94457 50230 73233,
297 3 0.86563 12023 87832, 0.75540 44083 55003, 0.61787 62444 02644,
298 4 0.45801 67776 57227, 0.28160 35507 79259, 0.09501 25098 37637/
299 result=0.
300 eps=0.
301 au=a
302 c=(b-a)/float(lambda)
303 DO 1 i = 1,lambda
304 fi=i
305 ao=a+fi*c
306 c1=0.5*(au+ao)
307 c2=c1-au
308 s8=0.
309 s16=0.
310 DO 2 j = 1,4
311 u=x(j)*c2
312 2 s8=s8+w(j)*(f(c1+u)+f(c1-u))
313 DO 3 j = 5,12
314 u=x(j)*c2
315 3 s16=s16+w(j)*(f(c1+u)+f(c1-u))
316 result=result+c2*s16
317 eps=eps+abs(c2*(s16-s8))
318 1 au=ao
319 RETURN
320 END
321
322*******************************************************************************
323* FUNCTIONS FOR CHANGE OF VARIABLES *******************************************
324*******************************************************************************
325 DOUBLE PRECISION FUNCTION f_change(X,F,AMS1,AMS2)
326 IMPLICIT NONE
327 EXTERNAL f
328 DOUBLE PRECISION F
329 DOUBLE PRECISION X, NEW_X
330 DOUBLE PRECISION AMS1,AMS2
331 DOUBLE PRECISION ALP1,ALP2,ALP
332 DOUBLE PRECISION AMRX,GAMRX
333 DATA amrx/0.77/,gamrx/1.8/
334
335 alp1 = atan((ams1-amrx**2)/amrx/gamrx) ! integration range for ALP
336 alp2 = atan((ams2-amrx**2)/amrx/gamrx) ! integration range for ALP
337 alp = alp1 + x*(alp2-alp1) ! change of variables
338
339 new_x = amrx**2+amrx*gamrx*tan(alp) ! second change of variables
340
341 f_change = f(new_x)
342
343! Jacobian for change of variables
344 f_change = f_change * (alp2-alp1) * ((new_x-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
345
346 RETURN
347 END
348
349 DOUBLE PRECISION FUNCTION f_change2(X,F,AMS1,AMS2)
350 IMPLICIT NONE
351 EXTERNAL f
352 DOUBLE PRECISION F
353 DOUBLE PRECISION X, NEW_X
354 DOUBLE PRECISION AMS1,AMS2
355 DOUBLE PRECISION ALP1,ALP2,ALP
356 DOUBLE PRECISION AMRX,GAMRX
357 DATA amrx/0.77/,gamrx/1.8/
358
359 alp1 = atan((ams1-amrx**2)/amrx/gamrx) ! integration range for ALP
360 alp2 = atan((ams2-amrx**2)/amrx/gamrx) ! integration range for ALP
361 alp = alp1 + x*(alp2-alp1) ! change of variables
362
363 new_x = amrx**2+amrx*gamrx*tan(alp) ! second change of variables
364
365 f_change2 = f(new_x)
366
367! Jacobian for change of variables
368 f_change2 = f_change2 * (alp2-alp1) * ((new_x-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
369
370 RETURN
371 END
372
373 DOUBLE PRECISION FUNCTION f_change3(X,F,AMS1,AMS2)
374 IMPLICIT NONE
375 EXTERNAL f
376 DOUBLE PRECISION F
377 DOUBLE PRECISION X, NEW_X
378 DOUBLE PRECISION AMS1,AMS2
379 DOUBLE PRECISION ALP1,ALP2,ALP
380 DOUBLE PRECISION AMRX,GAMRX
381 DATA amrx/0.77/,gamrx/1.8/
382
383 alp1 = atan((ams1-amrx**2)/amrx/gamrx) ! integration range for ALP
384 alp2 = atan((ams2-amrx**2)/amrx/gamrx) ! integration range for ALP
385 alp = alp1 + x*(alp2-alp1) ! change of variables
386
387 new_x = amrx**2+amrx*gamrx*tan(alp) ! second change of variables
388
389 f_change3 = f(new_x)
390
391! Jacobian for change of variables
392 f_change3 = f_change3 * (alp2-alp1) * ((new_x-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
393
394 RETURN
395 END
396
397*******************************************************************************
398* INTEGRATION WITH CHANGED VARIABLES ******************************************
399*******************************************************************************
400 SUBROUTINE changegauss(F,AMS1,AMS2,LAMBDA,RESULT,EPS)
401 implicit real*8 (a-h,o-z)
402 EXTERNAL f,f2
403 dimension w(12),x(12)
404 DATA w
405 1/0.10122 85362 90376, 0.22238 10344 53374, 0.31370 66458 77887,
406 2 0.36268 37833 78362, 0.02715 24594 11754, 0.06225 35239 38648,
407 3 0.09515 85116 82493, 0.12462 89712 55534, 0.14959 59888 16577,
408 4 0.16915 65193 95003, 0.18260 34150 44924, 0.18945 06104 55069/
409 DATA x
410 1/0.96028 98564 97536, 0.79666 64774 13627, 0.52553 24099 16329,
411 2 0.18343 46424 95650, 0.98940 09349 91650, 0.94457 50230 73233,
412 3 0.86563 12023 87832, 0.75540 44083 55003, 0.61787 62444 02644,
413 4 0.45801 67776 57227, 0.28160 35507 79259, 0.09501 25098 37637/
414 result=0.
415 eps=0.
416
417! Range is changed from [AMS1,AMS2] to [0,1]
418! F, AMS1 and AMS2 are parameters of F_CHANGE
419 a=0
420 b=1
421
422 au=a
423 c=(b-a)/float(lambda)
424 DO 1 i = 1,lambda
425 fi=i
426 ao=a+fi*c
427 c1=0.5*(au+ao)
428 c2=c1-au
429 s8=0.
430 s16=0.
431 DO 2 j = 1,4
432 u=x(j)*c2
433 2 s8=s8+w(j)*(f_change(c1+u,f,ams1,ams2)+f_change(c1-u,f,ams1,ams2))
434 DO 3 j = 5,12
435 u=x(j)*c2
436 3 s16=s16+w(j)*(f_change(c1+u,f,ams1,ams2)+f_change(c1-u,f,ams1,ams2))
437 result=result+c2*s16
438 eps=eps+abs(c2*(s16-s8))
439 1 au=ao
440 RETURN
441 END
442
443 SUBROUTINE changegauss2(F,AMS1,AMS2,LAMBDA,RESULT,EPS)
444 implicit real*8 (a-h,o-z)
445 EXTERNAL f,f2
446 dimension w(12),x(12)
447 DATA w
448 1/0.10122 85362 90376, 0.22238 10344 53374, 0.31370 66458 77887,
449 2 0.36268 37833 78362, 0.02715 24594 11754, 0.06225 35239 38648,
450 3 0.09515 85116 82493, 0.12462 89712 55534, 0.14959 59888 16577,
451 4 0.16915 65193 95003, 0.18260 34150 44924, 0.18945 06104 55069/
452 DATA x
453 1/0.96028 98564 97536, 0.79666 64774 13627, 0.52553 24099 16329,
454 2 0.18343 46424 95650, 0.98940 09349 91650, 0.94457 50230 73233,
455 3 0.86563 12023 87832, 0.75540 44083 55003, 0.61787 62444 02644,
456 4 0.45801 67776 57227, 0.28160 35507 79259, 0.09501 25098 37637/
457 result=0.
458 eps=0.
459
460! Range is changed from [AMS1,AMS2] to [0,1]
461! F, AMS1 and AMS2 are parameters of F_CHANGE
462 a=0
463 b=1
464
465 au=a
466 c=(b-a)/float(lambda)
467 DO 1 i = 1,lambda
468 fi=i
469 ao=a+fi*c
470 c1=0.5*(au+ao)
471 c2=c1-au
472 s8=0.
473 s16=0.
474 DO 2 j = 1,4
475 u=x(j)*c2
476 2 s8=s8+w(j)*(f_change2(c1+u,f,ams1,ams2)+f_change2(c1-u,f,ams1,ams2))
477 DO 3 j = 5,12
478 u=x(j)*c2
479 3 s16=s16+w(j)*(f_change2(c1+u,f,ams1,ams2)+f_change2(c1-u,f,ams1,ams2))
480 result=result+c2*s16
481 eps=eps+abs(c2*(s16-s8))
482 1 au=ao
483 RETURN
484 END
485
486 SUBROUTINE changegauss3(F,AMS1,AMS2,LAMBDA,RESULT,EPS)
487 implicit real*8 (a-h,o-z)
488 EXTERNAL f,f2
489 dimension w(12),x(12)
490 DATA w
491 1/0.10122 85362 90376, 0.22238 10344 53374, 0.31370 66458 77887,
492 2 0.36268 37833 78362, 0.02715 24594 11754, 0.06225 35239 38648,
493 3 0.09515 85116 82493, 0.12462 89712 55534, 0.14959 59888 16577,
494 4 0.16915 65193 95003, 0.18260 34150 44924, 0.18945 06104 55069/
495 DATA x
496 1/0.96028 98564 97536, 0.79666 64774 13627, 0.52553 24099 16329,
497 2 0.18343 46424 95650, 0.98940 09349 91650, 0.94457 50230 73233,
498 3 0.86563 12023 87832, 0.75540 44083 55003, 0.61787 62444 02644,
499 4 0.45801 67776 57227, 0.28160 35507 79259, 0.09501 25098 37637/
500 result=0.
501 eps=0.
502
503! Range is changed from [AMS1,AMS2] to [0,1]
504! F, AMS1 and AMS2 are parameters of F_CHANGE
505 a=0
506 b=1
507
508 au=a
509 c=(b-a)/float(lambda)
510 DO 1 i = 1,lambda
511 fi=i
512 ao=a+fi*c
513 c1=0.5*(au+ao)
514 c2=c1-au
515 s8=0.
516 s16=0.
517 DO 2 j = 1,4
518 u=x(j)*c2
519 2 s8=s8+w(j)*(f_change3(c1+u,f,ams1,ams2)+f_change3(c1-u,f,ams1,ams2))
520 DO 3 j = 5,12
521 u=x(j)*c2
522 3 s16=s16+w(j)*(f_change3(c1+u,f,ams1,ams2)+f_change3(c1-u,f,ams1,ams2))
523 result=result+c2*s16
524 eps=eps+abs(c2*(s16-s8))
525 1 au=ao
526 RETURN
527 END