corr

src/spefu.cpp

@@ -153,42 +153,43 @@ const double twopi = 6.28318530717959;
 
 double stretched_exponential_cos_transform( double w, double beta )
 {
-    static double tolerance = 1e-5;
+    static double tolerance = 1e-10;
+    static int max_iter = 150; // limited by gamma function overflow
     const int M = 8;
     const double Qinf[M] = { .0010, .030, .03, .15, .25, .45, .85, 2. };
     const double Qsup[M] = { .0001, .002, .01, .05, .1,  .25, .4,  .5 };
     const double Blim[M+1] = { .2, .3, .4, .5, .6, .7, .8, .9, 1. };
 
-    int i, isig;
+    int ilim, i, isig;
     double wlim, result, result1, x, rr, tt1, tt2;
 
     if( w<0 || beta<Blim[0] || beta>Blim[M] )
         throw string( "kww_cos: invalid argument" );
-    for( i=0; i<M-1; ++i )
-        if( beta<=Blim[1+i] )
+    for( ilim=0; ilim<M-1; ++ilim )
+        if( beta<=Blim[1+ilim] )
             break;
-
     wlim = w * gsl_sf_gamma(1/beta) / beta; // to use original limits {24jan00}
 
-    if      ( wlim>=Qinf[i] ) { // Dishon et al, Eq (5), convergent
+    if      ( wlim>=Qinf[ilim] ) { // Dishon et al, Eq (5), convergent
         result = 0;
         isig = 1;
-        for( i=1; i<150; ++i ){
+        for( i=1; i<max_iter; ++i ){
             x = i*beta+1;
             rr = gsl_sf_gamma(x)/gsl_sf_gamma(i+1.)/exp(x*log(w));
             result = result + isig*rr*sin(twopi/4*i*beta);
-            if( result<=0 )
-                throw string( "kww_cos_5: r<0" );
+            // printf( "%9.4g %3i %9.4g %9.4g\n", w, i, rr, result );
+//            if( i>7 && result<=0 )
+//                throw string( "kww_cos_5: r<0" );
             tt1=fabs(rr/result);
-            if( tt1<tolerance && i>40 )
+            if( tt1<tolerance ) //&& i>40 )
                 return result;
             isig = -isig;
         }
         throw string( "kww_cos_5: not converged" );
-    } else if ( wlim<=Qsup[i] ) { // ibid Eq (6), divergent
+    } else if ( wlim<=Qsup[ilim] ) { // ibid Eq (6), divergent
         result = 0;
         isig = 1;
-        for( i=1; i<150; ++i ){
+        for( i=1; i<max_iter; ++i ){
             result1 = result;
             x = 2*i-1;
             rr = gsl_sf_gamma(x/beta)*exp((x-1)*log(w))/gsl_sf_gamma(x)/beta;
@@ -196,7 +197,7 @@ double stretched_exponential_cos_transform( double w, double beta )
             if( result<=0 )
                 throw string( "kww_cos_6: r<0" );
             tt1=fabs(rr/result);
-            if( tt1<tolerance && i>40 )
+            if( tt1<tolerance ) //&& i>40 )
                 return result;
             if( tt1>tt2 && i>=3 )
                 return result1;