EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
gadap2.F
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file gadap2.F
1 
2 C **********************************************************************
3 
4  SUBROUTINE gadap2(A0,B0,FL,FU,F,EPS,SUM)
5 
6  IMPLICIT NONE
7 C
8 C PURPOSE - INTEGRATE A FUNCTION F(X,Y) OF TWO VARIABLES
9 C METHOD - ADAPTIVE GAUSSIAN IN BOTH DIRECTIONS
10 C USAGE - CALL GADAP2(A0,B0,FL,FU,F,EPS,SUM)
11 C PARAMETERS A0 - LOWER X-LIMIT (INPUT,REAL)
12 C B0 - UPPER X-LIMIT (INPUT,REAL)
13 C FL - USER SUPPLIED FUNCTION FL(X) GIVING THE LOWER
14 C Y-LIMIT FOR A GIVEN X-VALUE
15 C (INPUT,REAL FUNCTION)
16 C FU - USER SUPPLIED FUNCTION FU(X) GIVING THE UPPER
17 C Y-LIMIT FOR A GIVEN X-VALUE
18 C (INPUT,REAL FUNCTION)
19 C F - USER SUPPLIED FUNCTION F(X,Y) TO BE INTEGRATED
20 C (INPUT,REAL FUNCTION)
21 C EPS - DESIRED ACCURACY (INPUT,REAL)
22 C SUM - CALCULATED VALUE FOR THE INTEGRAL (OUTPUT,REAL)
23 C PRECISION - SINGLE
24 C REQ'D PROG'S - FL,FU,F,GADAPF
25 C AUTHOR - THOMAS JOHANSSON, LDC,1973
26 C
27  common/gadap1/ num,ifu
28  INTEGER num,ifu
29  SAVE /gadap1/
30 
31  INTEGER n,l,i
32  REAL aa,bb,f1f,f2f,f3f
33  REAL a0,b0,fl,fu,f,eps,sum,a,b,f1,f2,f3,s
34  REAL dsum,red,c,x,ay,by,w1,u2,ss,sold
35  REAL gadapf
36  EXTERNAL f,fl,fu
37  dimension a(300),b(300),f1(300),f2(300),f3(300),s(300),n(300)
38  1 FORMAT(16h gadap:i too big)
39  dsum(f1f,f2f,f3f,aa,bb)=5./18.*(bb-aa)*(f1f+1.6*f2f+f3f)
40  IF(eps.LT.1.0e-8) eps=1.0e-8
41  red=1.4
42  l=1
43  i=1
44  sum=0.
45  c=sqrt(15.)/5.
46  a(1)=a0
47  b(1)=b0
48  x=0.5*(1+c)*a0+0.5*(1-c)*b0
49  ay=fl(x)
50  by=fu(x)
51  f1(1)=gadapf(x,ay,by,f,eps)
52  x=0.5*(a0+b0)
53  ay=fl(x)
54  by=fu(x)
55  f2(1)=gadapf(x,ay,by,f,eps)
56  x=0.5*(1-c)*a0+0.5*(1+c)*b0
57  ay=fl(x)
58  by=fu(x)
59  f3(1)=gadapf(x,ay,by,f,eps)
60  ifu=3
61  s(1)= dsum(f1(1),f2(1),f3(1),a0,b0)
62  100 CONTINUE
63  l=l+1
64  n(l)=3
65  eps=eps*red
66  a(i+1)=a(i)+c*(b(i)-a(i))
67  b(i+1)=b(i)
68  a(i+2)=a(i)+b(i)-a(i+1)
69  b(i+2)=a(i+1)
70  a(i+3)=a(i)
71  b(i+3)=a(i+2)
72  w1=a(i)+(b(i)-a(i))/5.
73  u2=2.*w1-(a(i)+a(i+2))/2.
74  x=a(i)+b(i)-w1
75  ay=fl(x)
76  by=fu(x)
77  f1(i+1)=gadapf(x,ay,by,f,eps)
78  f2(i+1)=f3(i)
79  x=b(i)-a(i+2)+w1
80  ay=fl(x)
81  by=fu(x)
82  f3(i+1)=gadapf(x,ay,by,f,eps)
83  x=u2
84  ay=fl(x)
85  by=fu(x)
86  f1(i+2)=gadapf(x,ay,by,f,eps)
87  f2(i+2)=f2(i)
88  x=b(i+2)+a(i+2)-u2
89  ay=fl(x)
90  by=fu(x)
91  f3(i+2)=gadapf(x,ay,by,f,eps)
92  x=a(i)+a(i+2)-w1
93  ay=fl(x)
94  by=fu(x)
95  f1(i+3)=gadapf(x,ay,by,f,eps)
96  f2(i+3)=f1(i)
97  x=w1
98  ay=fl(x)
99  by=fu(x)
100  f3(i+3)=gadapf(x,ay,by,f,eps)
101  ifu=ifu+6
102  IF(ifu.GT.5000) goto 130
103  s(i+1)= dsum(f1(i+1),f2(i+1),f3(i+1),a(i+1),b(i+1))
104  s(i+2)= dsum(f1(i+2),f2(i+2),f3(i+2),a(i+2),b(i+2))
105  s(i+3)= dsum(f1(i+3),f2(i+3),f3(i+3),a(i+3),b(i+3))
106  ss=s(i+1)+s(i+2)+s(i+3)
107  i=i+3
108  IF(i.GT.300)goto 120
109  sold=s(i-3)
110  IF(abs(sold-ss).GT.eps*(1.+abs(ss))/2.) goto 100
111  sum=sum+ss
112  i=i-4
113  n(l)=0
114  l=l-1
115  110 CONTINUE
116  IF(l.EQ.1) goto 130
117  n(l)=n(l)-1
118  eps=eps/red
119  IF(n(l).NE.0) goto 100
120  i=i-1
121  l=l-1
122  goto 110
123  120 WRITE(6,1)
124  130 RETURN
125  END