help-gsl
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

## [Help-gsl] Discrete Hankel transform issue

 From: Michael Carley Subject: [Help-gsl] Discrete Hankel transform issue Date: Tue, 26 Jul 2011 11:09:28 +0100 (BST) User-agent: Alpine 2.00 (LNX 1167 2008-08-23)

```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 ;
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, g, 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/

```

reply via email to

 [Prev in Thread] Current Thread [Next in Thread]
• [Help-gsl] Discrete Hankel transform issue, Michael Carley <=