## 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 != 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"
      scale1 = sqrt (n+0.5);
    case "sch"
      scale1 = sqrt (2);
    case "unnorm"
      scale1 = 1;
    otherwise
      print_usage ();
  endswitch
  scale2 = x .* scale1;
  ## 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 
    LPM = zeros(n+2-m,length(x));
    LPM(1,:) = scale1;
    LPM(2,:) = (2*m-1) .* scale2;
    for l = m+1:n 
      a = (2*l-1) .* x .* LPM(l-m+1,:);
      b = (l+m-2) .* LPM(l-m,:);
      LPM(l-m+2,:) = ((2*l-1) .* x .* LPM(l-m+1,:)...
		      -(l+m-2) .* LPM(l-m,:))/(l-m+1); 
    endfor
    L(m,:) = LPM(n+2-m,:);
    if (strcmp (normalization, "unnorm")) 
      scale1 = -scale1 * (2*m-1); 
      scale2 = -scale2 * (2*m-1); 
    else
      ## normalization == "sch" or normalization == "norm" 
      scale1 = scale1 / sqrt ((n-m+1)*(n+m))*(2*m-1); 
      scale2 = scale2 / sqrt ((n-m+1)*(n+m))*(2*m-1); 
    endif
    scale1 = scale1 .* sqrt(1-x.^2); 
    scale2 = scale2 .* sqrt(1-x.^2);
  endfor
  L(n+1,:) = scale1; 
  if (strcmp (normalization, "sch")) 
    L(1,:) = L(1,:) / sqrt (2); 
  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);