[Help-gsl] Discrete Hankel transform issue

[Michael Carley]

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/