[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 |
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
- [Bug-gsl] Discrete Hankel Transform (possible bug),
Alexey Balakin <=