Re: [Help-gsl] legendre poly deriv

[Patrick Alken]

I think the older routine from 1.16 has the same issue. For Legendre
derivatives, I believe the current algorithm divides by sin(theta), so
at the poles there is a singularity, therefore the routine fails for x=1
or x=-1 as you found.

The following paper discusses the problem and proposes an algorithm to
fix it:

http://www.sciencedirect.com/science/article/pii/S1464189500001010

It is not currently implemented in GSL but if you're interested in
working on it I would be glad to accept a patch for this

Patrick

On 11/10/2017 07:09 AM, Li Dong wrote:
> Hi,
>
> First of all, thanks for this great library! This is my first post in this
> group. Sorry if there is any rule I fail to follow.
>
> I was using gsl 1.16 and recently upgraded to 2.3. I found
> gsl_sf_legendre_Plm_deriv_array is replaced by gsl_sf_legendre_deriv_array.
> Maybe in most cases, this replacement is just fine. However, I find that in
> certain cases there is no way to get a result from current
> gsl_sf_legendre_deriv_array.
>
> In 1.16, the following code works.
> double L[5], DL[5]; // here 5 is just an arbitrary large number to
> initialize the array
> int lmax = 3, m = 2;
> double x=-1.0;
> gsl_sf_legendre_Plm_deriv_array (lmax, m, x, L, DL); // If m=1, this line
> doesn't work since x=-1.
>
> In 2.3, it needs to be written as:
> double L[5], DL[5]; \\ 5 is again arbitrary
> int lmax = 3;
> double x = -1.0;
> gsl_sf_legendre_deriv_array (GSL_SF_LEGENDRE_NONE, lmax, x, L, DL); // this
> line won't work since we calculate all 0<=m<=lmax, it breaks when
> calculates m=1 and x=-1
>
> I wonder if there is a way around to implement this?
>
> Thanks,
> Li