[Top][All Lists]

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

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

From: Alexey Balakin
Subject: [Bug-gsl] Discrete Hankel Transform (possible bug)
Date: Tue, 2 Feb 2010 08:12:53 +0300


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];

        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

reply via email to

[Prev in Thread] Current Thread [Next in Thread]