On 08/27/2010 06:47 AM, Brian Gough wrote:
At Wed, 25 Aug 2010 23:41:14 -0400,
Alexey A. Illarionov wrote:
It seems that I find a bug in gsl_sf_coulomb_wave_FG_e. For large values
of lambda and small values of x, with eta = 0 it produces nan values
without raising of the flag.
Here the sample cases (see attachment with example)
gsl_sf_coulomb_wave_FG_e(0.,1.2693881947287221e-07, 37, 1, ...)
gsl_sf_coulomb_wave_FG_e(0.,5.911135077240691e-07, 39, 1, ...)
gsl_sf_coulomb_wave_FG_e(0.,6.118185507743884e-07, 40, 1, ...)
and lots of others.
Thanks for the bug report. I'm adding a test case for this. Can you
confirm if the following values are correct (I'm not working with
these functions often):
lam_F = 37.0;
eta = 0.0;
x = 1.2693881947287221e-07;
k_G = 1;
gsl_sf_coulomb_wave_FG_e(eta, x, lam_F, k_G,&F,&Fp,&G,&Gp,&Fe,&Ge);
s = 0;
message_buff[0] = 0;
s += test_sf_check_result(message_buff, F,
-2.020327525317343313380459069e299, TEST_TOL3);
s += test_sf_check_result(message_buff, Fp,
5.729672862364037942289798904e307, TEST_TOL3);
s += test_sf_check_result(message_buff, G,
4.397355593129492554472984282e299, TEST_TOL3);
s += test_sf_check_result(message_buff, Gp, -
1.247095270068213211088454516e308, TEST_TOL3);
I've computed them using Pari as
c(L,n) = { 2^L * exp(-Pi*n/2) * abs(gamma(L+1+I*n)) / gamma(2*L+2)}
cofac(L,n,r) = { I * exp(-I*r)*r^(-L)/(gamma(2*L+2)*c(L,n)) }
FplusIG(L,n,r) = { cofac(L,n,r) *
intnum(t=0,[[1],1],exp(-t)*t^(L-I*n)*(t+2*I*r)^(L+I*n)) }
FplusIG(37,0,1.2693881947287221e-07+x) # for F and F' (real part)
FplusIG(36,0,1.2693881947287221e-07+x) # for G and G' (imaginary part)
Hi Brian,
It looks like that these values are incorrect since both F and Fp should
be very small (correct me if I'm wrong since I'm not familiar with Pari
and gsl test system). If you need verify the values than its better to
use bessel function series.
For eta == 0. and integer lambda coulomb function are just
bessel-ricatti functions, i.e. spherical_bessel functions times argument
F_l(x) = j_l (x) * x
G_l(x) = n_l (x) * x
Thus
F_l(x) = 1/(2l+1)!! x^(l+1) [1 - x^2/(2l+3)/2 + ...]
G_l(x) = (2l+1)!!/(2l+1) x^(-l) [1 + x^2/(2l-1)/2 + ...]
I think Abramowitz, Stegun should have the complete expressions.
--
Best regards, Alexey Alexandrovich Illarionov
_______________________________________________
Bug-gsl mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/bug-gsl
