1      SUBROUTINE bratio(a,b,x,y,w,w1,ierr)
2C-----------------------------------------------------------------------
3C
4C            EVALUATION OF THE INCOMPLETE BETA FUNCTION IX(A,B)
5C
6C                     --------------------
7C
8C     IT IS ASSUMED THAT A AND B ARE NONNEGATIVE, AND THAT X .LE. 1
9C     AND Y = 1 - X.  BRATIO ASSIGNS W AND W1 THE VALUES
10C
11C                      W  = IX(A,B)
12C                      W1 = 1 - IX(A,B)
13C
14C     IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS.
15C     IF NO INPUT ERRORS ARE DETECTED THEN IERR IS SET TO 0 AND
16C     W AND W1 ARE COMPUTED. OTHERWISE, IF AN ERROR IS DETECTED,
17C     THEN W AND W1 ARE ASSIGNED THE VALUE 0 AND IERR IS SET TO
18C     ONE OF THE FOLLOWING VALUES ...
19C
20C        IERR = 1  IF A OR B IS NEGATIVE
21C        IERR = 2  IF A = B = 0
22C        IERR = 3  IF X .LT. 0 OR X .GT. 1
23C        IERR = 4  IF Y .LT. 0 OR Y .GT. 1
24C        IERR = 5  IF X + Y .NE. 1
25C        IERR = 6  IF X = A = 0
26C        IERR = 7  IF Y = B = 0
27C
28C--------------------
29C     WRITTEN BY ALFRED H. MORRIS, JR.
30C        NAVAL SURFACE WARFARE CENTER
31C        DAHLGREN, VIRGINIA
32C     REVISED ... NOV 1991
33C-----------------------------------------------------------------------
34C     .. Scalar Arguments ..
35      DOUBLE PRECISION a,b,w,w1,x,y
36      INTEGER ierr
37C     ..
38C     .. Local Scalars ..
39      DOUBLE PRECISION a0,b0,eps,lambda,t,x0,y0,z
40      INTEGER ierr1,ind,n
41C     ..
42C     .. External Functions ..
43      DOUBLE PRECISION apser,basym,bfrac,bpser,bup,fpser,spmpar
44      EXTERNAL apser,basym,bfrac,bpser,bup,fpser,spmpar
45C     ..
46C     .. External Subroutines ..
47      EXTERNAL bgrat
48C     ..
49C     .. Intrinsic Functions ..
50      INTRINSIC abs,dmax1,dmin1
51C     ..
52C     .. Executable Statements ..
53C-----------------------------------------------------------------------
54C
55C     ****** EPS IS A MACHINE DEPENDENT CONSTANT. EPS IS THE SMALLEST
56C            FLOATING POINT NUMBER FOR WHICH 1.0 + EPS .GT. 1.0
57C
58      eps = spmpar(1)
59C
60C-----------------------------------------------------------------------
61      w = 0.0D0
62      w1 = 0.0D0
63      IF (a.LT.0.0D0 .OR. b.LT.0.0D0) GO TO 270
64      IF (a.EQ.0.0D0 .AND. b.EQ.0.0D0) GO TO 280
65      IF (x.LT.0.0D0 .OR. x.GT.1.0D0) GO TO 290
66      IF (y.LT.0.0D0 .OR. y.GT.1.0D0) GO TO 300
67      z = ((x+y)-0.5D0) - 0.5D0
68      IF (abs(z).GT.3.0D0*eps) GO TO 310
69C
70      ierr = 0
71      IF (x.EQ.0.0D0) GO TO 210
72      IF (y.EQ.0.0D0) GO TO 230
73      IF (a.EQ.0.0D0) GO TO 240
74      IF (b.EQ.0.0D0) GO TO 220
75C
76      eps = dmax1(eps,1.D-15)
77      IF (dmax1(a,b).LT.1.D-3*eps) GO TO 260
78C
79      ind = 0
80      a0 = a
81      b0 = b
82      x0 = x
83      y0 = y
84      IF (dmin1(a0,b0).GT.1.0D0) GO TO 40
85C
86C             PROCEDURE FOR A0 .LE. 1 OR B0 .LE. 1
87C
88      IF (x.LE.0.5D0) GO TO 10
89      ind = 1
90      a0 = b
91      b0 = a
92      x0 = y
93      y0 = x
94C
95   10 IF (b0.LT.dmin1(eps,eps*a0)) GO TO 90
96      IF (a0.LT.dmin1(eps,eps*b0) .AND. b0*x0.LE.1.0D0) GO TO 100
97      IF (dmax1(a0,b0).GT.1.0D0) GO TO 20
98      IF (a0.GE.dmin1(0.2D0,b0)) GO TO 110
99      IF (x0**a0.LE.0.9D0) GO TO 110
100      IF (x0.GE.0.3D0) GO TO 120
101      n = 20
102      GO TO 140
103C
104   20 IF (b0.LE.1.0D0) GO TO 110
105      IF (x0.GE.0.3D0) GO TO 120
106      IF (x0.GE.0.1D0) GO TO 30
107      IF ((x0*b0)**a0.LE.0.7D0) GO TO 110
108   30 IF (b0.GT.15.0D0) GO TO 150
109      n = 20
110      GO TO 140
111C
112C             PROCEDURE FOR A0 .GT. 1 AND B0 .GT. 1
113C
114   40 IF (a.GT.b) GO TO 50
115      lambda = a - (a+b)*x
116      GO TO 60
117
118   50 lambda = (a+b)*y - b
119   60 IF (lambda.GE.0.0D0) GO TO 70
120      ind = 1
121      a0 = b
122      b0 = a
123      x0 = y
124      y0 = x
125      lambda = abs(lambda)
126C
127   70 IF (b0.LT.40.0D0 .AND. b0*x0.LE.0.7D0) GO TO 110
128      IF (b0.LT.40.0D0) GO TO 160
129      IF (a0.GT.b0) GO TO 80
130      IF (a0.LE.100.0D0) GO TO 130
131      IF (lambda.GT.0.03D0*a0) GO TO 130
132      GO TO 200
133
134   80 IF (b0.LE.100.0D0) GO TO 130
135      IF (lambda.GT.0.03D0*b0) GO TO 130
136      GO TO 200
137C
138C            EVALUATION OF THE APPROPRIATE ALGORITHM
139C
140   90 w = fpser(a0,b0,x0,eps)
141      w1 = 0.5D0 + (0.5D0-w)
142      GO TO 250
143C
144  100 w1 = apser(a0,b0,x0,eps)
145      w = 0.5D0 + (0.5D0-w1)
146      GO TO 250
147C
148  110 w = bpser(a0,b0,x0,eps)
149      w1 = 0.5D0 + (0.5D0-w)
150      GO TO 250
151C
152  120 w1 = bpser(b0,a0,y0,eps)
153      w = 0.5D0 + (0.5D0-w1)
154      GO TO 250
155C
156  130 w = bfrac(a0,b0,x0,y0,lambda,15.0D0*eps)
157      w1 = 0.5D0 + (0.5D0-w)
158      GO TO 250
159C
160  140 w1 = bup(b0,a0,y0,x0,n,eps)
161      b0 = b0 + n
162  150 CALL bgrat(b0,a0,y0,x0,w1,15.0D0*eps,ierr1)
163      w = 0.5D0 + (0.5D0-w1)
164      GO TO 250
165C
166  160 n = b0
167      b0 = b0 - n
168      IF (b0.NE.0.0D0) GO TO 170
169      n = n - 1
170      b0 = 1.0D0
171  170 w = bup(b0,a0,y0,x0,n,eps)
172      IF (x0.GT.0.7D0) GO TO 180
173      w = w + bpser(a0,b0,x0,eps)
174      w1 = 0.5D0 + (0.5D0-w)
175      GO TO 250
176C
177  180 IF (a0.GT.15.0D0) GO TO 190
178      n = 20
179      w = w + bup(a0,b0,x0,y0,n,eps)
180      a0 = a0 + n
181  190 CALL bgrat(a0,b0,x0,y0,w,15.0D0*eps,ierr1)
182      w1 = 0.5D0 + (0.5D0-w)
183      GO TO 250
184C
185  200 w = basym(a0,b0,lambda,100.0D0*eps)
186      w1 = 0.5D0 + (0.5D0-w)
187      GO TO 250
188C
189C               TERMINATION OF THE PROCEDURE
190C
191  210 IF (a.EQ.0.0D0) GO TO 320
192  220 w = 0.0D0
193      w1 = 1.0D0
194      RETURN
195C
196  230 IF (b.EQ.0.0D0) GO TO 330
197  240 w = 1.0D0
198      w1 = 0.0D0
199      RETURN
200C
201  250 IF (ind.EQ.0) RETURN
202      t = w
203      w = w1
204      w1 = t
205      RETURN
206C
207C           PROCEDURE FOR A AND B .LT. 1.E-3*EPS
208C
209  260 w = b/ (a+b)
210      w1 = a/ (a+b)
211      RETURN
212C
213C                       ERROR RETURN
214C
215  270 ierr = 1
216      RETURN
217
218  280 ierr = 2
219      RETURN
220
221  290 ierr = 3
222      RETURN
223
224  300 ierr = 4
225      RETURN
226
227  310 ierr = 5
228      RETURN
229
230  320 ierr = 6
231      RETURN
232
233  330 ierr = 7
234      RETURN
235
236      END
237