//  Copyright (c) CNES  2008
//
//  This software is part of CelestLab, a CNES toolbox for Scilab
//
//  This software is governed by the CeCILL  license under French law and
//  abiding by the rules of distribution of free software.  You can  use,
//  modify and/ or redistribute the software under the terms of the CeCILL
//  license as circulated by CEA, CNRS and INRIA at the following URL
//  'http://www.cecill.info'.

function [p, pdot, p2dot] = CL_legendrePoly(n)
// Legendre polynomials 
//
// Calling Sequence
// [p, pdot, p2dot] = CL_legendrePoly(n)
//
// Description
// <itemizedlist><listitem>
// <p>Computes Legendre polynomials.</p>
// <p>Legendre polynomials are defined as follows:</p>
// <p>P0(x) = 1  (degree 0) </p>
// <p>P1(x) = x  (degree 1) </p>
// <p>(k+1) * Pk+1(x) = (2*k+1)*Pk(x) - k * Pk-1(x)  (for k >= 1) </p>
// <p></p>
// <p>The function returns the polynomials of degree <b>n</b>.
// <b>p</b>(k) is the polynomial of degree n(k).</p>
// <p>The function also optionally returns the first and second derivatives of the polynomials.</p>
// <p></p></listitem>
// </itemizedlist>
//
// Parameters
// n: (integer) Degrees of the polynomials in [0,100] (PxQ)
// p: Legendre polynomials (same size as n) 
// pdot: First derivative of the polynomials (same size as p) 
// p2dot: Second derivative of the polynomials (same size as p) 
//
// Authors
// CNES - DCT/SB
//
// Examples
// // Generate Legendre polynomials (degree 0 to 5) and plot them
// n = 0:5; 
// p = CL_legendrePoly(n);
// t = linspace(-1,1,100); 
// scf(); 
// plot(t, horner(p',t)); 
// CL_g_legend(gca(), "degree: " + string(n));

// Declarations:

// Code:

// number of output args (left/right)
[lhs,rhs] = argn(); 

if (rhs <> 1)
  CL__error("Invalid number of input arguments"); 
end

if (find(n < 0 | n > 100 | round(n) <> n) <> [])
  CL__error("Invalid input argument: n"); 
end

// kmax: max number of computed polynomials (min=2)
kmax = max(max(n)+1,2); 
p = zeros(kmax, 1); 

x = poly(0,"x");

// Polynomial of degree 0
p(1) = 1 + 0*x;

// Polynomial of degree 1
p(2) = x; 

for (k = 1 : kmax-2)
  p(k+2) = ((2*k+1)*x*p(k+1) - k * p(k))/(k+1); 
end

// only degrees in n are returned
p = matrix(p(n+1), size(n)); 

if (lhs >= 2)
  pdot = derivat(p); 
end

if (lhs >= 3)
  p2dot = derivat(pdot); 
end

endfunction