00001 #include "cdflib.h"
00002 void bratio(double *a,double *b,double *x,double *y,double *w,
00003             double *w1,int *ierr)
00004 /*
00005 -----------------------------------------------------------------------
00006  
00007             EVALUATION OF THE INCOMPLETE BETA FUNCTION IX(A,B)
00008  
00009                      --------------------
00010  
00011      IT IS ASSUMED THAT A AND B ARE NONNEGATIVE, AND THAT X .LE. 1
00012      AND Y = 1 - X.  BRATIO ASSIGNS W AND W1 THE VALUES
00013  
00014                       W  = IX(A,B)
00015                       W1 = 1 - IX(A,B)
00016  
00017      IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS.
00018      IF NO INPUT ERRORS ARE DETECTED THEN IERR IS SET TO 0 AND
00019      W AND W1 ARE COMPUTED. OTHERWISE, IF AN ERROR IS DETECTED,
00020      THEN W AND W1 ARE ASSIGNED THE VALUE 0 AND IERR IS SET TO
00021      ONE OF THE FOLLOWING VALUES ...
00022  
00023         IERR = 1  IF A OR B IS NEGATIVE
00024         IERR = 2  IF A = B = 0
00025         IERR = 3  IF X .LT. 0 OR X .GT. 1
00026         IERR = 4  IF Y .LT. 0 OR Y .GT. 1
00027         IERR = 5  IF X + Y .NE. 1
00028         IERR = 6  IF X = A = 0
00029         IERR = 7  IF Y = B = 0
00030  
00031 --------------------
00032      WRITTEN BY ALFRED H. MORRIS, JR.
00033         NAVAL SURFACE WARFARE CENTER
00034         DAHLGREN, VIRGINIA
00035      REVISED ... NOV 1991
00036 -----------------------------------------------------------------------
00037 */
00038 {
00039 static int K1 = 1;
00040 static double a0,b0,eps,lambda,t,x0,y0,z;
00041 static int ierr1,ind,n;
00042 static double T2,T3,T4,T5;
00043 /*
00044      ..
00045      .. Executable Statements ..
00046 */
00047 /*
00048      ****** EPS IS A MACHINE DEPENDENT CONSTANT. EPS IS THE SMALLEST
00049             FLOATING POINT NUMBER FOR WHICH 1.0 + EPS .GT. 1.0
00050 */
00051     eps = spmpar(&K1);
00052     *w = *w1 = 0.0e0;
00053     if(*a < 0.0e0 || *b < 0.0e0) goto S270;
00054     if(*a == 0.0e0 && *b == 0.0e0) goto S280;
00055     if(*x < 0.0e0 || *x > 1.0e0) goto S290;
00056     if(*y < 0.0e0 || *y > 1.0e0) goto S300;
00057     z = *x+*y-0.5e0-0.5e0;
00058     if(fabs(z) > 3.0e0*eps) goto S310;
00059     *ierr = 0;
00060     if(*x == 0.0e0) goto S210;
00061     if(*y == 0.0e0) goto S230;
00062     if(*a == 0.0e0) goto S240;
00063     if(*b == 0.0e0) goto S220;
00064     eps = fifdmax1(eps,1.e-15);
00065     if(fifdmax1(*a,*b) < 1.e-3*eps) goto S260;
00066     ind = 0;
00067     a0 = *a;
00068     b0 = *b;
00069     x0 = *x;
00070     y0 = *y;
00071     if(fifdmin1(a0,b0) > 1.0e0) goto S40;
00072 /*
00073              PROCEDURE FOR A0 .LE. 1 OR B0 .LE. 1
00074 */
00075     if(*x <= 0.5e0) goto S10;
00076     ind = 1;
00077     a0 = *b;
00078     b0 = *a;
00079     x0 = *y;
00080     y0 = *x;
00081 S10:
00082     if(b0 < fifdmin1(eps,eps*a0)) goto S90;
00083     if(a0 < fifdmin1(eps,eps*b0) && b0*x0 <= 1.0e0) goto S100;
00084     if(fifdmax1(a0,b0) > 1.0e0) goto S20;
00085     if(a0 >= fifdmin1(0.2e0,b0)) goto S110;
00086     if(pow(x0,a0) <= 0.9e0) goto S110;
00087     if(x0 >= 0.3e0) goto S120;
00088     n = 20;
00089     goto S140;
00090 S20:
00091     if(b0 <= 1.0e0) goto S110;
00092     if(x0 >= 0.3e0) goto S120;
00093     if(x0 >= 0.1e0) goto S30;
00094     if(pow(x0*b0,a0) <= 0.7e0) goto S110;
00095 S30:
00096     if(b0 > 15.0e0) goto S150;
00097     n = 20;
00098     goto S140;
00099 S40:
00100 /*
00101              PROCEDURE FOR A0 .GT. 1 AND B0 .GT. 1
00102 */
00103     if(*a > *b) goto S50;
00104     lambda = *a-(*a+*b)**x;
00105     goto S60;
00106 S50:
00107     lambda = (*a+*b)**y-*b;
00108 S60:
00109     if(lambda >= 0.0e0) goto S70;
00110     ind = 1;
00111     a0 = *b;
00112     b0 = *a;
00113     x0 = *y;
00114     y0 = *x;
00115     lambda = fabs(lambda);
00116 S70:
00117     if(b0 < 40.0e0 && b0*x0 <= 0.7e0) goto S110;
00118     if(b0 < 40.0e0) goto S160;
00119     if(a0 > b0) goto S80;
00120     if(a0 <= 100.0e0) goto S130;
00121     if(lambda > 0.03e0*a0) goto S130;
00122     goto S200;
00123 S80:
00124     if(b0 <= 100.0e0) goto S130;
00125     if(lambda > 0.03e0*b0) goto S130;
00126     goto S200;
00127 S90:
00128 /*
00129             EVALUATION OF THE APPROPRIATE ALGORITHM
00130 */
00131     *w = fpser(&a0,&b0,&x0,&eps);
00132     *w1 = 0.5e0+(0.5e0-*w);
00133     goto S250;
00134 S100:
00135     *w1 = apser(&a0,&b0,&x0,&eps);
00136     *w = 0.5e0+(0.5e0-*w1);
00137     goto S250;
00138 S110:
00139     *w = bpser(&a0,&b0,&x0,&eps);
00140     *w1 = 0.5e0+(0.5e0-*w);
00141     goto S250;
00142 S120:
00143     *w1 = bpser(&b0,&a0,&y0,&eps);
00144     *w = 0.5e0+(0.5e0-*w1);
00145     goto S250;
00146 S130:
00147     T2 = 15.0e0*eps;
00148     *w = bfrac(&a0,&b0,&x0,&y0,&lambda,&T2);
00149     *w1 = 0.5e0+(0.5e0-*w);
00150     goto S250;
00151 S140:
00152     *w1 = bup(&b0,&a0,&y0,&x0,&n,&eps);
00153     b0 += (double)n;
00154 S150:
00155     T3 = 15.0e0*eps;
00156     bgrat(&b0,&a0,&y0,&x0,w1,&T3,&ierr1);
00157     *w = 0.5e0+(0.5e0-*w1);
00158     goto S250;
00159 S160:
00160     n = b0;
00161     b0 -= (double)n;
00162     if(b0 != 0.0e0) goto S170;
00163     n -= 1;
00164     b0 = 1.0e0;
00165 S170:
00166     *w = bup(&b0,&a0,&y0,&x0,&n,&eps);
00167     if(x0 > 0.7e0) goto S180;
00168     *w += bpser(&a0,&b0,&x0,&eps);
00169     *w1 = 0.5e0+(0.5e0-*w);
00170     goto S250;
00171 S180:
00172     if(a0 > 15.0e0) goto S190;
00173     n = 20;
00174     *w += bup(&a0,&b0,&x0,&y0,&n,&eps);
00175     a0 += (double)n;
00176 S190:
00177     T4 = 15.0e0*eps;
00178     bgrat(&a0,&b0,&x0,&y0,w,&T4,&ierr1);
00179     *w1 = 0.5e0+(0.5e0-*w);
00180     goto S250;
00181 S200:
00182     T5 = 100.0e0*eps;
00183     *w = basym(&a0,&b0,&lambda,&T5);
00184     *w1 = 0.5e0+(0.5e0-*w);
00185     goto S250;
00186 S210:
00187 /*
00188                TERMINATION OF THE PROCEDURE
00189 */
00190     if(*a == 0.0e0) goto S320;
00191 S220:
00192     *w = 0.0e0;
00193     *w1 = 1.0e0;
00194     return;
00195 S230:
00196     if(*b == 0.0e0) goto S330;
00197 S240:
00198     *w = 1.0e0;
00199     *w1 = 0.0e0;
00200     return;
00201 S250:
00202     if(ind == 0) return;
00203     t = *w;
00204     *w = *w1;
00205     *w1 = t;
00206     return;
00207 S260:
00208 /*
00209            PROCEDURE FOR A AND B .LT. 1.E-3*EPS
00210 */
00211     *w = *b/(*a+*b);
00212     *w1 = *a/(*a+*b);
00213     return;
00214 S270:
00215 /*
00216                        ERROR RETURN
00217 */
00218     *ierr = 1;
00219     return;
00220 S280:
00221     *ierr = 2;
00222     return;
00223 S290:
00224     *ierr = 3;
00225     return;
00226 S300:
00227     *ierr = 4;
00228     return;
00229 S310:
00230     *ierr = 5;
00231     return;
00232 S320:
00233     *ierr = 6;
00234     return;
00235 S330:
00236     *ierr = 7;
00237     return;
00238 } /* END */