//  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 [x, y, s, xdot, ydot, sdot] = CL__iers_xys2006A_cla(jd, cder)
// X,Y coordinates of the Celestial Intermediate Pole and CIO locator s, using classical angles (IAU 2006 precession and IAU 2000AR06 nutation)
//
// Calling Sequence
// [x, y, s, xdot, ydot, sdot] = CL__iers_xys2006A_cla(jd [,cder])
//
// Description
// <itemizedlist><listitem>
// <p>X,Y coordinates of the Celestial Intermediate Pole and CIO locator s 
// (using the IAU 2006 precession and IAU 2000AR06 nutation models)</p>
// Note: cla : "classic" <=> using classical angles. 
// <p></p></listitem>
// </itemizedlist>
//
// Parameters
// jd: Two-part julian day (Time scale: TT) (2xN)
// cder: (boolean, optional) Option to compute derivatives. If cder is %f, the derivatives will be set to []. Default is %t. (1x1)
// x: X coordinate of CIP [rad] (1xN)
// y: Y coordinate of CIP [rad] (1xN)
// s: CIO locator [rad] (1xN)
// xdot: Time derivative of X coordinate of CIP [rad/s] (1xN)
// ydot: Time derivative of Y coordinate of CIP [rad/s] (1xN)
// sdot: Time derivative of CIO locator [rad/s] (1xN)
//
// Authors
// CNES - DCT/SB
//
// Bibliography
// 1) Technical Note 36, IERS, 2010

// Declarations:


// Code:

if (~exists("cder","local")); cder = %t; end;
if (argn(1) <= 3); cder = %f; end;

// Frame bias and precession (IAU2006)
[gamb, phib, psib, epsa, gambdot, phibdot, psibdot, epsadot] = CL__iers_prec2006(jd, cder);

// Nutation components (IAU2000AR06)
[dp, de, dpdot, dedot] = CL__iers_nuta2000AR06(jd, cder);

// Extract the X,Y coordinates.
// This is the literal expression for x=NPB(3,1,:) and y=NPB(3,2,:)
// NPB is the equinox based nutation-precession-bias matrix.
// PB = CL_rot_angles2matrix([3,1,3,1], [ gamb; phib ; -psib ; -epsa]);
// N = CL_rot_angles2matrix([1,3,1],[epsa ; -dp ; -(epsa+de)]);
// NPB = N*PB;
// which is the same as:
// NPB = CL_rot_angles2matrix([3,1,3,1], [ gamb; phib ; -(psib+dp) ; -(epsa+de)]);
// 
// direct computation for efficiency: 
cos1 = cos(gamb);           
sin1 = sin(gamb);

cos2 = cos(phib);           
sin2 = sin(phib);

cos3 = cos(-(psib+dp));     
sin3 = sin(-(psib+dp));

cos4 = cos(-(epsa+de));     
sin4 = sin(-(epsa+de));

x = sin4.*(sin3.*cos1 + cos3.*cos2.*sin1) + cos4.*sin2.*sin1;
y = sin4.*(sin3.*sin1 - cos3.*cos2.*cos1) - cos4.*sin2.*cos1;

// ----------------------------
// Time derivatives of x and y
// ----------------------------
xdot = [];
ydot = [];
if (cder)       
  dot1 = gambdot;        
  dot2 = phibdot; 
  dot3 = -(psibdot+dpdot);
  dot4 = -(epsadot+dedot);
  
  xdot = dot4.*cos4.*(sin3.*cos1 + cos3.*cos2.*sin1) + ...
         sin4.*(dot3.*cos3.*cos1 + sin3.*dot1.*(-sin1) + dot3.*(-sin3).*cos2.*sin1 + cos3.*dot2.*(-sin2).*sin1 + cos3.*cos2.*dot1.*cos1) + ...
         dot4.*(-sin4).*sin2.*sin1 + cos4.*dot2.*cos2.*sin1 + cos4.*sin2.*dot1.*cos1;
  ydot = dot4.*cos4.*(sin3.*sin1 - cos3.*cos2.*cos1) + ...
         sin4.*(dot3.*cos3.*sin1 + sin3.*dot1.*cos1 - dot3.*(-sin3).*cos2.*cos1 - cos3.*dot2.*(-sin2).*cos1 - cos3.*cos2.*dot1.*(-sin1)) - ...
         dot4.*(-sin4).*sin2.*cos1 - cos4.*dot2.*cos2.*cos1 - cos4.*sin2.*dot1.*(-sin1);
end

// CIO locator s
[s, sdot] = CL__iers_s2006(jd, x, y, xdot, ydot, cder);

endfunction