//  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] = CL_fr_ter2J2000Mat(cjd, ut1_utc,xp,yp,dPsi,dEps,conv_iers)
// Terrestrial frame to EME2000 (J2000) transformation matrix - DEPRECATED
//
// Calling Sequence
// M = CL_fr_ter2J2000Mat(cjd [,ut1_utc,xp,yp,dPsi,dEps,conv_iers])
//
// Description
// <itemizedlist><listitem>
// <p>This function is deprecated. </p>
// <p>Replacement function: <link linkend="CL_fr_convertMat">CL_fr_convertMat</link></p>
// <p></p></listitem>
// <listitem>
// <p>Computes the frame transformation matrix <b>M</b> from the 
// terrestrial ("Earth fixed") reference frame to EME2000.</p>
// <p> By convention, multiplying <b>M</b> by coordinates relative to the terrestrial frame yields coordinates relative to EME2000. </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
// cjd: modified julian days from 1950.0 (UTC) (1xN)
// ut1_utc : (optional) ut1-utc [seconds] (default is 0) (1xN)
// xp : (optional) x polar coordinate [radians] (default is 0) (1xN)
// yp : (optional) y polar coordinate [radians] (default is 0) (1xN)
// dPsi : (optional) Nutation corrections [radians] (default is 0) (1xN)
// dEps : (optional) Nutation corrections [radians] (default is 0) (1xN)
// conv_iers : (optional) Convention IERS. Only iers 1996 (Lieske/Wahr) is implemented (default is "iers_1996")
// M: terrestrial frame to EME2000 transformation matrix (3x3xN)
//
// Authors
// CNES - DCT/SB
//
// Bibliography
// 1) IERS Conventions (1996), Dennis D. McCarthy
// 2) Explanatory Supplement to the Astronomical Almanac, Seidelman (1992)
//
// See also
// CL_mod_precessionMatrix
// CL_mod_nutationMatrix
// CL_mod_polarMotionMatrix
//
// Examples
// // Conversion terrestrial to J2000
// pos_ter = [ 7000.e3; 0; 0];
// cjd = [21010, 21011];
// M = CL_fr_ter2J2000Mat(cjd);
// pos_J2000 = M * pos_ter;


// Declarations:


// Code:

CL__warnDeprecated(); // deprecated function

[lhs,rhs]=argn();

if (~exists('ut1_utc','local')); ut1_utc = zeros(cjd); end
if (~exists('xp','local')); xp = zeros(cjd); end
if (~exists('yp','local')); yp = zeros(cjd); end
if (~exists('dPsi','local')); dPsi = zeros(cjd); end
if (~exists('dEps','local')); dEps = zeros(cjd); end
if (~exists('conv_iers','local')); conv_iers = "iers_1996"; end

if conv_iers~="iers_1996" then
      CL__error("only model iers_1996 is implemented, sorry...");
end

//verif dimension des args en entree
sz=size(cjd);
if (sz(1)~=1) then CL__error("check dim of input arguments (column vect.)"); end
sz=size(ut1_utc);
if (sz(1)~=1) then CL__error("check dim of input arguments (column vect.)"); end
sz=size(xp);
if (sz(1)~=1) then CL__error("check dim of input arguments (column vect.)"); end
sz=size(yp);
if (sz(1)~=1) then CL__error("check dim of input arguments (column vect.)"); end
sz=size(dPsi);
if (sz(1)~=1) then CL__error("check dim of input arguments (column vect.)"); end
sz=size(dEps);
if (sz(1)~=1) then CL__error("check dim of input arguments (column vect.)"); end


W=CL_mod_polarMotionMatrix(xp,yp);

//Conversion TUC --> TAI
[cjdtai]=CL_dat_utc2tai(cjd);
jj_tai = floor(cjdtai);
sec_tai = (cjdtai-jj_tai)*86400.0;

cjdut1 = cjd + ut1_utc/86400.0;

TempsSideralMoyen = CL_mod_sidTime(cjdut1);

eps0 = CL_mod_meanObliquity(jj_tai,sec_tai,"n");

[NUT,F] = CL_mod_nutationAngles(jj_tai,sec_tai,"p");

[equi,equip] = CL_mod_equinoxesEquation(NUT(1,:),eps0,F(5,:));
TempsSideralVrai = TempsSideralMoyen+equi;
R = CL_rot_angles2matrix(3,-TempsSideralVrai);

NUT=NUT+[dPsi;dEps];
[N,OM,OMP] = CL_mod_nutationMatrix([eps0; NUT;],[],[],"n");

[K,KP,KPP]=CL_mod_precessionAngles(jj_tai,sec_tai,"n");
[PREC,OM,OMP] = CL_mod_precessionMatrix(K,KP,KPP,"n");

M=PREC*N*R*W;

endfunction