Re: [Bug-gsl] bug in gsl_sf_coulomb_wave_FG_e (gsl 1.14 compiled with gcc 4.5.0)

[Alexey A. Illarionov]

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