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 }
