00001 00005 #include <bpm/bpm_messages.h> 00006 #include <bpm/bpm_nr.h> 00007 00011 int nr_gcf(double *gammcf, double a, double x, double *gln) { 00012 00013 int i; 00014 double an,b,c,d,del,h; 00015 00016 *gln = nr_gammln(a); 00017 if( *gln == -DBL_MAX ) { 00018 bpm_error( "nr_gammln failed in nr_gcf(...)", __FILE__, __LINE__ ); 00019 return BPM_FAILURE; 00020 } 00021 00022 b=x+1.0-a; 00023 c=1.0/GCF_FPMIN; 00024 d=1.0/b; 00025 h=d; 00026 for (i=1;i<=GCF_ITMAX;i++) 00027 { 00028 an = -i*(i-a); 00029 b +=2.0; 00030 d=an*d+b; 00031 if (ABS(d) < GCF_FPMIN) d=GCF_FPMIN; 00032 c=b+an/c; 00033 if (ABS(c) < GCF_FPMIN) c=GCF_FPMIN; 00034 d=1.0/d; 00035 del=d*c; 00036 h *= del; 00037 if (ABS(del-1.0) < GCF_EPS) break; 00038 } 00039 00040 if (i > GCF_ITMAX) { 00041 bpm_warning( "A too large, GCF_ITMAX too small in nr_gcf(...)", 00042 __FILE__, __LINE__ ); 00043 } 00044 00045 *gammcf = exp(-x+a*log(x)-(*gln))*h; 00046 00047 return BPM_SUCCESS; 00048 }