//  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 [M, omega] = CL__fr_tirs2itrs(jd, args, comega)
// TIRS to ITRS frame transformation matrix and angular velocity
//
// Calling Sequence
// [M, omega] = CL__fr_tirs2itrs(jd ,args [,comega])
//
// Description
// <itemizedlist><listitem>
// <p>Computes the frame transformation matrix <b>M</b> from TIRS to ITRS.
// By convention, multiplying <b>M</b> by coordinates relative to TIRS yields coordinates relative to ITRS.</p> 
// <p>Optionaly computes the angular velocity vector <b>omega</b> of ITRS relative to TIRS, with coordinates relative to TIRS.
// See <link linkend="Data types">Data types</link> for more details on the definition of angular velocity vectors and frame transformation matrix.</p>
// <p></p></listitem>
// <p>args is a structure containing the fields "xp" and "yp" which are the polar motion. (radians)</p>
// <p></p></listitem>
// <listitem>
// <p>See <link linkend="Reference frames">Reference frames</link> for more details on the definition of reference frames.</p> 
// </listitem>
// </itemizedlist>
//
// Parameters
// jd: Two-part julian day (Time scale: TT) (2xN)
// args: Structure of arguments. See description for more details.
// comega: (boolean, optional) Option to compute omega. If comega is %f, omega will be set to []. Default is %t. (1x1)
// M: Transformation matrix (3x3xN)
// omega: (optional) Angular velocity vector [rad/s] (3xN or 3x1)
//
// Authors
// CNES - DCT/SB
//
// Bibliography
// 1) Technical Note 36, IERS, 2010
//
// See also


// Declarations:


// Code:

if (~exists("comega","local")); comega = %t; end;
if (argn(1) <= 1); comega = %f; end;

// Resize xp,yp
[jd, args.xp, args.yp, N] = CL__checkInputs(jd,2, args.xp,1, args.yp,1);


[sp, spdot] = CL__iers_sp2000(jd, comega);
M = CL_rot_angles2matrix([3,2,1], [sp; -args.xp; -args.yp]);

// Angular velocity vector
omega = [];
if (comega)
  // NB: xpdot, ypdot = 0
  omega = [zeros(2,N); spdot];
end

endfunction