diff options
Diffstat (limited to 'gsl-1.9/dht/test.c')
| -rw-r--r-- | gsl-1.9/dht/test.c | 193 |
1 files changed, 193 insertions, 0 deletions
diff --git a/gsl-1.9/dht/test.c b/gsl-1.9/dht/test.c new file mode 100644 index 0000000..1ebbd54 --- /dev/null +++ b/gsl-1.9/dht/test.c @@ -0,0 +1,193 @@ +/* dht/test_dht.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman + * + * 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: G. Jungman + */ +#include <config.h> +#include <stdlib.h> +#include <stdio.h> +#include <math.h> +#include <gsl/gsl_ieee_utils.h> +#include <gsl/gsl_test.h> +#include <gsl/gsl_dht.h> + + +/* Test exact small transform. + */ +int +test_dht_exact(void) +{ + int stat = 0; + double f_in[3] = { 1.0, 2.0, 3.0 }; + double f_out[3]; + gsl_dht * t = gsl_dht_new(3, 1.0, 1.0); + gsl_dht_apply(t, f_in, f_out); + + /* Check values. */ + if(fabs( f_out[0]-( 0.375254649407520))/0.375254649407520 > 1.0e-14) stat++; + if(fabs( f_out[1]-(-0.133507872695560))/0.133507872695560 > 1.0e-14) stat++; + if(fabs( f_out[2]-( 0.044679925143840))/0.044679925143840 > 1.0e-14) stat++; + + + /* Check inverse. + * We have to adjust the normalization + * so we can use the same precalculated transform. + */ + gsl_dht_apply(t, f_out, f_in); + f_in[0] *= 13.323691936314223*13.323691936314223; /* jzero[1,4]^2 */ + f_in[1] *= 13.323691936314223*13.323691936314223; + f_in[2] *= 13.323691936314223*13.323691936314223; + + /* The loss of precision on the inverse + * is a little surprising. However, this + * thing is quite tricky since the band-limited + * function represented by the samples {1,2,3} + * need not be very nice. Like in any spectral + * application, you really have to have some + * a-priori knowledge of the underlying function. + */ + if(fabs( f_in[0]-1.0)/1.0 > 2.0e-05) stat++; + if(fabs( f_in[1]-2.0)/2.0 > 2.0e-05) stat++; + if(fabs( f_in[2]-3.0)/3.0 > 2.0e-05) stat++; + + gsl_dht_free(t); + + return stat; +} + + + +/* Test the transform + * Integrate[x J_0(a x) / (x^2 + 1), {x,0,Inf}] = K_0(a) + */ +int +test_dht_simple(void) +{ + int stat = 0; + int n; + double f_in[128]; + double f_out[128]; + gsl_dht * t = gsl_dht_new(128, 0.0, 100.0); + + for(n=0; n<128; n++) { + const double x = gsl_dht_x_sample(t, n); + f_in[n] = 1.0/(1.0+x*x); + } + + gsl_dht_apply(t, f_in, f_out); + + /* This is a difficult transform to calculate this way, + * since it does not satisfy the boundary condition and + * it dies quite slowly. So it is not meaningful to + * compare this to high accuracy. We only check + * that it seems to be working. + */ + if(fabs( f_out[0]-4.00)/4.00 > 0.02) stat++; + if(fabs( f_out[5]-1.84)/1.84 > 0.02) stat++; + if(fabs(f_out[10]-1.27)/1.27 > 0.02) stat++; + if(fabs(f_out[35]-0.352)/0.352 > 0.02) stat++; + if(fabs(f_out[100]-0.0237)/0.0237 > 0.02) stat++; + + gsl_dht_free(t); + + return stat; +} + + +/* Test the transform + * Integrate[ x exp(-x) J_1(a x), {x,0,Inf}] = a F(3/2, 2; 2; -a^2) + */ +int +test_dht_exp1(void) +{ + int stat = 0; + int n; + double f_in[128]; + double f_out[128]; + gsl_dht * t = gsl_dht_new(128, 1.0, 20.0); + + for(n=0; n<128; n++) { + const double x = gsl_dht_x_sample(t, n); + f_in[n] = exp(-x); + } + + gsl_dht_apply(t, f_in, f_out); + + /* Spot check. + * Note that the systematic errors in the calculation + * are quite large, so it is meaningless to compare + * to a high accuracy. + */ + if(fabs( f_out[0]-0.181)/0.181 > 0.02) stat++; + if(fabs( f_out[5]-0.357)/0.357 > 0.02) stat++; + if(fabs(f_out[10]-0.211)/0.211 > 0.02) stat++; + if(fabs(f_out[35]-0.0289)/0.0289 > 0.02) stat++; + if(fabs(f_out[100]-0.00221)/0.00211 > 0.02) stat++; + + gsl_dht_free(t); + + return stat; +} + + +/* Test the transform + * Integrate[ x^2 (1-x^2) J_1(a x), {x,0,1}] = 2/a^2 J_3(a) + */ +int +test_dht_poly1(void) +{ + int stat = 0; + int n; + double f_in[128]; + double f_out[128]; + gsl_dht * t = gsl_dht_new(128, 1.0, 1.0); + + for(n=0; n<128; n++) { + const double x = gsl_dht_x_sample(t, n); + f_in[n] = x * (1.0 - x*x); + } + + gsl_dht_apply(t, f_in, f_out); + + /* Spot check. This function satisfies the boundary condition, + * so the accuracy should be ok. + */ + if(fabs( f_out[0]-0.057274214)/0.057274214 > 1.0e-07) stat++; + if(fabs( f_out[5]-(-0.000190850))/0.000190850 > 1.0e-05) stat++; + if(fabs(f_out[10]-0.000024342)/0.000024342 > 1.0e-04) stat++; + if(fabs(f_out[35]-(-4.04e-07))/4.04e-07 > 1.0e-03) stat++; + if(fabs(f_out[100]-1.0e-08)/1.0e-08 > 0.25) stat++; + + gsl_dht_free(t); + + return stat; +} + + +int main() +{ + gsl_ieee_env_setup (); + + gsl_test( test_dht_exact(), "Small Exact DHT"); + gsl_test( test_dht_simple(), "Simple DHT"); + gsl_test( test_dht_exp1(), "Exp J1 DHT"); + gsl_test( test_dht_poly1(), "Poly J1 DHT"); + + exit (gsl_test_summary()); +} |