//  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, xdot, ydot] = CL__iers_xy2006A_ser(jd, cder)
// X,Y coordinates of celestial intermediate pole from series (IAU 2006 precession and IAU 2000AR06 nutation).
//
// Calling Sequence
// [x, y, xdot, ydot] = CL__iers_xy2006A_ser(jd [,cder])
//
// Description
// <itemizedlist><listitem>
// <p>Compute the X,Y coordinates of CIP from series based on IAU 2006 precession 
// and IAU 2000A nutation</p>
// The X,Y coordinates are those of the unit vector of the
// celestial intermediate pole.  They represent the combined effects
// of frame bias, precession and nutation.
// <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)
// xdot: Time derivative of X coordinate of CIP [rad/s] (1xN)
// ydot: Time derivative of Y coordinate of CIP [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) <= 2); cder = %f; end;

// Conversion factors
ARCSEC_TO_RAD = %pi / (180*3600);
CENT_TO_SEC = 86400 * 36525; 

// Julian centuries since J2000 (12h00)
t = ((jd(1,:) - 2451545) + jd(2,:)) / 36525;

// Fundamental arguments (units: rad and rad/s)
[Fls, Flsdot] = CL__iers_fundArgsLS(jd, cder);
[Fpl, Fpldot] = CL__iers_fundArgsPL(jd, cder);


// loads series and compute X
fname = fullfile(CL_home(), "data", "frame", "iers_series_x_2006_2000AR06.dat"); 
if (~isfile(fname))
  CL__error("File: " + fname + ": not found"); 
end
load(fname, "data"); 

// x and xdot: initially arcseconds and arcseconds/century, 
// converted to rad and rad/s
// Note: fund args derivatives converted to rad/century
// because unit of t = century

[x, xdot] = CL__iers_seriesEval(data, t, [Fls; Fpl], [Flsdot; Fpldot]*CENT_TO_SEC, cder); 
x = x * ARCSEC_TO_RAD; 
xdot = xdot * ARCSEC_TO_RAD / CENT_TO_SEC; 


// loads series and compute X
fname = fullfile(CL_home(), "data", "frame", "iers_series_y_2006_2000AR06.dat"); 
if (~isfile(fname))
  CL__error("File: " + fname + ": not found"); 
end
load(fname, "data"); 

// y and ydot: initially arcseconds and arcseconds/century, 
// converted to rad and rad/s
// Note: fund args derivatives converted to rad/century
// because unit of t = century

[y, ydot] = CL__iers_seriesEval(data, t, [Fls; Fpl], [Flsdot; Fpldot]*CENT_TO_SEC, cder); 
y = y * ARCSEC_TO_RAD; 
ydot = ydot * ARCSEC_TO_RAD / CENT_TO_SEC; 


endfunction