Doxygen Source Code Documentation

Main Page Alphabetical List Data Structures File List Data Fields Globals Search
cdf_08.c

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      ASYMPTOTIC EXPANSION FOR IX(A,B) WHEN A IS LARGER THAN B.
00007      THE RESULT OF THE EXPANSION IS ADDED TO W. IT IS ASSUMED
00008      THAT A .GE. 15 AND B .LE. 1.  EPS IS THE TOLERANCE USED.
00009      IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS.
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      .. 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 */
AFNI/NIfTI Server

Sections

Personal tools

Navigation

Doxygen Source Code Documentation

cdf_08.c
