## Copyright (C) 2000, 2006, 2007 Kai Habel

## 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} address@hidden =} legendre (@var{n}, @var{X})

## @deftypefnx {Function File} address@hidden =} legendre (@var{n}, @var{X}, "unnorm")

##

## Legendre Function of degree n and order m

## where all values for m = address@hidden are returned.

## @var{n} must be a scalar in the range [0..255].

## 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} address@hidden =} 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} address@hidden =} 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 <address@hidden>



function L = legendre (n, x, normalization)



  if (nargin < 2 || nargin > 3)

    print_usage ();

  endif

  if nargin == 3

    normalization = lower (normalization);

  else

    normalization = "unnorm";

  endif



  if (! isscalar (n) || n < 0 || n > 255 || n != fix (n))

    error ("n must be a integer between 0 and 255]");

  endif



  if (! isvector (x) || any (x < -1 || x > 1))

    error ("x must be vector in range -1 <= x <= 1");

  endif



  switch normalization

    case "norm"

      scale = sqrt (n+0.5);

    case "sch"

      scale = sqrt (2);

    case "unnorm"

      scale = 1;

    otherwise

      print_usage ();

  endswitch

  ## Based on the recurrence relation below

  ##            m                 m              m

  ## (n-m+1) * P (x) = (2*n+1)*x*P (x)  - (n+1)*P (x)

  ##            n+1               n              n-1

  ## http://en.wikipedia.org/wiki/Associated_Legendre_function 

  underflow = false;

  for m = 1:n 

    LP = ones (n+1, length (x)); 

    LP(m,:) = scale .* LP(m,:); 

    LP(m+1,:) = (2*m-1) .* x .* LP(m,:); 

    for l = m+1:n 

      a = (2*l-1) .* x .* LP(l,:);

      b = (l+m-2) .* LP(l-1,:);

      if (any((a-b) == a & b != 0) || any((a-b) == -b & a != 0))

	a-b

        underflow = true;

      endif

      LP(l+1,:) = (a-b)./(l-m+1); 

      ##LP(l+1,:) = ((2*l-1).*x.*LP(l,:)-(l+m-2).*LP(l-1,:))./(l-m+1); 

    endfor

    L(m,:) = LP(n+1,:); 

    if (strcmp (normalization, "unnorm")) 

      scale = -scale * (2*m-1) .* sqrt (1-x.^2); 

      if (max(abs(scale)) == Inf)

	underflow = true;

      endif

    else

      ## normalization == "sch" or normalization == "norm" 

      scale = scale / sqrt ((n-m+1)*(n+m))*(2*m-1) .* sqrt (1-x.^2); 

    endif

  endfor

  L(n+1,:) = scale; 

  if (strcmp (normalization, "sch")) 

    L(1,:) = L(1,:) / sqrt (2); 

  endif

  if (underflow)

    warning ("legendre: results may be  unstable for high orders.");

  endif



endfunction



%!test

%! result=legendre(3,[-1.0 -0.9 -0.8]);

%! expected = [

%!    -1.00000  -0.47250  -0.08000

%!     0.00000  -1.99420  -1.98000

%!     0.00000  -2.56500  -4.32000

%!     0.00000  -1.24229  -3.24000

%! ];

%! assert(result,expected,1e-5);



%!test

%! result=legendre(3,[-1.0 -0.9 -0.8], "sch");

%! expected = [

%!    -1.00000  -0.47250  -0.08000

%!     0.00000   0.81413   0.80833

%!    -0.00000  -0.33114  -0.55771

%!     0.00000   0.06547   0.17076

%! ];

%! assert(result,expected,1e-5);



%!test

%! result=legendre(3,[-1.0 -0.9 -0.8], "norm");

%! expected = [

%!    -1.87083  -0.88397  -0.14967

%!     0.00000   1.07699   1.06932

%!    -0.00000  -0.43806  -0.73778

%!     0.00000   0.08661   0.22590

%! ];

%! assert(result,expected,1e-5);