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