Re: [Help-gsl] legendre poly deriv

[Li Dong]

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 <address@hidden> 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/>*