diff options
Diffstat (limited to 'gsl-1.9/dht/dht.c')
| -rw-r--r-- | gsl-1.9/dht/dht.c | 225 |
1 files changed, 225 insertions, 0 deletions
diff --git a/gsl-1.9/dht/dht.c b/gsl-1.9/dht/dht.c new file mode 100644 index 0000000..2ba6fb1 --- /dev/null +++ b/gsl-1.9/dht/dht.c @@ -0,0 +1,225 @@ +/* dht/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 <gsl/gsl_errno.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_sf_bessel.h> +#include <gsl/gsl_dht.h> + + +gsl_dht * +gsl_dht_alloc (size_t size) +{ + gsl_dht * t; + + if(size == 0) { + GSL_ERROR_VAL("size == 0", GSL_EDOM, 0); + } + + t = (gsl_dht *)malloc(sizeof(gsl_dht)); + + if(t == 0) { + GSL_ERROR_VAL("out of memory", GSL_ENOMEM, 0); + } + + t->size = size; + + t->xmax = -1.0; /* Make it clear that this needs to be calculated. */ + t->nu = -1.0; + + t->j = (double *)malloc((size+2)*sizeof(double)); + + if(t->j == 0) { + free(t); + GSL_ERROR_VAL("could not allocate memory for j", GSL_ENOMEM, 0); + } + + t->Jjj = (double *)malloc(size*(size+1)/2 * sizeof(double)); + + if(t->Jjj == 0) { + free(t->j); + free(t); + GSL_ERROR_VAL("could not allocate memory for Jjj", GSL_ENOMEM, 0); + } + + t->J2 = (double *)malloc((size+1)*sizeof(double)); + + if(t->J2 == 0) { + free(t->Jjj); + free(t->j); + free(t); + GSL_ERROR_VAL("could not allocate memory for J2", GSL_ENOMEM, 0); + } + + return t; +} + +/* Handle internal calculation of Bessel zeros. */ +static int +dht_bessel_zeros(gsl_dht * t) +{ + unsigned int s; + gsl_sf_result z; + int stat_z = 0; + t->j[0] = 0.0; + for(s=1; s < t->size + 2; s++) { + stat_z += gsl_sf_bessel_zero_Jnu_e(t->nu, s, &z); + t->j[s] = z.val; + } + if(stat_z != 0) { + GSL_ERROR("could not compute bessel zeroes", GSL_EFAILED); + } + else { + return GSL_SUCCESS; + } +} + +gsl_dht * +gsl_dht_new (size_t size, double nu, double xmax) +{ + int status; + + gsl_dht * dht = gsl_dht_alloc (size); + + if (dht == 0) + return 0; + + status = gsl_dht_init(dht, nu, xmax); + + if (status) + return 0; + + return dht; +} + +int +gsl_dht_init(gsl_dht * t, double nu, double xmax) +{ + if(xmax <= 0.0) { + GSL_ERROR ("xmax is not positive", GSL_EDOM); + } else if(nu < 0.0) { + GSL_ERROR ("nu is negative", GSL_EDOM); + } + else { + size_t n, m; + int stat_bz = GSL_SUCCESS; + int stat_J = 0; + double jN; + + if(nu != t->nu) { + /* Recalculate Bessel zeros if necessary. */ + t->nu = nu; + stat_bz = dht_bessel_zeros(t); + } + + jN = t->j[t->size+1]; + + t->xmax = xmax; + t->kmax = jN / xmax; + + t->J2[0] = 0.0; + for(m=1; m<t->size+1; m++) { + gsl_sf_result J; + stat_J += gsl_sf_bessel_Jnu_e(nu + 1.0, t->j[m], &J); + t->J2[m] = J.val * J.val; + } + + /* J_nu(j[n] j[m] / j[N]) = Jjj[n(n-1)/2 + m - 1], 1 <= n,m <= size + */ + for(n=1; n<t->size+1; n++) { + for(m=1; m<=n; m++) { + double arg = t->j[n] * t->j[m] / jN; + gsl_sf_result J; + stat_J += gsl_sf_bessel_Jnu_e(nu, arg, &J); + t->Jjj[n*(n-1)/2 + m - 1] = J.val; + } + } + + if(stat_J != 0) { + GSL_ERROR("error computing bessel function", GSL_EFAILED); + } + else { + return stat_bz; + } + } +} + + +double gsl_dht_x_sample(const gsl_dht * t, int n) +{ + return t->j[n+1]/t->j[t->size+1] * t->xmax; +} + + +double gsl_dht_k_sample(const gsl_dht * t, int n) +{ + return t->j[n+1] / t->xmax; +} + + +void gsl_dht_free(gsl_dht * t) +{ + free(t->J2); + free(t->Jjj); + free(t->j); + free(t); +} + + +int +gsl_dht_apply(const gsl_dht * t, double * f_in, double * f_out) +{ + const double jN = t->j[t->size + 1]; + const double r = t->xmax / jN; + size_t m; + size_t i; + + for(m=0; m<t->size; m++) { + double sum = 0.0; + double Y; + for(i=0; i<t->size; i++) { + /* Need to find max and min so that we + * address the symmetric Jjj matrix properly. + * FIXME: we can presumably optimize this + * by just running over the elements of Jjj + * in a deterministic manner. + */ + size_t m_local; + size_t n_local; + if(i < m) { + m_local = i; + n_local = m; + } + else { + m_local = m; + n_local = i; + } + Y = t->Jjj[n_local*(n_local+1)/2 + m_local] / t->J2[i+1]; + sum += Y * f_in[i]; + } + f_out[m] = sum * 2.0 * r*r; + } + + return GSL_SUCCESS; +} + |