//  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 [zetaa,za,thetaa, zetaadot,zadot,thetaadot] = CL__iers_prec1976(jd, cder)
// IAU 1976 precession model.
//
// Calling Sequence
// [zetaa,za,thetaa, [zetaadot,zadot,thetaadot]] = CL__iers_prec1976(jd [,cder])
//
// Description
// <itemizedlist><listitem>
// <p>IAU 1976 precession model.</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)
// zetaa,za,thetaa: IAU 1976 precession [rad] (1xN)
// zetaadot,zadot,thetaadot: Time derivatives of IAU 1976 precession [rad/s] (1xN)
//
// Authors
// CNES - DCT/SB
//
// Bibliography
// 1) Technical Note 21, IERS, 1996


// Declarations:


// Code:

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

ARCSEC_TO_RAD = %pi / (180*3600);

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

// p22, chapter 5
zetaa = (2306.2181 + (0.30188  + 0.017998 .* t) .* t) .* t * ARCSEC_TO_RAD;
za = (2306.2181 + (1.09468  + 0.018203 .* t) .* t) .* t * ARCSEC_TO_RAD;
thetaa = (2004.3109 + (-0.42665 - 0.041833 .* t) .* t) .* t * ARCSEC_TO_RAD;
          
   
// ---------------
// Time derivatives     
// ---------------
zetaadot = [];
zadot = [];
thetaadot = []; 
if (cder)
  zetaadot = (2306.2181 + (2*0.30188  + 3*0.017998 .* t) .* t) * ARCSEC_TO_RAD / (36525*86400);
  zadot = (2306.2181 + (2*1.09468  + 3*0.018203 .* t) .* t) * ARCSEC_TO_RAD / (36525*86400);
  thetaadot = (2004.3109 + (2*-0.42665 - 3*0.041833 .* t) .* t) * ARCSEC_TO_RAD / (36525*86400);   
end 


endfunction