//  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 [v, vdot] = CL__iers_seriesEval(st, t, arg, argdot, cder)
// Evaluates a series (polynomial + non polynomial parts) 
//
// Calling Sequence
// [v, vdot] = CL__iers_seriesEval(st, t, arg, argdot, cder)
//
// Description
// <itemizedlist><listitem>
// <p>st is a structure containing the fields: </p>
// <p>- cpol: polynomial coefficients (deg 0, 1, ...) (1xNpol)</p>
// <p>- l = list of structures of "harmonic" terms </p>
// <p>l(k) is a structure containing the fields: </p>
// <p>ccos: coefficients of cosine terms [Qx1]</p>
// <p>csin: coefficients of sine terms [Qx1] </p>
// <p>carg: coefficients of args in sin() or cos() expressions [QxP] (number of columns of carg = number of rows of arg) </p>
// <p>The result is evaluated by: </p>
// <p>s = sum_over_k>=0 { cpol_k * t^k }  </p>
// <p>  + sum_over_k>=0 { t^k * (sum(ccos_k * cos(carg_k.*arg)) + sum(csin_k * sin(carg_k.*arg)) } </p>
// <p>with: </p>
// <p>ccos_k = l(k).ccos </p>
// <p>csin_k = l(k).csin </p>
// <p>carg_k = l(k).carg </p>
// <p></p>
// <p>All units in the structure must be consistent (and consistent with the arguments t and arg). </p>
// <p></p></listitem>
// </itemizedlist>
//
// Parameters
// t: Time (1xN)
// series: (struct) Structure describing the series (see above) (1x1)
// arg: argument for sin and cos terms (rad) (PxN)
// argdot: derivative of arg (rad/unit of t) (PxN)
// cder: (boolean, optional) Option to compute derivatives. If cder is %f, the derivatives will be set to []. Default is %t. (1x1)
// v: result of evaluation (1xN)
// vdot: time derivative (1xN)
//
// Authors
// CNES - DCT/SB
//

// Declarations:

// Code:
if (~exists("cder", "local")); cder=%t; end
if (argn(1) <= 1); cder = %f; end

// tp: powers of t: [t.^0; t.^1; t.^2 ...]
// n: number of time coef in polynom or series
n = max(size(st.cpol,2), size(st.l)); 
tp = (ones(n,1) * t) .^ ((0:n-1)' * ones(t));  

// Default output values
v = zeros(t);
vdot = [];  
if (cder)
  vdot = zeros(t);
end

// Polynomial part

for k = 1 : size(st.cpol,2)
  v = v + st.cpol(k) * tp(k,:);  
  
  if (cder & k > 1)
    vdot = vdot + (k-1) * st.cpol(k) * tp(k-1,:); 
  end
end


// Non polynomial part
ones1N = ones(t); 

for k = 1 : size(st.l)

  // loop over blocks (optimization)
  for (i = 1 : size(st.l(k).ind))
    I = st.l(k).ind(i).I; 
    J = st.l(k).ind(i).J; 

    ARG = st.l(k).carg(I,J) * arg(J,:);
    A = sum( (st.l(k).ccos(I)*ones1N) .* cos(ARG) + ..
             (st.l(k).csin(I)*ones1N) .* sin(ARG), "r"); 
    v = v + A .* tp(k,:);

    if (cder)
      ARGDOT = st.l(k).carg(I,J) * argdot(J,:); 
      Adot = sum( -(st.l(k).ccos(I)*ones1N) .* ARGDOT .* sin(ARG) + ..
                   (st.l(k).csin(I)*ones1N) .* ARGDOT .* cos(ARG), "r"); 
      vdot = vdot + Adot .* tp(k,:); 
      if (k > 1)
        vdot = vdot + (k-1) * A .* tp(k-1,:); 
      end
    end
  
  end
  
end


endfunction