//  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 [psi, eps, psidot, epsdot] = CL__iers_nuta2000AR06(jd, cder)
// Nutation, IAU 2000A R06 model
//
// Calling Sequence
// [psi, eps, psidot, epsdot] = CL__iers_nuta2000AR06(jd [,cder])
//
// Description
// <itemizedlist><listitem>
// <p>Compute the longitude (psi) and obliquity (eps) components of the 
// nutation using IAU2000A R06 nutation model (series psi/eps).</p>
// <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)
// psi: Nutation angle in longitude [rad] (1xN)
// eps: Nutation angle in obliquity [rad] (1xN)
// psidot: Time derivative of nutation angle in longitude [rad/s](1xN)
// epsdot: Time derivative of nutation angle in obliquity [rad/s] (1xN)
//
// Authors
// CNES - DCT/SB
//
// Bibliography
// 1) Technical Note 36, IERS, 2010
// 2) ftp://tai.bipm.org/iers/conv2010/chapter5/ 
//


// 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 (lunisolar and planetary) (units: rad and rad/s)
[Fls, Flsdot] = CL__iers_fundArgsLS(jd, cder);
[Fpl, Fpldot] = CL__iers_fundArgsPL(jd, cder);

// Note: 
// series for psi/eps and psidot/epsdot: initially arcseconds and arcseconds/century, 
// converted to rad and rad/s
// Note: fund args derivatives converted to rad/century
// because unit of t = century


// 1) psi/psidot
fname = fullfile(CL_home(), "data", "frame", "iers_series_psi_2006_2000AR06.dat"); 
if (~isfile(fname))
  CL__error("File: " + fname + ": not found"); 
end
load(fname, "data"); 

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


// 2) eps/epsdot
fname = fullfile(CL_home(), "data", "frame", "iers_series_eps_2006_2000AR06.dat"); 
if (~isfile(fname))
  CL__error("File: " + fname + ": not found"); 
end
load(fname, "data"); 

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


endfunction