//  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 [pos_G50,vel_G50,jacob] = CL_fr_ter2G50(cjd,pos_ter, vel_ter,ut1_utc)
// Terrestrial frame to Gamma50 (Veis) vector transformation - DEPRECATED
//
// Calling Sequence
// [pos_G50,vel_G50,jacob] = CL_fr_ter2G50(cjd,pos_ter [,vel_ter,ut1_utc])
//
// Description
// <itemizedlist><listitem>
// <p>This function is deprecated. </p>
// <p>Replacement function: <link linkend="CL_fr_convert">CL_fr_convert</link></p>
// <p></p></listitem>
// <listitem>
// <p>Converts position and (optionally) velocity vectors relative to the terrestrial ("Earth fixed") reference frame to 
// position and (optionally) velocity vectors relative to Gamma50 (Veis). </p>
// <p>The jacobian of the transformation is optionally computed.</p>
// <p></p></listitem>
// <listitem>
// <p><b>Notes:</b></p>
// <p> Transformation from Gamma50 (Veis) to terrestrial frame consists in a single rotation 
// of the Veis sideral time (see <link linkend="CL_mod_sidTimeG50">CL_mod_sidTimeG50</link>) around the Z axis.</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 day from 1950.0 (UTC) (1xN or 1x1)
// pos_ter: position vector relative to terrestrial frame [m] (3xN or 3x1)
// vel_ter: (optional) velocity vector relative to terrestrial frame [m/s] (3xN or 3x1)
// ut1_utc : (optional) ut1-utc [seconds] (default is 0) (1xN or 1x1)
// pos_G50: position vector relative to Gamma50 (Veis) [m] (3xN)
// vel_G50: (optional) velocity vector relative to Gamma50 (Veis) [m/s] (3xN)
// jacob: (optional) jacobian of the transformation [d(x,y,z,vx,vy,vz)_G50/d(x,y,z,vx,vy,vz)_ter] (6x6xN)
//
// Authors
// CNES - DCT/SB
//
// Bibliography
// 1) CNES - MSLIB FORTRAN 90, Volume R (mr_TerVrai_veis)
//
// See also
// CL_fr_ter2J2000Mat
// CL_fr_G502terMat
// CL_fr_G502ter
// CL_mod_sidTimeG50
//
// Examples
// // Conversion of position ter to G50
// pos_ter = [[3500.e3;2500.e3;5800.e3] , [4500.e3;2100.e3;6800.e3]];
// cjd = [21010 , 21011];
// pos_G50=CL_fr_ter2G50(cjd,pos_ter);
//
// // Conversion of position and velocity : ter to G50  + jacobian
// pos_ter = [[3500.e3;2500.e3;5800.e3] , [4500.e3;2100.e3;6800.e3]];
// vel_ter = [[1.e3;3.e3;7.e3] , [2.e3;3.e3;6.e3]];
// cjd = [21010 , 21011];
// [pos_G50,vel_G50,jacob]=CL_fr_ter2G50(cjd,pos_ter,vel_ter);


// Declarations:


// Code:

CL__warnDeprecated(); // deprecated function

vel_present = %t;
if (~exists('vel_ter','local')); vel_present = %f; end
if (~exists('ut1_utc','local')); ut1_utc=zeros(cjd); end

compute_vel = %f;
compute_jacob = %f;
vel_G50 = [];
jacob = [];

[lhs,rhs] = argn()
if (rhs<2 | rhs > 4) then CL__error('check number of input arguments'); end


if (lhs >=2)
  compute_vel = %t;
end
if (lhs ==3)
  compute_jacob = %t;
end

if (~vel_present & compute_vel)
  CL__error("Input velocity missing");
end

// Check argument sizes, and resize if necessary: 
if (vel_present)
  [cjd,pos_ter,vel_ter,ut1_utc] = CL__checkInputs(cjd,1, pos_ter,3, vel_ter,3, ut1_utc,1);
else
  [cjd,pos_ter,ut1_utc] = CL__checkInputs(cjd,1, pos_ter,3, ut1_utc,1);
end

// Position computation :
cjd_ut1 = cjd + ut1_utc/86400.;
[tsid,tsidt] = CL_mod_sidTimeG50(cjd_ut1);

costsid = cos(tsid);
sintsid = sin(tsid);

pos_G50 = [ costsid .* pos_ter(1,:) - sintsid .* pos_ter(2,:);
            sintsid .* pos_ter(1,:) + costsid .* pos_ter(2,:);
            pos_ter(3,:) ];


// Velocity/jacobian computation
if (compute_vel)
  n = size(pos_ter,2)

  // omega G50/ter (coordinates in either frame)
  omega = [ zeros(1,n); zeros(1,n); -tsidt.*ones(1,n) ];

  // next line: vel_G50 = vel/ter, coordinates in G50
  vel_G50 = [ costsid .* vel_ter(1,:) - sintsid .* vel_ter(2,:);
              sintsid .* vel_ter(1,:) + costsid .* vel_ter(2,:);
              vel_ter(3,:) ];

  // next line: vel_G50 = vel/G50, coordinates in G50
  vel_G50 = vel_G50 - CL_cross(omega, pos_G50); 


  // jacobian computation
  if (compute_jacob)
 
    // transf matrix ter -> G50 
    // X_G50 = mat * X_ter 
    mat = zeros(3,3,n); 
    mat(1,1,1:n) = costsid; 
    mat(1,2,1:n) = -sintsid; 
    mat(2,1,1:n) = sintsid; 
    mat(2,2,1:n) = costsid; 
    mat(3,3,1:n) = 1; 

    jacob = CL_rot_pvJacobian(mat, omega); 

  end

end


endfunction