//  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 [Fls, Flsdot] = CL__iers_fundArgsLS(jd, cder)
// Compute the fundamental arguments for lunisolar nutation (IERS2010).
//
// Calling Sequence
// [Fls, Flsdot] = CL__iers_fundArgsLS(jd [,cder])
//
// Description
// <itemizedlist><listitem>
// Compute the fundamental arguments of lunisolar nutation <b>Fls</b> and their time derivatives
// - Mean Anomaly of the Moon (l)
// - Mean Anomaly of the Sun (l')
// - Mean longitude of the Moon minus that of the ascending node of the Moon (F=L-omega)
// - Mean Elongation of the Moon from the Sun (D)
// - Mean longitude of the ascending node of the Moon (omega)
// 
// Note : used in nutation and series (X,Y,s) computations
//
// </itemizedlist>
//
// Parameters
// jd: Two-part julian day (Time scale: TT) (2xN)
// cder: (boolean, optional) Option to compute derivative. If cder is %f, the derivative will be set to []. Default is %t. (1x1)
// Fls: Arguments of the lunisolar nutation [rad] (5xN)
// Flsdot: Derivatives of arguments of the lunisolar nutation [rad/s] (5xN)
//
// Authors
// CNES - DCT/SB
//
// Bibliography
// 1) Technical Note 36, IERS, 2010

// Declarations:


// Code:

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

ARCSEC_TO_RAD = %pi / (180*3600);
DEG_TO_RAD = %pi/180;

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

N = length(t);

// ---------------------------------------------
// Fundamental argument for the lunisolar nutation (eq 5.43, p67, paragraph 5.7.2)
// ---------------------------------------------
Fls = zeros(5,N);
   
// Mean Anomaly of the Moon
Fls(1,:) = 134.96340251 * DEG_TO_RAD  + t .* (1717915923.2178 + t .* (31.8792 + t .* (0.051635 + t .* (-0.00024470 )))) * ARCSEC_TO_RAD;

// Mean Anomaly of the Sun
Fls(2,:) = 357.52910918 * DEG_TO_RAD  + t .* (129596581.0481 + t .* (-0.5532 + t .* (0.000136 + t .* (-0.00001149 )))) * ARCSEC_TO_RAD;

// Mean longitude of the Moon minus that of the ascending node (L-omega)
Fls(3,:) = 93.27209062 * DEG_TO_RAD + t .* (1739527262.8478 + t .* (-12.7512 + t .* (-0.001037 + t .* (0.00000417 )))) * ARCSEC_TO_RAD;

// Mean Elongation of the Moon from the Sun
Fls(4,:) = 297.85019547 * DEG_TO_RAD  + t .* (1602961601.2090 + t .* (-6.3706 + t .* (0.006593 + t .* (-0.00003169 )))) * ARCSEC_TO_RAD;

// Mean longitude of the ascending node of the Moon
Fls(5,:) = 125.04455501 * DEG_TO_RAD  + t .* (-6962890.5431 + t .* (7.4722 + t .* (0.007702 + t .* (-0.00005939 )))) * ARCSEC_TO_RAD;

Fls = CL_rMod(Fls, 2*%pi);


// ------------------
// Time derivatives
// ------------------
Flsdot = [];
if (cder)
  Flsdot = zeros(5,N);
   
  // Mean Anomaly of the Moon
  Flsdot(1,:) = 1717915923.2178  + t .* (2*31.8792 + t .* (3*0.051635 + t .* (4*-0.00024470 )));

  // Mean Anomaly of the Sun
  Flsdot(2,:) = 129596581.0481  + t .* (2*-0.5532 + t .* (3*0.000136 + t .* (4*-0.00001149 )));

  // Mean longitude of the Moon minus that of the ascending node
  Flsdot(3,:) = 1739527262.8478 + t .* (2*-12.7512 + t .* (3*-0.001037 + t .* (4*0.00000417 )));

  // Mean Elongation of the Moon from the Sun
  Flsdot(4,:) = 1602961601.2090  + t .* (2*-6.3706 + t .* (3*0.006593 + t .* (4*-0.00003169 )));

  // Mean longitude of the ascending node of the Moon
  Flsdot(5,:) = -6962890.5431  + t .* (2*7.4722 + t .* (3*0.007702 + t .* (4*-0.00005939 )));

  // Convert from arcseconds/century to rad/s
  Flsdot = Flsdot *ARCSEC_TO_RAD / (36525 * 86400);
end

endfunction