00001
00005 #include <bpm/bpm_messages.h>
00006 #include <bpm/bpm_nr.h>
00007
00011 int nr_gser(double *gamser, double a, double x, double *gln ) {
00012
00013 int n;
00014 double sum, del, ap;
00015
00016
00017 if( a == 0. ) {
00018 bpm_error( "a equals 0 in nr_gser(...)", __FILE__, __LINE__ );
00019 return BPM_FAILURE;
00020 }
00021
00022 *gln = nr_gammln(a);
00023 if( *gln = -DBL_MAX ) {
00024 bpm_error( "nr_gammln failed in nr_gser(...)", __FILE__, __LINE__ );
00025 return BPM_FAILURE;
00026 }
00027
00028 if (x <= 0.0) {
00029 if (x < 0.0 )
00030 bpm_warning("x less than 0 in routine nr_gser(...)", __FILE__, __LINE__ );
00031
00032 *gamser = 0.0;
00033 return BPM_SUCCESS;
00034 } else {
00035 ap = a;
00036 del = sum = 1.0 / a;
00037 for (n=1;n<=GSER_ITMAX;n++) {
00038 ++ap;
00039 del *= x/ap;
00040 sum += del;
00041 if (ABS(del) < ABS(sum) * GSER_EPS ) {
00042
00043 *gamser = sum*exp(-x+a*log(x)-(*gln));
00044 return BPM_SUCCESS;
00045 }
00046 }
00047
00048 bpm_error( "a too large, GSER_ITMAX too small in nr_gser(...)",
00049 __FILE__, __LINE__ );
00050 return BPM_FAILURE;
00051 }
00052
00053
00054 return BPM_FAILURE;
00055 }