From fdfc18a8c74e085a15a91ed455eb1c39ab3b578a Mon Sep 17 00:00:00 2001
From: "Joachim Wuttke (laptop)" <j.wuttke@fz-juelich.de>
Date: Fri, 6 Mar 2009 15:33:34 +0100
Subject: [PATCH] corr
---
src/spefu.cpp | 27 ++++++++++++++-------------
1 file changed, 14 insertions(+), 13 deletions(-)
diff --git a/src/spefu.cpp b/src/spefu.cpp
index 009734d2..51507261 100644
--- a/src/spefu.cpp
+++ b/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;
--
GitLab