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
