//  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 [K,KP,KPP] = CL_mod_precessionAngles(jj_tai,sec_tai, flag)
// Precession angles (Lieske 1976) - DEPRECATED
//
// Calling Sequence
// [K,KP,KPP] = CL_mod_precessionAngles(jj_tai,sec_tai [,flag])
//
// Description
// <itemizedlist><listitem>
// <p>This function is deprecated. </p>
// <p></p></listitem>
// <listitem>
// <p>Computes the precession angles according to Lieske precession theory.</p>
// <p>if P0 is the pole of J2000 equator, and P is the pole at time t;
// if G0 is the equinox of J2000 and G is the equinox at time t:</p>
// <p>zetaa: angle between the meridians (pole P0) P0-P and P0-G0</p>
// <p>thetaa: angle between the equateur J2000 and the equateur t</p>
// <p>za: angle between the meridians (pole P) P-P0 and P-G</p>
// </listitem>
// </itemizedlist>
//
// Parameters
// jj_tai : (integer) Modified (1950.0) Julian day (1xN)
// ss_tai : seconds  (TAI time scale) (1xN)
// flag   : (optional) "n" -> K; "p" -> K and Kp; "s" -> K, Kp, Kpp (default depends on number of outputs)
// K      : precession angles K = [zetaa;tetaa;za] [rad] (3xN)
// KP : (optional) first time derivative of K [rad/s] (3xN)
// KPP : (optional) second time derivative of K [rad/s^2] (3xN)
//
// Authors
// CNES - DCT/SB
//
// Bibliography
// 1) IERS Conventions (1996), Dennis D. McCarthy
// 2) Explanatory Supplement to the Astronomical Almanac, Seidelman (1992)
//
// See also
// CL_mod_nutationAngles
//
// Examples
// jj_tai = [19500:1:20500];
// ss_tai = 43200*ones(jj_tai);
// [K,KP,KPP]=CL_mod_precessionAngles(jj_tai,ss_tai,"s");
// plot(jj_tai'+ss_tai'/86400,K');
//

// Declarations:


// Code:

CL__warnDeprecated(); // deprecated function

[lhs rhs]=argn(0);
if (~exists('flag','local'))
  select lhs
    case 1
      flag='n';
    case 2
      flag='p';
    case 3
      flag='s';
    else
      CL__error("check number of output arguments");
  end
end

if (rhs>3)|(rhs==1) CL__error("check number of input arguments"); end,

if (flag~="s" & flag~="p" & flag~="n") then
  CL__error("value of parameter ''flag'' = "+flag+" unknown");
end

J2000 = 36525. / 2.;
TT_MOINS_TAI = 32.184;
t = (jj_tai - J2000 + (TT_MOINS_TAI+sec_tai)/86400) / 36525;
t = t(:)';

COEFF  = [0. 2306.2181  0.30188  0.017998 0. 0.;
          0. 2004.3109 -0.42665 -0.041833 0. 0.;
          0. 2306.2181  1.09468  0.018203 0. 0.;];

SEC2RAD = (%pi / 180) / 3600;
COEFF=COEFF*SEC2RAD;

//evaluation du polynom
K=CL_evalPoly(COEFF,t,0);
K=modulo(K,2*%pi);
KP=[]; KPP=[];
if (flag == "p" | flag =="s") then
  KP=CL_evalPoly(COEFF,t,1);
  KP=KP/(86400*36525);
end,
if (flag =="s") then
  KPP=CL_evalPoly(COEFF,t,2);
  KPP=KPP/((86400*36525)^2);
end

endfunction