Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
cdf_08.c
Go to the documentation of this file.00001 #include "cdflib.h"
00002 void bgrat(double *a,double *b,double *x,double *y,double *w,
00003 double *eps,int *ierr)
00004
00005
00006
00007
00008
00009
00010
00011
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
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
00033
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
00073
00074 *ierr = 0;
00075 *w += (u*sum);
00076 return;
00077 S70:
00078
00079
00080
00081 *ierr = 1;
00082 return;
00083 }