//  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 [val] = CL_evalPoly(coeffs,t, nd)
// Evaluation of polynomial or polynomial derivatives
//
// Calling Sequence
// val = CL_evalPoly(coeffs,t [,nd])
//
// Description
// <itemizedlist><listitem>
// <p>Computes the value of a polynomial or its derivatives.</p> 
// <p> The polynomials are defined by the coefficients <b>coeffs</b>. 
// Each row (number i) in coeffs contains the coefficients of the ith polynom. 
// coeffs(i,k) is the coefficient of the term t^(k-1). </p>
// <p> The evaluation is done for each value in <b>t</b>. </p>  
// <p> <b>nd</b> is the derivation order:  </p>
// <p> - nd == 1 => first derivative </p>
// <p> - nd == 2 => second derivative </p>
// <p> - etc... </p>
// </listitem>
// </itemizedlist>
//
// Parameters
// coeffs : Polynomial coefficients. Each row corresponds to a different polynomial. (PxM)
// t : Values for the evaluation (1xN)
// nd : (optional) Derivation order. Default is 0.
//
// Authors
// CNES - DCT/SB
//
// Examples
// // Example 1:
// coeffs = [0,1,2];  // P(x) = x + 2*x^2
// CL_evalPoly(coeffs,[0,1])
// CL_evalPoly(coeffs,[0,1],1) // P'(x) = 1 + 4*x
// CL_evalPoly(coeffs,[0,1],2) // P''(x) = 4;
// CL_evalPoly(coeffs,[0,1],3) // P'''(x) = 0;
//
// // Example 2:
// t=-0.05:0.001:0.05;
// coeffs = rand(500,3,"normal"); // 1000 2nd degree polynoms
// x = CL_evalPoly(coeffs,t);     // polynom values
// x2 = CL_evalPoly(coeffs,t,2);  // second derivative 
//

// Declarations:


// Code:

if (~exists('nd', 'local'))
   nd = 0; 
end

if (size(t,1) <> 1) 
  CL__error("Wrong size for t");
end

if (nd < 0) 
  CL__error("Wrong value for nd");
end

[P,nn]=size(coeffs);

if (nd >= nn) 
   val = zeros(P,size(t,2));
   return; 
end


A = coeffs(:,nd+1:nn);
K = ones(P,1)*[1:nd];

l=1;
if nd > 0 then
  for c=nd+1:nn
    A(:,l) = A(:,l)*prod(c-[1:nd]);
    l=l+1;
  end,
end,

FPOLY=[];
for I=1:P
  FPOLY=[FPOLY; poly(A(I,:),"t","coeff");];
end

val = horner(FPOLY,t);

endfunction