//  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_chebyshevPoly(n, kind)
// Chebyshev polynomials (of first or second kind)
//
// Calling Sequence
// [p, pdot, p2dot] = CL_chebyshevPoly(n [,kind])
//
// Description
// <itemizedlist><listitem>
// <p>Computes Chebyshev polynomials of first or second kind.</p>
// <p>The polynomials are defined as follows:</p>
// <p>P0(x) = 1  (degree 0)</p>
// <p>P1(x) = kind * x, with kind = 1 or 2  (degree 1)</p>
// <p>Pn(x) = 2 * x * Pn-1(x) - Pn-2(x), for n >= 2 </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)
// kind: (optional) 1 for first kind, 2 for second kind. Default is 1 (1x1)
// p: Chebyshev 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
//
// See also
// CL_chebyshevEval
//
// Examples
// // Generate Chebyshev polynomials (degree 0 to 5, first kind) and plot them
// n = 0:5; 
// p = CL_chebyshevPoly(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 | rhs > 2)
  CL__error("Invalid number of input arguments"); 
end

if (~exists("kind", "local")); kind = 1; end

if ((kind <> 1 & kind <> 2) | size(kind,"*") <> 1)
  CL__error("Invalid input argument: kind"); 
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) = kind * x;  

// Recurrence formula
for k = 3:kmax
  p(k) = 2 * x * p(k-1) - p(k-2); 
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