//  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_gcrs2mod(jd, args, comega)
// GCRS to Mean Of Date transformation matrix and angular velocity vector
//
// Calling Sequence
// [M, omega] = CL__fr_gcrs2mod(jd, args [,comega])
//
// Description
// <itemizedlist><listitem>
// <p>Computes the frame transformation matrix <b>M</b> from GCRS to MOD : Precession + frame bias</p>
// <p>By convention, multiplying <b>M</b> by coordinates relative to GCRS yields coordinates relative to MOD.</p> 
// <p>Optionaly computes the angular velocity vector <b>omega</b> of MOD relative to GCRS, with coordinates relative to GCRS.
// 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 one field: "precession_model"</p> 
// <p>Available precession models are :</p>
// <p>- "1976" which uses Lieske angles for precession. (IAU1976)</p>
// <p>- "2006" which uses Fukushima-Williams angles for precession. (IAU2006)</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)
//
// 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;

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


// ---------------------------------------------------------
// IAU2006 : Fukushima-Williams 2006 precession (includes bias)
// ---------------------------------------------------------
if (args.precession_model == "2006")
  
  [gamb,phib,psib,epsa, gambdot,phibdot,psibdot,epsadot] = CL__iers_prec2006(jd, comega);
  
  M = CL_rot_angles2matrix([3,1,3,1], [gamb; phib; -psib; -epsa]);
    
  // Angular velocity vector
  omega = [];
  if (comega) 
    omega = CL_rot_angVelocity([3,1,3,1],[gamb;phib;-psib;-epsa], [gambdot;phibdot;-psibdot;-epsadot]);
  end
  
  
// ---------------------------------------------------------
// IAU1976 : Lieske precession
// ---------------------------------------------------------
elseif (args.precession_model == "1976")
  
  B = CL__fr_gcrs2eme2000(jd);
  
  [zetaa,za,thetaa, zetaadot,zadot,thetaadot] = CL__iers_prec1976(jd, comega);
  
  P = CL_rot_angles2matrix([3,2,3],[-zetaa;thetaa;-za]);
  M = P*B;
  
  // Angular velocity vector
  omega = [];
  if (comega)
    // omega = omega(MOD/EME2000) with coordinates in EME2000 
    //       = omega(MOD/GCRS) with coordinates in EME2000 because omega (GCRS/EME2000)=0
    omega = CL_rot_angVelocity([3,2,3],[-zetaa;thetaa;-za], [-zetaadot;thetaadot;-zadot]);
    // We want omega = omega(MOD/GCRS) with coordinates in GCRS
    omega = B'*omega;
  end
  
else
  CL__error("Unkown precession model");
end
  
endfunction