//  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 [dpsi, deps, dpsidot, depsdot] = CL__iers_nuta1980(jd, cder)
// Nutation angles (IAU 1980)
//
// Calling Sequence
// [dpsi, deps, dpsidot, depsdot] = CL__iers_nuta1980(jd [,cder])
//
// Description
// <itemizedlist><listitem>
// <p>Compute the longitude (dpsi) and obliquity (deps) of the nutation using
// IAU1980 nutation 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)
// dpsi: Longitude nutation angle [rad] (1xN)
// deps: Obliquity nutation angle [rad] (1xN)
// dpsidot: Time derivative longitude nutation angle [rad/s] (1xN)
// depsdot: Time derivative obliquity nutation angle [rad/s] (1xN)
//
// Authors
// CNES - DCT/SB
//
// Bibliography
// 1) 1980 IAU theory of nutation: The final report of 
//    the IAU working group on nutation, P.K.Seidelmann, 1981
//


// Declarations:


// Code:

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


ARCSEC_TO_RAD = %pi / (180*3600);

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

// ------------------
// LUNI-SOLAR NUTATION
// ------------------

// Luni-Solar nutation series
[LS,sp,spt,ce,cet] = CL__iers_nuta1980_LS()

// --------------------- 
// Fundamental arguments (ref1, p98, paragraph 10)
// --------------------- 

// Mean longitude of Moon minus mean longitude of Moon's perigee. 
// (original equation: first term = 134deg 57' 46.733", second term = 198deg 52' 02.633")
el = (485866.733 + (715922.633 + (31.310 + 0.064 .* t) .* t) .* t) * ARCSEC_TO_RAD + CL_rMod(1325.0 * t, 1.0) * 2*%pi;

// Mean longitude of Sun minus mean longitude of Sun's perigee.
// (original equation: first term = 357deg 31' 39.804", second term = 359deg 03' 01.224")
elp = (1287099.804 + (1292581.224 + (-0.577 - 0.012 .* t) .* t) .* t) * ARCSEC_TO_RAD + CL_rMod(99.0 * t, 1.0) * 2*%pi;

// Mean longitude of Moon minus mean longitude of Moon's node.
// (original equation: first term = 93deg 16' 18.877", second term = 82deg 01' 03.137")
f = (335778.877 + (295263.137 + (-13.257 + 0.011 .* t) .* t) .* t) * ARCSEC_TO_RAD + CL_rMod(1342.0 * t, 1.0) * 2*%pi;

// Mean elongation of Moon from Sun.
// (original equation: first term = 297deg 51' 01.307", second term = 307deg 06' 41.328")
d = (1072261.307 + (1105601.328 + (-6.891 + 0.019 .* t) .* t) .* t) * ARCSEC_TO_RAD + CL_rMod(1236.0 * t, 1.0) * 2*%pi;

// Longitude of the mean ascending node of the lunar orbit on the 
// ecliptic, measured from the mean equinox of date.
// (original equation: first term = 125deg 02' 40.280", second term = 134deg 08' 10.539")
om = (450160.280 + (-482890.539 + (7.455 + 0.008 .* t) .* t) .* t) * ARCSEC_TO_RAD + CL_rMod(-5.0 * t, 1.0) * 2*%pi;

Fls = CL_rMod([el; elp; f; d; om], 2*%pi);

// arg = linear combination of fundamental arguments
arg = LS * Fls ;
sarg = sin(arg);
carg = cos(arg);

// Summation of luni-solar nutation series
dpsi = sum( (sp*ones(t) + spt * t) .* sarg , "r");
deps = sum( (ce*ones(t) + cet * t) .* carg , "r");

// Convert from 0.1 milliarcsec units to radians.
dpsi = dpsi * ARCSEC_TO_RAD * 1.e-4;
deps = deps * ARCSEC_TO_RAD * 1.e-4;



// ----------------
// Time derivatives
// ----------------
dpsidot = [];
depsdot = [];
if (cder)
  // rad/century
  eldot = (715922.633 + (2*31.310 + 3*0.064 .* t) .* t) * ARCSEC_TO_RAD + 1325.0 * 2*%pi;
  elpdot = (1292581.224 + (2*-0.577 - 3*0.012 .* t) .* t) * ARCSEC_TO_RAD + 99.0 * 2*%pi;
  fdot = (295263.137 + (2*-13.257 + 3*0.011 .* t) .* t) * ARCSEC_TO_RAD + 1342.0 * 2*%pi;
  ddot = (1105601.328 + (2*-6.891 + 3*0.019 .* t) .* t) * ARCSEC_TO_RAD + 1236.0 * 2*%pi;
  omdot = (-482890.539 + (2*7.455 + 3*0.008 .* t) .* t) * ARCSEC_TO_RAD + -5.0 * 2*%pi;

  argdot = LS * [eldot; elpdot; fdot; ddot; omdot];
  sargdot = argdot .* carg;
  cargdot = -argdot .* sarg;

  dpsidot = sum( (sp*ones(t) + spt * t) .* sargdot + (spt*ones(t)) .* sarg , "r");
  depsdot = sum( (ce*ones(t) + cet * t) .* cargdot + (cet*ones(t)) .* carg , "r");

   // Convert from 0.1 milliarcsec/century to radians/sec
  dpsidot = dpsidot * ARCSEC_TO_RAD * 1.e-4 / (36525*86400);
  depsdot = depsdot * ARCSEC_TO_RAD * 1.e-4 / (36525*86400);
