Hello all,
I am using the DHT functions in GSL and having some trouble with
the inversion. Since I actually want the coefficients of the expansion, I
am testing my code with a `manual inversion' (see source code below). The
problem is that this gives me a version of the original input scaled on
(X^2+1), except at the first point, where the scaling factor is (X^2).
Using the DHT function to invert the result gives what I expect. Any ideas
what I am doing wrong?
#include <math.h>
#include <stdio.h>
#include <gsl/gsl_dht.h>
#include <gsl/gsl_sf_bessel.h>
int hankel_invert(gsl_dht *H, int m, double X, double *g, int n,
double *f)
{
double sc, k, x ;
int i, j ;
for ( j = 0 ; j < n ; j ++ ) f[0] = 0.0 ;
for ( i = 0 ; i < n ; i ++ ) {
k = gsl_dht_k_sample(H, i) ;
sc = gsl_sf_bessel_Jn(m+1, k*X) ;
for ( j = 0 ; j < n ; j ++ ) {
x = gsl_dht_x_sample(H, j) ;
f[j] += g[i]*2.0*gsl_sf_bessel_Jn(m, k*x)/(sc*sc) ;
}
}
return 0 ;
}
int main(int argc, char **argv)
{
double f[1024], g[1024], x, X ;
int n = 32, i, p ;
gsl_dht *H ;
H = gsl_dht_alloc(n) ;
p = 2 ; X = 4.5 ;
gsl_dht_init(H, (double)p, X) ;
for ( i = 0 ; i < n ; i ++ ) {
x = gsl_dht_x_sample(H, i) ;
f[i] = x*(X-x) ;
}
gsl_dht_apply(H, f, g) ;
hankel_invert(H, p, X, g, n, f) ;
for ( i = 0 ; i < n ; i ++ ) {
fprintf(stdout, "%d %1.16e %1.16e\n", i, gsl_dht_x_sample(H, i), f[i]);
}
return 0 ;
}
--
`To tell the truth, let us be honest at least, it is some considerable
time since I last knew what I was talking about.'
No MS attachments: http://www.gnu.org/philosophy/no-word-attachments.html
Home page: http://people.bath.ac.uk/ensmjc/
_______________________________________________
Help-gsl mailing list
Help-gsl@gnu.org
https://lists.gnu.org/mailman/listinfo/help-gsl