[Bug-gsl] Discrete Hankel Transform (possible bug)

[Alexey Balakin]

Hi,

I think a formula for discrete Hankel transform (DHT) contain a bug.
The following code

void identity(int n, double *a, double *b)
{
        gsl_dht *dht = gsl_dht_new(n,0,1);
        double *t=new double[n];
        gsl_dht_apply(dht,a,t);
        gsl_dht_apply(dht,t,b);

        double f = gsl_sf_bessel_zero_J0(n+1);
        for(long i=0;i<n;i++)   b[i]*=f*f;      // note here!!!

        delete []t;     gsl_dht_free(dht);
}

should give the same arrays b[i]=a[i] (at least with accuracy of
numeric error). But it is possible only if one introduce additional
multiplier f=j_{\nu, M} (in notation of GSL documentation). Note, here
I use f*f because I do 2 transforms (like forward and backward).

So, I suggest to note this in documentation (as GSL "feature"), or
change the formula and the code of gsl_dht_apply() function to
additional multiplication by j_{\nu, M}.

One more comment. I think that it will be nice if gsl_dht_apply() can
handle the same arrays for input and output :) .

-- 
All the best,
Alexey Balakin