end


endfunction


// Luni-solar nutation series (ref1, Table 1, p81-82-83, paragraph 2)
// LS = coefficient of l, lp, ... (Px5)
function [LS,sp,spt,ce,cet] = CL__iers_nuta1980_LS()

  // Description of columns : 
  // coef for: l lp F D omega + values of: sp spt ce cet
  // sp = coeff for sin(psi)
  // spt = coeff for t*sin(psi)
  // ce = coeff for cos(eps)
  // cet = coeff for t*cos(eps)
  // NB: 
  // psi = longitude
  // eps = Obliquity
  data = [

   // 1-10 
        0   0   0   0   1  -171996.0  -174.2   92025.0     8.9 
        0   0   0   0   2     2062.0     0.2    -895.0     0.5 
       -2   0   2   0   1       46.0     0.0     -24.0     0.0 
        2   0  -2   0   0       11.0     0.0       0.0     0.0 
       -2   0   2   0   2       -3.0     0.0       1.0     0.0 
        1  -1   0  -1   0       -3.0     0.0       0.0     0.0 
        0  -2   2  -2   1       -2.0     0.0       1.0     0.0 
        2   0  -2   0   1        1.0     0.0       0.0     0.0 
        0   0   2  -2   2   -13187.0    -1.6    5736.0    -3.1 
        0   1   0   0   0     1426.0    -3.4      54.0    -0.1 

   // 11-20 
        0   1   2  -2   2     -517.0     1.2     224.0    -0.6 
        0  -1   2  -2   2      217.0    -0.5     -95.0     0.3 
        0   0   2  -2   1      129.0     0.1     -70.0     0.0 
        2   0   0  -2   0       48.0     0.0       1.0     0.0 
        0   0   2  -2   0      -22.0     0.0       0.0     0.0 
        0   2   0   0   0       17.0    -0.1       0.0     0.0 
        0   1   0   0   1      -15.0     0.0       9.0     0.0 
        0   2   2  -2   2      -16.0     0.1       7.0     0.0 
        0  -1   0   0   1      -12.0     0.0       6.0     0.0 
       -2   0   0   2   1       -6.0     0.0       3.0     0.0 

   // 21-30 
        0  -1   2  -2   1       -5.0     0.0       3.0     0.0 
        2   0   0  -2   1        4.0     0.0      -2.0     0.0 
        0   1   2  -2   1        4.0     0.0      -2.0     0.0 
        1   0   0  -1   0       -4.0     0.0       0.0     0.0 
        2   1   0  -2   0        1.0     0.0       0.0     0.0 
        0   0  -2   2   1        1.0     0.0       0.0     0.0 
        0   1  -2   2   0       -1.0     0.0       0.0     0.0 
        0   1   0   0   2        1.0     0.0       0.0     0.0 
       -1   0   0   1   1        1.0     0.0       0.0     0.0 
        0   1   2  -2   0       -1.0     0.0       0.0     0.0 

   // 31-40 
        0   0   2   0   2    -2274.0    -0.2     977.0    -0.5 
        1   0   0   0   0      712.0     0.1      -7.0     0.0 
        0   0   2   0   1     -386.0    -0.4     200.0     0.0 
        1   0   2   0   2     -301.0     0.0     129.0    -0.1 
        1   0   0  -2   0     -158.0     0.0      -1.0     0.0 
       -1   0   2   0   2      123.0     0.0     -53.0     0.0 
        0   0   0   2   0       63.0     0.0      -2.0     0.0 
        1   0   0   0   1       63.0     0.1     -33.0     0.0 
       -1   0   0   0   1      -58.0    -0.1      32.0     0.0 
       -1   0   2   2   2      -59.0     0.0      26.0     0.0 

   // 41-50 
        1   0   2   0   1      -51.0     0.0      27.0     0.0 
        0   0   2   2   2      -38.0     0.0      16.0     0.0 
        2   0   0   0   0       29.0     0.0      -1.0     0.0 
        1   0   2  -2   2       29.0     0.0     -12.0     0.0 
        2   0   2   0   2      -31.0     0.0      13.0     0.0 
        0   0   2   0   0       26.0     0.0      -1.0     0.0 
       -1   0   2   0   1       21.0     0.0     -10.0     0.0 
       -1   0   0   2   1       16.0     0.0      -8.0     0.0 
        1   0   0  -2   1      -13.0     0.0       7.0     0.0 
       -1   0   2   2   1      -10.0     0.0       5.0     0.0 

   // 51-60 
        1   1   0  -2   0       -7.0     0.0       0.0     0.0 
        0   1   2   0   2        7.0     0.0      -3.0     0.0 
        0  -1   2   0   2       -7.0     0.0       3.0     0.0 
        1   0   2   2   2       -8.0     0.0       3.0     0.0 
        1   0   0   2   0        6.0     0.0       0.0     0.0 
        2   0   2  -2   2        6.0     0.0      -3.0     0.0 
        0   0   0   2   1       -6.0     0.0       3.0     0.0 
        0   0   2   2   1       -7.0     0.0       3.0     0.0 
        1   0   2  -2   1        6.0     0.0      -3.0     0.0 
        0   0   0  -2   1       -5.0     0.0       3.0     0.0 

   // 61-70 
        1  -1   0   0   0        5.0     0.0       0.0     0.0 
        2   0   2   0   1       -5.0     0.0       3.0     0.0 
        0   1   0  -2   0       -4.0     0.0       0.0     0.0 
        1   0  -2   0   0        4.0     0.0       0.0     0.0 
        0   0   0   1   0       -4.0     0.0       0.0     0.0 
        1   1   0   0   0       -3.0     0.0       0.0     0.0 
        1   0   2   0   0        3.0     0.0       0.0     0.0 
        1  -1   2   0   2       -3.0     0.0       1.0     0.0 
       -1  -1   2   2   2       -3.0     0.0       1.0     0.0 
       -2   0   0   0   1       -2.0     0.0       1.0     0.0 

   // 71-80 
        3   0   2   0   2       -3.0     0.0       1.0     0.0 
        0  -1   2   2   2       -3.0     0.0       1.0     0.0 
        1   1   2   0   2        2.0     0.0      -1.0     0.0 
       -1   0   2  -2   1       -2.0     0.0       1.0     0.0 
        2   0   0   0   1        2.0     0.0      -1.0     0.0 
        1   0   0   0   2       -2.0     0.0       1.0     0.0 
        3   0   0   0   0        2.0     0.0       0.0     0.0 
        0   0   2   1   2        2.0     0.0      -1.0     0.0 
       -1   0   0   0   2        1.0     0.0      -1.0     0.0 
        1   0   0  -4   0       -1.0     0.0       0.0     0.0 

   // 81-90 
       -2   0   2   2   2        1.0     0.0      -1.0     0.0 
       -1   0   2   4   2       -2.0     0.0       1.0     0.0 
        2   0   0  -4   0       -1.0     0.0       0.0     0.0 
        1   1   2  -2   2        1.0     0.0      -1.0     0.0 
        1   0   2   2   1       -1.0     0.0       1.0     0.0 
       -2   0   2   4   2       -1.0     0.0       1.0     0.0 
       -1   0   4   0   2        1.0     0.0       0.0     0.0 
        1  -1   0  -2   0        1.0     0.0       0.0     0.0 
        2   0   2  -2   1        1.0     0.0      -1.0     0.0 
        2   0   2   2   2       -1.0     0.0       0.0     0.0 

   // 91-100 
        1   0   0   2   1       -1.0     0.0       0.0     0.0 
        0   0   4  -2   2        1.0     0.0       0.0     0.0 
        3   0   2  -2   2        1.0     0.0       0.0     0.0 
        1   0   2  -2   0       -1.0     0.0       0.0     0.0 
        0   1   2   0   1        1.0     0.0       0.0     0.0 
       -1  -1   0   2   1        1.0     0.0       0.0     0.0 
        0   0  -2   0   1       -1.0     0.0       0.0     0.0 
        0   0   2  -1   2       -1.0     0.0       0.0     0.0 
        0   1   0   2   0       -1.0     0.0       0.0     0.0 
        1   0  -2  -2   0       -1.0     0.0       0.0     0.0 

   // 101-106 
        0  -1   2   0   1       -1.0     0.0       0.0     0.0 
        1   1   0  -2   1       -1.0     0.0       0.0     0.0 
        1   0  -2   2   0       -1.0     0.0       0.0     0.0 
        2   0   0   2   0        1.0     0.0       0.0     0.0 
        0   0   2   4   2       -1.0     0.0       0.0     0.0 
        0   1   0   1   0        1.0     0.0       0.0     0.0 
   ];

  sp = data(:,6);      // coeff for sin(psi)
  spt = data(:,7);     // coeff for sin(psi)*t
  ce = data(:,8);      // coeff for cos(eps)
  cet = data(:,9);     // coeff for cos(eps)*t 
  LS = data(:,1:5);    // coefficients of l l' F D Omega 
endfunction