Drexel dragonThe Math ForumDonate to the Math Forum



Search All of the Math Forum:

Views expressed in these public forums are not endorsed by Drexel University or The Math Forum.


Math Forum » Discussions » sci.math.* » sci.math.num-analysis.independent

Topic: Legendre Polynomials normalised to Sin[theta]^m
Replies: 13   Last Post: Jul 1, 2013 12:13 PM

Advanced Search

Back to Topic List Back to Topic List Jump to Tree View Jump to Tree View   Messages: [ Previous | Next ]
Guido Walter Pettinari

Posts: 29
From: Italy & UK
Registered: 3/24/10
Re: Legendre Polynomials normalised to Sin[theta]^m
Posted: Jul 1, 2013 12:13 PM
  Click to see the message monospaced in plain text Plain Text   Click to reply to this topic Reply

Dear Wolfgang,

I modified my Legendre Polynomial routine by setting f=1, and... it worked! Now the routine correctly computes P_lm(x) divided by (1-x*x)^(m/2). Thank you very much for your precious suggestion.

Please find attached the C code to estimate the special function.

Thanks again,
Guido


/**
* Return the value of
* P_lm (x) / (1-x^2)^(m/2)
* where P_lm(x) is the associated Legendre polynomial, computed using a recurrence relation
* algorithm. The returned polynomials include the same normalisation constant found in the
* definition of the spherical harmonics, that is sqrt( (2l+1)/(4pi) * (l-m)!/(l+m)! )
*
* The credits for this function go to Press et al. 2002. ("Numerical Recipes, Third edition",
* pag.294) for the computation of the Legendre polynomial, and to Wolfgang Ehrhardt
* (https://groups.google.com/forum/#!topic/sci.math.num-analysis/13ZmcOtpWzg) for the
* trick to rescale the function by (1-x^2)^(m/2) analitically.
*/
double plegendre_lm_rescaled_analytically (int l, int m, double x) {

int i,ll;
double fact,oldfact,pll,pmm,pmmp1,omx2;

if (m < 0 || m > l || fabs(x) > 1.0) {
printf("ERROR: Bad arguments in routine plegendre_lm\n");
}

pmm=1.0;
/* Compute P_m^m/(1-x*x)^(m/2) */
if (m > 0) {
omx2=1;
/* Uncomment the following line to compute the usual associate Legendre polynomials */
// omx2=(1.0-x)*(1.0+x);
fact=1.0;
for (i=1;i<=m;i++) {
pmm *= omx2*fact/(fact+1.0);
fact += 2.0;
}
}
pmm=sqrt((2*m+1)*pmm/(4.0*_PI_));
if (m & 1)
pmm=-pmm;
if (l == m)
return pmm;
else {
/* Compute P_m+1^m */
pmmp1=x*sqrt(2.0*m+3.0)*pmm;
if (l == (m+1))
return pmmp1;
else {
/* Compute P_l^m, l>m+1 */
oldfact=sqrt(2.0*m+3.0);
for (ll=m+2;ll<=l;ll++) {
fact=sqrt((4.0*ll*ll-1.0)/(ll*ll-m*m));
pll=(x*pmmp1-pmm/oldfact)*fact;
oldfact=fact;
pmm=pmmp1;
pmmp1=pll;
}
return pll;
}
}
}



Il giorno lunedì 1 luglio 2013 12:13:23 UTC+2, Wolfgang Ehrhardt ha scritto:
> On Sat, 29 Jun 2013 13:04:27 -0700 (PDT), Guido Walter Pettinari wrote:
>

> >Thank you for the reference. I have read it but I am afraid that it does no=
>
> >t help me. I have already got a routine to compute the associated Legendre =
>
> >Polynomials from Numerical Recipes; I have also already checked Abramowitz =
>
> >and Stegun, but I could not find a reference to the function that I need.
>
>
>
> You can use the Legendre recurrence relation for varying degree from
>
> {A&S 8.5.3} together with the starting value p0 = P_mm = (-1)^m *
>
> (2m-1)!! * f, {NumRec 6.8.8}, where f = (1-x^2)^(m/2). If you set
>
> f=1 this is IMHO the function you want, in my Pascal library
>
> <http://www.wolfgang-ehrhardt.de/misc_en.html#specfun_units> there is
>
> a function for this where f is externally supplied. If it helps here
>
> is the source code (the used functions should be obvious, HMF =
>
> Abramowitz/Stegun, NR=Numerical recipes in C, 2nd ed., [19] = Boost
>
> library).
>
>
>
> Hope that helps
>
> Wolfgang
>
>
>
> ------------------------
>
> function legendre_plmf(l,m: integer; x,f: extended): extended;
>
> {-Associated Legendre polynomial P_lm(x), f=|1-x^2|^(|m|/2)
>
> calculated externally}
>
> var
>
> p0, p1, pk: extended;
>
> k: integer;
>
> begin
>
> {Ref: HMF [1], ch.8; NR [13] ch.6.8; [19] legendre.hpp/legendre_imp}
>
> if l<0 then l := -l-1; {negative degree: HMF[1] 8.2.1}
>
> if m<0 then begin
>
> {Negative order: P(l,-m) = s*(l-m)!/(l+m)! P(l,m), for s see
>
> below}
>
> k := l+m;
>
> if k<0 then p1 := 0
>
> else p1 := legendre_plmf(l, -m, x, f);
>
> if p1<>0.0 then begin
>
> {calculate Gamma(l+m+1)/Gamma(l-m+1) via sfc_fac if possible}
>
> l := l-m;
>
> if l<MAXGAMX then p0 := sfc_fac(k)/sfc_fac(l)
>
> else p0 := sfc_gamma_ratio(k+1, l+1);
>
> {If |x| <= 1 there is an additional factor s=(-1)^m, but no
>
> factor}
>
> {if |x| > 1, see NIST [30], 14.9.3 vs 14.9.13}
>
> if odd(m) and (abs(x)<=1.0) then p1 := -p1*p0
>
> else p1 := p1*p0;
>
> end;
>
> legendre_plmf := p1;
>
> end
>
> else if m>l then legendre_plmf := 0.0
>
> else if m=0 then legendre_plmf := sfc_legendre_p(l,x)
>
> else begin
>
> {Starting value p0 = P_mm = (-1)^m * (2m-1)!! * f, NR[13] 6.8.8}
>
> p0 := sfc_dfac(2*m-1)*f;
>
> if odd(m) and (x<=1.0) then p0 := -p0;
>
> if m=l then begin
>
> legendre_plmf := p0;
>
> exit;
>
> end;
>
> p1 := x*p0*(2*m + 1);
>
> k := succ(m);
>
> while k<l do begin
>
> {Recurrence relation for varying degree: HMF[1] 8.5.3}
>
> pk := ((2*k+1)*x*p1 - (k+m)*p0)/(k+1-m);
>
> p0 := p1;
>
> p1 := pk;
>
> inc(k);
>
> end;
>
> legendre_plmf := p1;
>
> end;
>
> end;





Point your RSS reader here for a feed of the latest messages in this topic.

[Privacy Policy] [Terms of Use]

© Drexel University 1994-2014. All Rights Reserved.
The Math Forum is a research and educational enterprise of the Drexel University School of Education.