//  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_bodyConvertMat(body,frame1,frame2,cjd, tt_tref)
// Transformation matrix and angular velocity from a body frame to another
//
// Calling Sequence
// [M,omega] = CL_fr_bodyConvertMat(body,frame1,frame2,cjd [,tt_tref])
//
// Description
// <itemizedlist><listitem>
// <p>Computes the frame transformation from <b>frame1</b> to <b>frame2</b>, which is characterized by 
// a frame transformation matrix <b>M</b> and an angular velocity vector <b>omega</b>. </p>
// <p>By convention: </p>
// <p>- multiplying <b>M</b> by coordinates relative to <b>frame1</b> yields coordinates relative to <b>frame2</b>. </p>
// <p>- <b>omega</b> is the angular velocity vector of <b>frame2</b> wrt <b>frame1</b>, with coordinates in <b>frame1</b>. </p>
// <p></p></listitem>
// <listitem>
// <p>Available frames are : </p>
// <p>    'ICRS': International Celestial Reference System</p>
// <p>    'BCI': Body Centered Inertial</p>
// <p>    'BCF': Body Centered body Fixed</p>
// <p></p>
// <p>See <link linkend="Reference frames">Reference frames</link> for more details on the definition of reference frames.</p> 
// <p></p>
// <p>Available bodies are: "Mercury","Venus","Mars","Jupiter","Saturn","Uranus", "Neptune", "Sun" and "Moon" </p>
// <p></p></listitem>
// <listitem>
// <p>Notes : </p>
// <p>- The date "cjd" is relative to the TREF time scale. </p>
// <p>- The frame names and body names are case sensitive</p>
// <p>- Earth is not a valid body for this function. See CL_fr_convert.</p>
// <p>- GCRS (identical to ICRS in CelestLab) is also accepted as a frame name.</p>
// </listitem>
// </itemizedlist>
//
// Parameters
// body: (string) Name of the body. ("Mercury","Venus","Mars","Jupiter","Saturn","Uranus", "Neptune", "Sun" or "Moon") (1x1)
// frame1: (string) Name of the initial frame. (1x1)
// frame2: (string) Name of the final frame. (1x1)
// cjd: Modified (1950.0) julian day (Time scale: TREF) (1xN)
// tt_tref: (optional) TT-TREF [seconds]. Default is %CL_TT_TREF. (1xN or 1x1)
// M: Frame transformation matrix from frame1 to frame2 (3x3xN)
// omega: Angular velocity vector of frame2 wrt frame1 with coordinates in frame1 [rad/s] (3xN)
//
// Authors
// CNES - DCT/SB
//
// Bibliography
// Report of the IAU/IAG working group on cartographic coordinates and rotational elements: 2009
//
// See also
// CL_mod_IAUBodyAngles
// CL_fr_bodyConvert
// CL_fr_convertMat
//
// Examples
// cjd = CL_dat_cal2cjd(2010,02,03,05,35,25);
//
// // ICRS to BCF
// [M,omega] = CL_fr_bodyConvertMat("Mars","ICRS","BCI",cjd);


// Declarations:


// Code:
if (argn(2) <> 4 & argn(2) <> 5)
  CL__error("Invalid number of input arguments"); 
end

if (~exists("tt_tref", "local")); tt_tref = CL__dataGetEnv("TT_TREF"); end
if (typeof(body) <> "string"); CL__error("Invalid input argument ''body''"); end;
if (typeof(frame1) <> "string"); CL__error("Invalid input argument ''frame1''"); end;
if (typeof(frame2) <> "string"); CL__error("Invalid input argument ''frame2''"); end;

// Substitute GCRS by ICRS
if (frame1 == "GCRS") frame1 = "ICRS"; end;
if (frame2 == "GCRS") frame2 = "ICRS"; end;

// Validity of frame name
frame_names = ["ICRS", "BCI", "BCF"];
if (find(frame1 == frame_names) == [] | find(frame2 == frame_names) == [])
  CL__error("Invalid frame name (case sensitive!)");
end

// Validity of body:
body_names = ["Mercury", "Venus", "Mars", "Jupiter", "Saturn", "Uranus", "Neptune", "Sun", "Moon"];
if (find(body == body_names) == [])
  CL__error("Invalid body name (case sensitive!)");
end

// Check/resize input arguments
[cjd, tt_tref, N] = CL__checkInputs(cjd,1, tt_tref,1);

// ang = [RA;DEC;W]
// angdot = [RAdot;DECdot;Wdot]
[ang,angdot] = CL_mod_IAUBodyAngles(cjd, body, tt_tref);

// ICRS-->BCI
M1 = CL_rot_angles2matrix([3,1], [ang(1,:)+%pi/2; %pi/2-ang(2,:)]);
omega1 = CL_rot_angVelocity([3,1], [ang(1,:)+%pi/2; %pi/2-ang(2,:)], [angdot(1,:); -angdot(2,:)]);
// BCI-->BCF                                
M2 = CL_rot_angles2matrix(3, ang(3,:));
omega2 = CL_rot_angVelocity(3, ang(3,:), angdot(3,:));

if (frame1 == "ICRS" & frame2 == "BCI")
  M = M1;
  omega = omega1;
  
elseif (frame1 == "ICRS" & frame2 == "BCF")
  [M,omega] = CL_rot_compose(M1,omega1,1, M2,omega2,1);
  
elseif (frame1 == "BCI" & frame2 == "BCF")
  M = M2;
  omega = omega2;
  
elseif (frame1 == "BCI" & frame2 == "ICRS")
  [M,omega] = CL_rot_compose(M1,omega1,-1);
  
elseif (frame1 == "BCF" & frame2 == "ICRS")
  [M,omega] = CL_rot_compose(M2,omega2,-1, M1,omega1,-1);
  
elseif (frame1 == "BCF" & frame2 == "BCI")
  [M,omega] = CL_rot_compose(M2,omega2,-1);
  
elseif (frame1 == frame2)
  M = repmat(eye(3,3),1,1,N);
  omega = zeros(3,N);
end

endfunction