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 
24 C-----------------------------------------------------------------------------
25 C NEW INTEGRATION METHOD - constant step size (provided by fourth parameter)
26 C comment out following 5 lines to return to old solutions of advantages too
27 C-----------------------------------------------------------------------------
28  result=0
29  eps=0
30 c CALL STEPGAUSS(F,A,B,6,RESULT,EPS)
31  CALL changegauss(f,a,b,2,result,eps)
32  gaus=result
33  RETURN
34 C-----------------------------------------------------------------------------
35 C PREVIOUS INTEGRATION METHOD - step size depends on EPSSQ
36 C-----------------------------------------------------------------------------
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 
90 C-----------------------------------------------------------------------------
91 C NEW INTEGRATION METHOD - constant step size (provided by fourth parameter)
92 C comment out following 5 lines to return to old solutions of advantages too
93 C-----------------------------------------------------------------------------
94  result=0
95  eps=0
96 c CALL STEPGAUSS2(F,A,B,6,RESULT,EPS)
97  CALL changegauss2(f,a,b,2,result,eps)
98  gaus2=result
99  RETURN
100 C-----------------------------------------------------------------------------
101 C PREVIOUS INTEGRATION METHOD - step size depends on EPSSQ
102 C-----------------------------------------------------------------------------
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 
155 C-----------------------------------------------------------------------------
156 C NEW INTEGRATION METHOD - constant step size (provided by fourth parameter)
157 C comment out following 5 lines to return to old solutions of advantages too
158 C-----------------------------------------------------------------------------
159  result=0
160  eps=0
161 c CALL STEPGAUSS3(F,A,B,6,RESULT,EPS)
162  CALL changegauss3(f,a,b,2,result,eps)
163  gaus3=result
164  RETURN
165 C-----------------------------------------------------------------------------
166 C PREVIOUS INTEGRATION METHOD - step size depends on EPSSQ
167 C-----------------------------------------------------------------------------
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 
203 C--------------------------------------------------------------------------------------------------
204 C--------------------------------------------------------------------------------------------------
205 C--------------------------------------------------------------------------------------------------
206 C OBSOLETE 1.13 901108 3.40 CERN PROGRAM LIBRARY OBSOLETE PAM
207 C CERN OBSOLETE PROGRAMS AND SUBROUTINES PAM-FILE.
208 C MANY OF THE ROUTINES ARE FOR ONE FORTRAN DIALECT ONLY.
209 C THE FLAG PATCH F77 SHOULD BE USED FOR THE FORTRAN 77 VERSIONS
210 C SHOULD THEY EXIST - IF NOT FORTRAN 4 WILL BE OBTAINED.
211 C THE CERN LIBRARIES GENLIB AND KERNLIB ARE USUALLY REQUIRED.
212 C OBSOLETE ROUTINES OF THE PROGRAM LIBRARY ARE OFTEN REQUIRED.
213 C
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
Generated by doxygen 1.9.1