Doxygen Source Code Documentation
cdf_08.c File Reference
#include "cdflib.h"Go to the source code of this file.
Functions | |
| void | bgrat (double *a, double *b, double *x, double *y, double *w, double *eps, int *ierr) |
Function Documentation
|
||||||||||||||||||||||||||||||||
|
Definition at line 2 of file cdf_08.c. References a, algdiv(), alnrel(), c, gam1(), grat1(), i, l, n2, p, q, r, and v. Referenced by bratio().
00012 {
00013 static double bm1,bp2n,cn,coef,dj,j,l,lnx,n2,nu,p,q,r,s,sum,t,t2,u,v,z;
00014 static int i,n,nm1;
00015 static double c[30],d[30],T1;
00016 /*
00017 ..
00018 .. Executable Statements ..
00019 */
00020 bm1 = *b-0.5e0-0.5e0;
00021 nu = *a+0.5e0*bm1;
00022 if(*y > 0.375e0) goto S10;
00023 T1 = -*y;
00024 lnx = alnrel(&T1);
00025 goto S20;
00026 S10:
00027 lnx = log(*x);
00028 S20:
00029 z = -(nu*lnx);
00030 if(*b*z == 0.0e0) goto S70;
00031 /*
00032 COMPUTATION OF THE EXPANSION
00033 SET R = EXP(-Z)*Z**B/GAMMA(B)
00034 */
00035 r = *b*(1.0e0+gam1(b))*exp(*b*log(z));
00036 r *= (exp(*a*lnx)*exp(0.5e0*bm1*lnx));
00037 u = algdiv(b,a)+*b*log(nu);
00038 u = r*exp(-u);
00039 if(u == 0.0e0) goto S70;
00040 grat1(b,&z,&r,&p,&q,eps);
00041 v = 0.25e0*pow(1.0e0/nu,2.0);
00042 t2 = 0.25e0*lnx*lnx;
00043 l = *w/u;
00044 j = q/r;
00045 sum = j;
00046 t = cn = 1.0e0;
00047 n2 = 0.0e0;
00048 for(n=1; n<=30; n++) {
00049 bp2n = *b+n2;
00050 j = (bp2n*(bp2n+1.0e0)*j+(z+bp2n+1.0e0)*t)*v;
00051 n2 += 2.0e0;
00052 t *= t2;
00053 cn /= (n2*(n2+1.0e0));
00054 c[n-1] = cn;
00055 s = 0.0e0;
00056 if(n == 1) goto S40;
00057 nm1 = n-1;
00058 coef = *b-(double)n;
00059 for(i=1; i<=nm1; i++) {
00060 s += (coef*c[i-1]*d[n-i-1]);
00061 coef += *b;
00062 }
00063 S40:
00064 d[n-1] = bm1*cn+s/(double)n;
00065 dj = d[n-1]*j;
00066 sum += dj;
00067 if(sum <= 0.0e0) goto S70;
00068 if(fabs(dj) <= *eps*(sum+l)) goto S60;
00069 }
00070 S60:
00071 /*
00072 ADD THE RESULTS TO W
00073 */
00074 *ierr = 0;
00075 *w += (u*sum);
00076 return;
00077 S70:
00078 /*
00079 THE EXPANSION CANNOT BE COMPUTED
00080 */
00081 *ierr = 1;
00082 return;
00083 } /* END */
|
