//  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 [NUT,OM,OMP] = CL_mod_nutationMatrix(K,KP,KPP, flag)
// Nutation matrix and rotation vectors (IERS 1996) - DEPRECATED
//
// Calling Sequence
// [NUT,OM,OMP] = CL_mod_nutationMatrix(K,KP,KPP [,flag])
//
// Description
// <itemizedlist><listitem>
// <p>This function is deprecated. </p>
// <p></p></listitem>
// <listitem>
// <p>Computes the nutation matrix and the rotation vectors (standard IERS 1996).</p> 
// </listitem>
// </itemizedlist>
//
// Parameters
// K   : [EPS0  ; DPSI  ; DEPS ] (3xN)
// KP  : [EPS0P ; DPSIP ; DEPSP] (3xN)
// KPP : [EPS0PP; DPSIPP; DEPSPP] (3xN)
// flag : (optional) "n" -> matrix only; "p" -> matrix and rotation vector; "s" -> matrix, rotation vector and its time derivative (default depends on number of outputs)
// NUT : Nutation 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_nutationAngles
// CL_mod_nutationArg
// CL_mod_meanObliquity
// CL_rot_angles2matrix
// CL_rot_angularVelocity
//
// Examples
// jj_tai = [19500:1:20500];
// ss_tai = 43200*ones(jj_tai);
// // Nutation angles (dPsi,dEps) :
// [NUT,F,NUTP,FP,NUTPP,FPP]=CL_mod_nutationAngles(jj_tai,ss_tai,"s");
//
// // Mean obliquity :
// [eps0,eps0p,eps0pp]=CL_mod_meanObliquity(jj_tai,ss_tai,"s");
// K=[eps0; NUT]; Kp=[eps0p; NUTP]; Kpp=[eps0pp; NUTPP];
//
// [NUT,OM,OMP] = CL_mod_nutationMatrix(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 nutation

// premiere  rotation EPS0 + DEPS axe 1
// deuxieme  rotation DPSI        axe 3
// troisieme rotation -EPS0       axe 1

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

//premiere rotation EPS0 + DEPS obliquite
tab(1,:) = K(1,:) + K(3,:);
tabp(1,:) = KP(1,:) + KP(3,:);
tabpp(1,:) = KPP(1,:) + KPP(3,:);

//deuxieme rotation DPSI longitude
tab(2,:) = K(2,:);
tabp(2,:) = KP(2,:);
tabpp(2,:) = KPP(2,:);

//troisieme rotation -EPS0 obliquite moyenne
tab(3,:) = -K(1,:);
tabp(3,:) = -KP(1,:);
tabpp(3,:) = -KPP(1,:);

NUT=CL_rot_angles2matrix([1,3,1],tab);
OM=[]; OMP=[];

if flag=="p" then
  OM = CL_rot_angularVelocity([1,3,1],tab,tabp,tabpp);
elseif flag=="s" then
  [OM,OMP] = CL_rot_angularVelocity([1,3,1],tab,tabp,tabpp);
end

endfunction