00001 #include "cdflib.h"
00002 double basym(double *a,double *b,double *lambda,double *eps)
00003 /*
00004 -----------------------------------------------------------------------
00005      ASYMPTOTIC EXPANSION FOR IX(A,B) FOR LARGE A AND B.
00006      LAMBDA = (A + B)*Y - B  AND EPS IS THE TOLERANCE USED.
00007      IT IS ASSUMED THAT LAMBDA IS NONNEGATIVE AND THAT
00008      A AND B ARE GREATER THAN OR EQUAL TO 15.
00009 -----------------------------------------------------------------------
00010 */
00011 {
00012 static double e0 = 1.12837916709551e0;
00013 static double e1 = .353553390593274e0;
00014 static int num = 20;
00015 /*
00016 ------------------------
00017      ****** NUM IS THE MAXIMUM VALUE THAT N CAN TAKE IN THE DO LOOP
00018             ENDING AT STATEMENT 50. IT IS REQUIRED THAT NUM BE EVEN.
00019             THE ARRAYS A0, B0, C, D HAVE DIMENSION NUM + 1.
00020 ------------------------
00021      E0 = 2/SQRT(PI)
00022      E1 = 2**(-3/2)
00023 ------------------------
00024 */
00025 static int K3 = 1;
00026 static double basym,bsum,dsum,f,h,h2,hn,j0,j1,r,r0,r1,s,sum,t,t0,t1,u,w,w0,z,z0,
00027     z2,zn,znm1;
00028 static int i,im1,imj,j,m,mm1,mmj,n,np1;
00029 static double a0[21],b0[21],c[21],d[21],T1,T2;
00030 /*
00031      ..
00032      .. Executable Statements ..
00033 */
00034     basym = 0.0e0;
00035     if(*a >= *b) goto S10;
00036     h = *a/ *b;
00037     r0 = 1.0e0/(1.0e0+h);
00038     r1 = (*b-*a)/ *b;
00039     w0 = 1.0e0/sqrt(*a*(1.0e0+h));
00040     goto S20;
00041 S10:
00042     h = *b/ *a;
00043     r0 = 1.0e0/(1.0e0+h);
00044     r1 = (*b-*a)/ *a;
00045     w0 = 1.0e0/sqrt(*b*(1.0e0+h));
00046 S20:
00047     T1 = -(*lambda/ *a);
00048     T2 = *lambda/ *b;
00049     f = *a*rlog1(&T1)+*b*rlog1(&T2);
00050     t = exp(-f);
00051     if(t == 0.0e0) return basym;
00052     z0 = sqrt(f);
00053     z = 0.5e0*(z0/e1);
00054     z2 = f+f;
00055     a0[0] = 2.0e0/3.0e0*r1;
00056     c[0] = -(0.5e0*a0[0]);
00057     d[0] = -c[0];
00058     j0 = 0.5e0/e0*erfc1(&K3,&z0);
00059     j1 = e1;
00060     sum = j0+d[0]*w0*j1;
00061     s = 1.0e0;
00062     h2 = h*h;
00063     hn = 1.0e0;
00064     w = w0;
00065     znm1 = z;
00066     zn = z2;
00067     for(n=2; n<=num; n+=2) {
00068         hn = h2*hn;
00069         a0[n-1] = 2.0e0*r0*(1.0e0+h*hn)/((double)n+2.0e0);
00070         np1 = n+1;
00071         s += hn;
00072         a0[np1-1] = 2.0e0*r1*s/((double)n+3.0e0);
00073         for(i=n; i<=np1; i++) {
00074             r = -(0.5e0*((double)i+1.0e0));
00075             b0[0] = r*a0[0];
00076             for(m=2; m<=i; m++) {
00077                 bsum = 0.0e0;
00078                 mm1 = m-1;
00079                 for(j=1; j<=mm1; j++) {
00080                     mmj = m-j;
00081                     bsum += (((double)j*r-(double)mmj)*a0[j-1]*b0[mmj-1]);
00082                 }
00083                 b0[m-1] = r*a0[m-1]+bsum/(double)m;
00084             }
00085             c[i-1] = b0[i-1]/((double)i+1.0e0);
00086             dsum = 0.0e0;
00087             im1 = i-1;
00088             for(j=1; j<=im1; j++) {
00089                 imj = i-j;
00090                 dsum += (d[imj-1]*c[j-1]);
00091             }
00092             d[i-1] = -(dsum+c[i-1]);
00093         }
00094         j0 = e1*znm1+((double)n-1.0e0)*j0;
00095         j1 = e1*zn+(double)n*j1;
00096         znm1 = z2*znm1;
00097         zn = z2*zn;
00098         w = w0*w;
00099         t0 = d[n-1]*w*j0;
00100         w = w0*w;
00101         t1 = d[np1-1]*w*j1;
00102         sum += (t0+t1);
00103         if(fabs(t0)+fabs(t1) <= *eps*sum) goto S80;
00104     }
00105 S80:
00106     u = exp(-bcorr(a,b));
00107     basym = e0*t*u*sum;
00108     return basym;
00109 } /* END */