diff --git a/src/kww.c b/src/kww.c
index c23e3d8b6b90cdcad371097783d9f4c4e494f53d..e52a2fc94866d4057adf90202c23b4192975f626 100644
--- a/src/kww.c
+++ b/src/kww.c
@@ -1,73 +1,68 @@
-//**************************************************************************//
-//* FRIDA: fast reliable interactive data analysis *//
-//* spefu.cpp: special functions *//
-//* (C) Joachim Wuttke 1990-, v2(C++) 2001-, spefu 2009- *//
-//* http://frida.sourceforge.net *//
-//**************************************************************************//
-
-#include <string>
+/* kww.c
+ *
+ * Copyright (C) 1998, 2009 Joachim Wuttke
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or (at
+ * your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+/* Author:
+ * Joachim Wuttke
+ * www.messen-und-deuten.de
+ *
+ * Purpose:
+ * Calculate Fourier cos or sin transform of Kohlrausch's stretched
+ * exponential function exp(-t^beta). This problem occurs frequently
+ * in the dynamics of disordered systems.
+ *
+ * References:
+ * Dishon et al. J. Res. Natl. Bur. Stand. 90, 27 (1985).
+ * Chung, ..
+ */
+
#include <math.h>
#include <gsl/gsl_sf.h>
-#include <gsl/gsl_integration.h>
-using namespace std;
-
-const double twopi = 6.28318530717959;
-
-//**************************************************************************//
-//* Fourier transform of stretched exponential function *//
-//**************************************************************************//
+#include <gsl/gsl_errno.h>
+#include "kww.h"
-// References:
-// - Dishon et al. J. Res. Natl. Bur. Stand. 90, 27 (1985).
-// - Chung, Stevens Am. J. Phys. 59(11)1024 (1991).
-// Implementation:
-// - qvi.c by SH Chung, revised August, 1995;
-// - in F77, for Ida ("Frida1"), JWu 24aug98
-// - in C++, for Frida2, JWu 5-mar09
+#define twopi 6.28318530717959
-#define stretched_exponential_nwork 20000
-gsl_integration_workspace *stretched_exponential_work =
- gsl_integration_workspace_alloc( stretched_exponential_nwork );
-typedef struct {
- double w, beta;
-} stretched_exponential_integrand_parameter_type;
-
-double strechted_exponential_cos_integrand( double t, void *P )
-{
- // Chung&Stevens Eq 30
- stretched_exponential_integrand_parameter_type *p;
- p = (stretched_exponential_integrand_parameter_type*) P;
- return exp(-exp(p->beta*log(tan(t))))*cos(p->w*tan(t))/cos(t)/cos(t);
-}
-
-double strechted_exponential_sin_integrand( double t, void *P )
-{
- // Chung&Stevens Eq 31
- stretched_exponential_integrand_parameter_type *p =
- (stretched_exponential_integrand_parameter_type*) P;
- return exp(-exp(p->beta*log(tan(t))))*sin(p->w*tan(t))/cos(t)/cos(t);
-}
-
-double stretched_exponential_cos_transform( double w, double beta )
+double kohlrausch_cos_spectrum( const double w, const double beta )
{
static double tolerance = 1e-10;
static int max_iter = 150; // limited by gamma function overflow
- const int M = 8;
+#define 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, ilim, isig;
- double wlim, result, result1, x, rr, tt1, tt2, err;
- gsl_function F;
- stretched_exponential_integrand_parameter_type P;
+ int ilim, i, isig;
+ double wlim, result, result1, x, rr, tt1, tt2;
+
+ if( w<0 )
+ GSL_ERROR( "w<0", GSL_EDOM );
+ if( beta<Blim[0] )
+ GSL_ERROR( "beta too small", GSL_EDOM );
+ if( beta>Blim[M] )
+ GSL_ERROR( "beta too large", GSL_EDOM );
- if( w<0 || beta<Blim[0] || beta>Blim[M] )
- throw string( "kww_cos: invalid argument" );
+ /* determine auxiliary parameter range */
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[ilim] ) { // Dishon et al, Eq (5), convergent
@@ -77,15 +72,12 @@ double stretched_exponential_cos_transform( double w, double beta )
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);
- // 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 )
return result;
isig = -isig;
}
- throw string( "kww_cos_5: not converged" );
+ GSL_ERROR( "large-argument expansion not converged", GSL_EMAXITER );
} else if ( wlim<=Qsup[ilim] ) { // ibid Eq (6), divergent
result = 0;
isig = 1;
@@ -94,33 +86,27 @@ double stretched_exponential_cos_transform( double w, double beta )
x = 2*i-1;
rr = gsl_sf_gamma(x/beta)*exp((x-1)*log(w))/gsl_sf_gamma(x)/beta;
result = result + isig*rr;
- if( result<=0 )
- throw string( "kww_cos_6: r<0" );
tt1=fabs(rr/result);
- if( tt1<tolerance ) //&& i>40 )
+ if( tt1<tolerance )
return result;
if( tt1>tt2 && i>=3 )
return result1;
tt2 = tt1;
isig = -isig;
}
- throw string( "kww_cos_6: not converged" );
+ GSL_ERROR( "small-argument expansion not converged", GSL_EMAXITER );
} else { // numeric integration
- F.function = strechted_exponential_cos_integrand;
- F.params = (void*) &P;
- P.w = w;
- P.beta = beta;
- gsl_integration_qag( &F, 0, twopi/4,
- 1e-14, 1e-12, stretched_exponential_nwork,
- GSL_INTEG_GAUSS15, stretched_exponential_work,
- &result, &err );
- return result;
+
+
+
+ return 0;
}
}
-double stretched_exponential_sin_transform( double w, double beta )
+
+double kohlrausch_sin_spectrum( const double w, const double beta )
{
- const int M = 8;
+#define M 8
const double Vinf[M] = { .0021, .01, .025, .06, .15, .35, .65, 2. };
const double Vsup[M] = { .001, .001, .01, .045, .1, .25, .5, 1. };
const double Blim[M] = { .3, .4, .5, .6, .7, .8, .9, 1. };
diff --git a/src/kww.h b/src/kww.h
index 765ff5cde445a7e791d6a8b11dd4b312ef86ada1..01031745aa25852ec4950a5805d665b3eefe6842 100644
--- a/src/kww.h
+++ b/src/kww.h
@@ -1,2 +1,50 @@
-double stretched_exponential_cos_transform( double w, double beta );
-double stretched_exponential_sin_transform( double w, double beta );
+/* kww.h
+ *
+ * Copyright (C) 1998, 2009 Joachim Wuttke
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or (at
+ * your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+/* Author: Joachim Wuttke */
+
+#ifndef __KWW_H__
+#define __KWW_H__
+
+#undef __BEGIN_DECLS
+#undef __END_DECLS
+#ifdef __cplusplus
+# define __BEGIN_DECLS extern "C" {
+# define __END_DECLS }
+#else
+# define __BEGIN_DECLS /* empty */
+# define __END_DECLS /* empty */
+#endif
+
+__BEGIN_DECLS
+
+/* Cosine Fourier transform of stretched exponential function:
+ * \Integral_{-\infty}^{+\infty} dt cos(x*t) exp(-t^beta).
+ */
+double kohlrausch_cos_spectrum( const double w, const double beta );
+
+/* Sine Fourier transform of stretched exponential function:
+ * \Integral_{-\infty}^{+\infty} dt sin(x*t) exp(-t^beta).
+ * Also called Kohlrausch-Williams-Watts function.
+ */
+double kohlrausch_sin_spectrum( const double w, const double beta );
+
+__END_DECLS
+
+#endif /* __KWW_H__ */