EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
gadap.F
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file gadap.F
1 
2 C#######################################################################
3 C
4 C One- and two-dimensional adaptive Gaussian integration routines.
5 C
6 C **********************************************************************
7 
8  SUBROUTINE gadap(A0,B0,F,EPS,SUM)
9 
10  IMPLICIT NONE
11 C
12 C PURPOSE - INTEGRATE A FUNCTION F(X)
13 C METHOD - ADAPTIVE GAUSSIAN
14 C USAGE - CALL GADAP(A0,B0,F,EPS,SUM)
15 C PARAMETERS A0 - LOWER LIMIT (INPUT,REAL)
16 C B0 - UPPER LIMIT (INPUT,REAL)
17 C F - FUNCTION F(X) TO BE INTEGRATED. MUST BE
18 C SUPPLIED BY THE USER. (INPUT,REAL FUNCTION)
19 C EPS - DESIRED RELATIVE ACCURACY. IF SUM IS SMALL EPS
20 C WILL BE ABSOLUTE ACCURACY INSTEAD. (INPUT,REAL)
21 C SUM - CALCULATED VALUE FOR THE INTEGRAL (OUTPUT,REAL)
22 C PRECISION - SINGLE
23 C REQ'D PROG'S - F
24 C AUTHOR - T. JOHANSSON, LUND UNIV. COMPUTER CENTER, 1973
25 C REFERENCE(S) - THE AUSTRALIAN COMPUTER JOURNAL,3 P.126 AUG. -71
26 C
27  common/gadap1/ num,ifu
28  INTEGER num,ifu
29  SAVE /gadap1/
30 
31  INTEGER n,i,l
32  REAL aa,bb,f1f,f2f,f3f
33  REAL a0,b0,f,eps,sum,a,b,f1,f2,f3,s
34  REAL dsum,c,red,w1,u2,ss,sold
35  EXTERNAL f
36  dimension a(300),b(300),f1(300),f2(300),f3(300),s(300),n(300)
37  1 FORMAT(16h gadap:i too big)
38  dsum(f1f,f2f,f3f,aa,bb)=5./18.*(bb-aa)*(f1f+1.6*f2f+f3f)
39  IF(eps.LT.1.0e-8) eps=1.0e-8
40  red=1.3
41  l=1
42  i=1
43  sum=0.
44  c=sqrt(15.)/5.
45  a(1)=a0
46  b(1)=b0
47  f1(1)=f(0.5*(1+c)*a0+0.5*(1-c)*b0)
48  f2(1)=f(0.5*(a0+b0))
49  f3(1)=f(0.5*(1-c)*a0+0.5*(1+c)*b0)
50  ifu=3
51  s(1)= dsum(f1(1),f2(1),f3(1),a0,b0)
52  100 CONTINUE
53  l=l+1
54  n(l)=3
55  eps=eps*red
56  a(i+1)=a(i)+c*(b(i)-a(i))
57  b(i+1)=b(i)
58  a(i+2)=a(i)+b(i)-a(i+1)
59  b(i+2)=a(i+1)
60  a(i+3)=a(i)
61  b(i+3)=a(i+2)
62  w1=a(i)+(b(i)-a(i))/5.
63  u2=2.*w1-(a(i)+a(i+2))/2.
64  f1(i+1)=f(a(i)+b(i)-w1)
65  f2(i+1)=f3(i)
66  f3(i+1)=f(b(i)-a(i+2)+w1)
67  f1(i+2)=f(u2)
68  f2(i+2)=f2(i)
69  f3(i+2)=f(b(i+2)+a(i+2)-u2)
70  f1(i+3)=f(a(i)+a(i+2)-w1)
71  f2(i+3)=f1(i)
72  f3(i+3)=f(w1)
73  ifu=ifu+6
74  IF(ifu.GT.5000) goto 130
75  s(i+1)= dsum(f1(i+1),f2(i+1),f3(i+1),a(i+1),b(i+1))
76  s(i+2)= dsum(f1(i+2),f2(i+2),f3(i+2),a(i+2),b(i+2))
77  s(i+3)= dsum(f1(i+3),f2(i+3),f3(i+3),a(i+3),b(i+3))
78  ss=s(i+1)+s(i+2)+s(i+3)
79  i=i+3
80  IF(i.GT.300)goto 120
81  sold=s(i-3)
82  IF(abs(sold-ss).GT.eps*(1.+abs(ss))/2.) goto 100
83  sum=sum+ss
84  i=i-4
85  n(l)=0
86  l=l-1
87  110 CONTINUE
88  IF(l.EQ.1) goto 130
89  n(l)=n(l)-1
90  eps=eps/red
91  IF(n(l).NE.0) goto 100
92  i=i-1
93  l=l-1
94  goto 110
95  120 WRITE(6,1)
96  130 RETURN
97  END