// 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