//  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_cirs2teme(jd, args, comega)
// CIRS to TEME frame transformation matrix and angular velocity
//
// Calling Sequence
// [M, omega] = CL__fr_cirs2teme(jd ,args [,comega])
//
// Description
// <itemizedlist><listitem>
// <p>Computes the frame transformation matrix <b>M</b> from CIRS to TEME reference frame.</p>
// <p> By convention, <b>M</b> multiplied by coordinates relative to CIRS yields coordinates relative to Gamma50 Veis frame. </p>
// <p>Optionaly computes the angular velocity vector <b>omega</b> of TEME relative to CIRS, with coordinates relative to CIRS.
// 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>
// <listitem>
// <p>See <link linkend="Reference frames">Reference frames</link> for more details on the definition of reference frames.</p> 
// <p>args is a structure containing the field "ut1_tt" which is the difference between
// UT1 time scale and TT time scale in seconds.</p>
// <p>Note : The transformation is a single rotation around Z axis of (era+sprime)-gmst82</p>
// <p>sprime should be computed with a TT time scale and era and gmst82 with a UT1 time scale 
// but in practice, using TT for both leads to a 100 microarseconds error.</p>
// <p></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)
//
// Authors
// CNES - DCT/SB
//

// Declarations:


// Code:

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

if ~isstruct(args); CL__error("args should be a structure"); end;
if (~isfield(args,"ut1_tt")); CL__error("Field(s) missing in args"); end;


// Convert from TT to UT1
jd_ut1 = [jd(1,:) ; ...
          jd(2,:) + args.ut1_tt/86400];

// Note :
// Rigorously, the date used for the sidereal times (ERA and GMST82) should be UT1 dates.
// In practice, (ERA+ s')-GMST82  = -0.12''+46.1''t , with t = number of days (TT) from J2000
// Using TT instead of UT1, the error is about 60 seconds on the date, which gives:
// 60/(86400*365)*46.1'' = 0.000087 arcseconds = 87 microarcseconds
// So the error is small but not negligible.

[gmst82,gmst82dot] = CL__iers_gmst1982(jd_ut1, comega);
[era, eradot] = CL__iers_era2000(jd_ut1, comega);
[sprime, sprimedot] = CL__iers_sp2000(jd, comega);

M = CL_rot_angles2matrix(3, (era+sprime)-gmst82);

// Angular velocity vector
omega = [];
if (comega)
  N = size(jd,2);
  omega = zeros(3,N);
  omega(3,:) = (eradot+sprimedot)-gmst82dot;
end

endfunction