Hi Patrick, Thanks for your reply.
I think the older legendre derivative routine from 1.16 is slightly better than the one from 2.3. The older routine only fails at pole when m=1. Since you can specify m's value in the older routine, so you can get some results even at the pole. For the new routine from 2.3, m can't be specified and I believe it's calculated for 0<=m<=lmax. So no results at the pole can be obtained. I got these points simply from the error messages I got from GSL and I didn't read much about the current algorithm yet. By the way thanks for the article about the alternative way to implement a better legendre routine, I'll read it and see if I can implement in GSL. Thanks, Li On Fri, Nov 10, 2017 at 2:24 AM, Patrick Alken <[email protected]> wrote: > 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 > > > > -- *Postdoctoral* *Institute for Computational Engineering and Sciences (ICES)* *University of Texas, Austin* *http://users.ices.utexas.edu/~ldong/ <http://users.ices.utexas.edu/~ldong/>*
