//  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 [PREC,OM,OMP] = CL_mod_precessionMatrix(K,KP,KPP, flag)
// Precession matrix and rotation vectors (IERS 1996) - DEPRECATED
//
// Calling Sequence
// [PREC,OM,OMP] = CL_mod_precessionMatrix(K,KP,KPP [,flag])
//
// Description
// <itemizedlist><listitem>
// <p>This function is deprecated. </p>
// <p></p></listitem>
// <listitem>
// <p>Computes the precession matrix, rotation vector and its time derivatives.</p>
// <p>Matrix from mean equator/equinox at date t to mean equator/equinox at J2000.</p>
// <p>if P0 is the pole of J2000 equator, and P is the pole of the equator at date t;
// if G0 is the equinox of J2000 and G is the equinox at date 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
// K : [zetaa;tetaa;za] (3xN)
// KP : [zetaap;tetaap;zap] (3xN)
// KPP : [zetaapp;tetaapp;zapp] (3xN)
// flag : (optional) "n" -> matrix only; "p" -> matrix and rotation vector; "s" -> matrix, rotation vector and its time derivatives (defatul depends on number of outputs)
// PREC: precession matrix (3x3xN)
// OM: (optional) rotation vectors (3xN)
// OMP: (optional) first time derivative of OM (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_precessionAngles
// CL_rot_angles2matrix
// CL_rot_angularVelocity
//
// Examples
// jj_tai = [19500:1:20500];
// ss_tai = 43200*ones(jj_tai);
// [K,KP,KPP]=CL_mod_precessionAngles(jj_tai,ss_tai,"s");
// [PREC,OM,OMP] = CL_mod_precessionMatrix(K,KP,KPP,"s");
//

// 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>4)|(rhs<3) 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

//construction de  la matrice de precession

// premiere  rotation zA
// deuxieme  rotation -thetaA
// troisieme rotation zetaA

tab=zeros(K);
tabp=zeros(KP)
tabpp=zeros(KPP);

tab(1,:)=K(3,:);
tabp(1,:)=KP(3,:)
tabpp(1,:)=KPP(3,:)

tab(2,:)  =-K(2,:);
tabp(2,:) =-KP(2,:);
tabpp(2,:)=-KPP(2,:);

tab(3,:)=K(1,:);
tabp(3,:)=KP(1,:);
tabpp(3,:)=KPP(1,:);

PREC = CL_rot_angles2matrix([3,2,3],tab);

//calcul du vecteur vitesse de rotation (OM) et de la derive (OMP)
OM=[]; OMP=[];
if flag=="p" then
 OM = CL_rot_angularVelocity([3,2,3],tab,tabp,tabpp);
elseif flag=="s" then
 [OM,OMP] = CL_rot_angularVelocity([3,2,3],tab,tabp,tabpp);
end

endfunction