normalized ALF (Assotiated Legendre Function)

[Marco Caliari]

Dear all,

please consider the cleaner enclosed script.

I tried to run Ben's example in Maple for ten minutes, with no success. I 
don't know how to give a vector in input to LegendreP.

Best regards,

Marco

On Tue, 12 Feb 2008, Ben Abbott wrote:

>
> On Feb 12, 2008, at 1:44 AM, Dmitri A. Sergatskov wrote:
>
>> FWIW: Marko's script seems to agree with "legendre_Plm"
>> wrapper of GSL's function (gsl package from octave-forge).
>> 
>> Sincerely,
>> 
>> Dmitri.
>> --
>
>
> That is good news!
>
> It would be nice to do a direct comparison for the example below.
>
> 	result_matlab = legendre (80, [-1:0.1:1]);
>
> Any chance anyone who knows c/c++ has the inclination to write a short 
> program using gsl to calculate that example and post the results?
>
> In addition, I noticed the comment below in the reference manual for gsl,
>
> "The following functions compute the associated Legendre Polynomials 
> P_l^m(x). Note that this function grows combinatorially with l and can 
> overflow for l larger than about 150. There is no trouble for small m, but 
> overflow occurs when m and l are both large. Rather than allow overflows, 
> these functions refuse to calculate P_l^m(x) and return GSL_EOVRFLW when they 
> can sense that l and m are too big."
>
> It would be more proper to restrict the issuing of the warning for the cases 
> where overflow actually occurs. For example, even this trivial calculation 
> warns of overflow.
>
> 	octave:73> legendre(0,0)
> 	warning: legendre is unstable for higher orders
> 	ans =  1
>
> I'll take a look at modifying the script to catch overflows and issue the 
> warning only when such occurs.
>
> Ben
-------------- next part --------------
## Copyright (C) 2008 Marco Caliari
##
## This file is part of Octave.
##
## Octave is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or (at
## your option) any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING.  If not, see
## <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn {Function File} {@var{L} =} legendre (@var{n}, @var{X})
## @deftypefnx {Function File} {@var{L} =} legendre (@var{n}, @var{X}, "unnorm")
##
## Legendre Function of degree n and order m
## where all values for m = 0.. at var{n} are returned.
## @var{n} must be an integer.
## The return value has one dimension more than @var{x}.
##
## @example
## The Legendre Function of degree n and order m
##
## @group
##  m        m       2  m/2   d^m
## P(x) = (-1) * (1-x  )    * ----  P (x)
##  n                         dx^m   n
## @end group
##
## with:
## Legendre polynomial of degree n
##
## @group
##           1     d^n   2    n
## P (x) = ------ [----(x - 1)  ] 
##  n      2^n n!  dx^n
## @end group
##
## legendre(3,[-1.0 -0.9 -0.8]) returns the matrix
##
## @group
##  x  |   -1.0   |   -0.9   |  -0.8
## ------------------------------------
## m=0 | -1.00000 | -0.47250 | -0.08000
## m=1 |  0.00000 | -1.99420 | -1.98000
## m=2 |  0.00000 | -2.56500 | -4.32000
## m=3 |  0.00000 | -1.24229 | -3.24000 
## @end group
## @end example
##
## @deftypefnx {Function File} {@var{L} =} legendre (@var{n}, @var{X}, "sch")
##
## Computes the Schmidt semi-normalized associated Legendre function.
## The Schmidt semi-normalized associated Legendre function is related
## to the unnormalized Legendre functions by
##
## @example
## For Legendre functions of degree n and order 0
##
## @group
##   0       0
## SP (x) = P (x)
##   n       n
## @end group
##
## For Legendre functions of degree n and order m
##
## @group
##   m       m          m    2(n-m)! 0.5
## SP (x) = P (x) * (-1)  * [-------]
##   n       n               (n+m)!
## @end group
## @end example
##
## @deftypefnx {Function File} {@var{L} =} legendre (@var{n}, @var{X}, "norm")
##
## Computes the fully normalized associated Legendre function.
## The fully normalized associated Legendre function is related
## to the unnormalized Legendre functions by
##
## @example
## For Legendre functions of degree n and order m
##
## @group
##   m       m          m    (n+0.5)(n-m)! 0.5
## NP (x) = P (x) * (-1)  * [-------------]
##   n       n                   (n+m)!    
## @end group
## @end example
##
## @end deftypefn

## Author:	Marco Caliari <marco.caliari at univr.it>
function L = legendre(n,x,varargin)
x = x(:).';
L = zeros(n+1,length(x));
if (nargin == 2)
  normalization = "unnorm";
else
  normalization = varargin{1};
end
if (strcmp(normalization,"norm"))
  scale = sqrt(n+0.5);
elseif (strcmp(normalization,"sch"))
  scale = sqrt(2);
elseif (strcmp(normalization,"unnorm"))
  scale = 1;
else
  error("Normalization option not recognized.")
end
% Based on the first recurrence relation found in 
% http://en.wikipedia.org/wiki/Associated_Legendre_function
for m = 1:n
  LP = ones(n+1,length(x))*scale;
  LP(m,:) = prod(1:2:2*(m-1)-1)*LP(m,:).*sqrt(1-x.^2).^(m-1);
  if (m == n+1)
    L(m,:) = LP(m,:);
    break
  end  
  LP(m+1,:) = (2*m-1).*x.*LP(m,:);
  for l = m+1:n
    LP(l+1,:) = ((2*(l-1)+1).*x.*LP(l,:)-(l-1+m-1).*LP(l-1,:))./(l-m+1);
  end
  L(m,:) = LP(n+1,:);
  if (strcmp(normalization,"norm"))
    scale = scale/sqrt((n-m+1)*(n+m));
  elseif (strcmp(normalization,"sch"))
    scale = scale/sqrt((n-m+1)*(n+m));
  else
    scale = -scale;
  end
end
L(n+1,:) = prod(1:2:2*n-1)*scale.*sqrt(1-x.^2).^n;
if (strcmp(normalization,"sch"))
  L(1,:) = L(1,:)/sqrt(2);
end