//  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 [N, omega] = CL__fr_mod2tod(jd, args, comega)
// MOD (Mean of Date) to TOD (True of Date) vector transformation
//
// Calling Sequence
// [N, omega] = CL__fr_mod2tod(jd, args [,comega])

//
// Description
// <itemizedlist><listitem>
// <p>Computes the frame transformation matrix <b>N</b> from MOD to TOD : nutation. </p>
// By convention, multiplying <b>N</b> by coordinates relative to MOD yields coordinates relative to TOD.</p> 
// <p>Optionaly computes the angular velocity vector <b>omega</b> of TOD relative to MOD, with coordinates relative to MOD.
// 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 structure containing the following fields:</p> 
// <p>The main field of the structure used in this function is "nutation_model". Available models are : </p>
// <p>- "2000AR06": IAU2000AR06 Model. This is the default value. 
//             In that case "dx06" and "dy06" (CIP offsets wrt IAU 2006/2000A expressed in radians) should also be provided.</p>
// <p>- "1980": IAU1980 Model. In that case "ddp80" and "dde80" (corrections to the nutation longitude angle expressed in radians)
//            should also be provided.</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)
// N: Transformation matrix (3x3xN)
// omega: (optional) Angular velocity vector [rad/s] (3xN)
//
// Authors
// CNES - DCT/SB
//
// Bibliography
// 1) Technical Note 36, IERS, 2010


// Declarations:


// Code:

// Internal function :
// Transforms CIP offsets (dx06,dy06) to nutation corrections wrt IAU 2006/2000A
function [ddp06,dde06] = CL__iers_cip2nutCorr(jd, dx06, dy06, epsa, dpsi, deps)
  if (dx06 <> 0 | dy06 <> 0)
    args = struct();
    args.precession_model = "2006";
    PB = CL__fr_gcrs2mod(jd, args);
    N = CL_rot_angles2matrix([1,3,1],[epsa ; -dpsi ; -(epsa+deps)]);
    V1 = [ dx06 ; dy06 ; zeros(dx06) ];
    V2 = N*PB * V1;
    ddp06 = V2(1,:) ./ sin(epsa);
    dde06 = V2(2,:);
  else
    ddp06 = 0;
    dde06 = 0;
  end
endfunction

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,"nutation_model"));  CL__error("Field(s) missing in args"); end;

// ---------------------------------------------------------
// Nutation, IAU2000AR06 
// ---------------------------------------------------------
if (args.nutation_model == "2000AR06")
  if (~isfield(args,"dx06") | ~isfield(args,"dy06")) 
    CL__error("Field(s) missing in args");
  end

  [dpsi, deps, dpsidot, depsdot] = CL__iers_nuta2000AR06(jd, comega);
  [epsa, epsadot] = CL__iers_obli2006(jd, comega);
  
  // Apply nutation corrections (obtained from CIP offsets : dx06,dy06)
  [ddp06,dde06] = CL__iers_cip2nutCorr(jd, args.dx06, args.dy06, epsa, dpsi, deps);
  dpsi = dpsi + ddp06;
  deps = deps + dde06;
  
  
  
// ---------------------------------------------------------
// Nutation, IAU1980
// ---------------------------------------------------------
elseif (args.nutation_model == "1980")
  if (~isfield(args,"ddp80") | ~isfield(args,"dde80")) 
    CL__error("Field(s) missing in args");
  end
  
  [dpsi, deps, dpsidot, depsdot] = CL__iers_nuta1980(jd, comega);
  [epsa, epsadot] = CL__iers_obli1980(jd, comega);
  
  // Apply nutation corrections
  dpsi = dpsi + args.ddp80;
  deps = deps + args.dde80;
  
else
  CL__error("Unkown nutation model");
end


N = CL_rot_angles2matrix([1,3,1],[epsa ; -dpsi ; -(epsa+deps)]);

// Angular velocity
omega = [];
if (comega) 
  omega = CL_rot_angVelocity([1,3,1],[epsa ; -dpsi ; -(epsa+deps)],[epsadot ; -dpsidot ; -(epsadot+depsdot)]);
end


endfunction