//  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 [gamb,phib,psib,epsa, gambdot,phibdot,psibdot,epsadot] = CL__iers_prec2006(jd, cder)
// Precession angles (IAU2006)
//
// Calling Sequence
// [gamb,phib,psib,epsa, [gambdot,phibdot,psibdot,epsadot]] = CL__iers_prec2006(jd [,cder])
//
// Description
// Compute the 4 angles of Fukushima-Williams precession:
// - epsa is the mean obliquity of date
// - gamb (gamma bar) is the GCRS right ascension of the intersection of the
//   ecliptic of date with the GCRS equator
// - phib (phi bar) is the obliquity of the ecliptic of date on the GCRS equator 
// - psib (psi bar) is the precession angle plus bias in longitude along
//   the ecliptic of date
//
// 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)
// gamb, phib, psib, epsa: precession angles [rad] (1xN)
// gambdot, phibdot, psibdot, epsadot: Time derivatives of precession angles [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) <= 4); cder = %f; end;

ARCSEC_TO_RAD = %pi / (180*3600);

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

// (eq 5.40, p65, paragraph 5.6.4)
gamb = (-0.052928 + (10.556378 + (0.4932044 + (-0.00031238 + (-0.000002788 + ...
       (0.0000000260) .* t) .* t) .* t) .* t) .* t) * ARCSEC_TO_RAD;
           
phib = (84381.412819 + (-46.811016 + (0.0511268 + (0.00053289 + (-0.000000440 + ...
       (-0.0000000176) .* t) .* t) .* t) .* t) .* t) * ARCSEC_TO_RAD;
           
psib = (-0.041775 + (  5038.481484 + (1.5584175 + (-0.00018522 + (-0.000026452 + ...
       (-0.0000000148) .* t) .* t) .* t) .* t) .* t) * ARCSEC_TO_RAD;

// Mean obliquity of the ecliptic, IAU 2006 precession model.
[epsa, epsadot] = CL__iers_obli2006(jd,cder);       
      
// ---------------
// Time derivatives
// ---------------
gambdot = [];
phibdot = [];
psibdot = [];
if (cder)
  gambdot = (10.556378 + (2* 0.4932044 + (3*-0.00031238 + (4*-0.000002788 + ...
            (5* 0.0000000260) .* t) .* t) .* t) .* t) * ARCSEC_TO_RAD / (36525*86400);
             
  phibdot = (-46.811016 + (2* 0.0511268 + (3* 0.00053289 + (4*-0.000000440 + ...
            (5*-0.0000000176) .* t) .* t) .* t) .* t) * ARCSEC_TO_RAD / (36525*86400);
             
  psibdot =  (5038.481484 + (2* 1.5584175 + (3*-0.00018522 + (4*-0.000026452 + ...
             (5*-0.0000000148) .* t) .* t) .* t) .* t) * ARCSEC_TO_RAD / (36525*86400);
end 

endfunction