//  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_cirs2tod(jd, args, comega)
// CIRS to TOD transformation matrix and angular velocity vector
//
// Calling Sequence
// [M, omega] = CL__fr_cirs2tod(jd [,args, comega])
//
// Description
// <itemizedlist><listitem>
// <p>Computes the frame transformation matrix <b>M</b> from CIRS to TOD</p>
// <p>By convention, multiplying <b>M</b> by coordinates relative to CIRS yields coordinates relative to TOD.</p> 
// <p>Optionaly computes the angular velocity vector <b>omega</b> of TOD 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>
// <p>The transformation consists of a transformation to GCRS, 
//    then to MOD and to TOD using IAU2006 precession and IAU2006A nutation</p>
// <p>Note : args appears in the calling sequence but is not used.</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: Not used.
// 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
//
// Bibliography
// 1) Technical Note 36, IERS, 2010

// Declarations:

// Code:

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


args = struct();
// X,Y,s obtained by IAU2006 precession and IAU2000AR06 nutation 
// (interpolation or using classical angles if outside [2000,2100])
args.model = "interp";     
args.dx06 = 0;
args.dy06 = 0;
args.precession_model = "2006";    // IAU2006 precession
args.nutation_model = "2000AR06";  // IAU2000AR06 nutation


// om1 = omega(CIRS/GCRS), coordinates in GCRS
[M1,om1] = CL__fr_gcrs2cirs(jd, args, comega);
// om2 = omega(MOD/GCRS), coordinates in GCRS
[M2,om2] = CL__fr_gcrs2mod(jd, args, comega);
// om3 = omega(TOD/MOD), coordinates in MOD
[M3,om3] = CL__fr_mod2tod(jd, args, comega);
  
// CIRS-->GCRS-->MOD-->TOD
omega = [];
if (comega)
  [M,omega] = CL_rot_compose(M1,om1,-1, M2,om2,1, M3,om3,1);
else
  // CIRS-->GCRS-->MOD-->TOD
  M = M3*M2*M1';
end


endfunction