//  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__iers_xys2Mat(x, y, s, xdot, ydot, sdot, comega)
// GCRS to CIRS frame transformation matrix given the CIP coordinates X,Y and the CIO locator s.
//
// Calling Sequence
// [M,omega] = CL__iers_xys2Mat(x, y, s, [xdot, ydot, sdot, comega])
//
// Description
// <itemizedlist><listitem>
// <p>Computes the frame transformation matrix <b>M</b> from GCRS to CIRS</p>
// <p>By convention, multiplying <b>M</b> by coordinates relative to GCRS yields coordinates 
// relative to CIRS.</p> 
// <p>The function also optionnaly computes the angular velocity vector <b>omega</b></p> 
// <p></p></listitem>
// </itemizedlist>
//
// Parameters
// x,y: Coordinates of the Celestial Intermediate Pole (CIP) [unitless] (1xN)
// s: CIO locator [unitless] (1xN)
// comega: (boolean, optional) Option to compute omega. If comega is %f, omega will be set to []. Default is %t. (1x1)
// xdot,ydot: (optional) Time derivatives of coordinates of the CIP [unitless/s] (1xN)
// sdot: (optional) Time derivative of CIO locator [unitless/s] (1xN)
// M: GCRS to CIRS frame 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;

// Compute the spherical angles E and d
// Use eq 5.7, p48, paragraph 5.4.4 : X = sind*cosE, Y = sind*sinE, Z = cosd
// E <=> "right ascension"
// d <=> "%pi/2 - declination"
// rather than eq 5.10: Q(t) = mat * R3(s) 
// in order to simplify the computation of the derivatives
sind2 = x.*x + y.*y; // = (sind)^2

if (find(sind2 > 1))
  CL__error("Invalid values for x,y");
end

E = atan(y, x); // note if x and y == 0 : E = 0 (impossible in practice)
d = asin(sqrt(sind2)); // d >= 0 because z is assumed positive

// Form the matrix (eq 5.6, p48, paragraph 5.4.4)
// (M = Q' [Q: see IERS doc])
M = CL_rot_angles2matrix( [3,2,3] , [E; d ; -(E+s)]);

// Angular velocity vector
omega = [];
if (comega)
  Edot = (ydot .* x - y .* xdot) ./ sind2;
  ddot = (x.*xdot + y.*ydot) ./ sqrt(sind2 .* (1.0 - sind2));
  omega = CL_rot_angVelocity( [3,2,3] , [E; d ; -(E+s)], [Edot; ddot ; -(Edot+sdot)]);
end

endfunction