//  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 [NUT,F,NUTP,FP,NUTPP,FPP] = CL_mod_nutationAngles(jj_tai,sec_tai, flag)
// Nutation angles (Wahr IAU 1980) - DEPRECATED
//
// Calling Sequence
// [NUT,F,NUTP,FP,NUTPP,FPP] = CL_mod_nutationAngles(jj_tai,sec_tai [,flag])
//
// Description
// <itemizedlist><listitem>
// <p>This function is deprecated. </p>
// <p></p></listitem>
// <listitem>
// <p>Computes the nutation angles (delta longitude, delta obliquity)
// according to Wahr nutation theory. </p>
// </listitem>
// </itemizedlist>
//
// Parameters
// jj_tai : (integer) Modified (1950.0) Julian day (1xN)
// sec_tai : seconds  (TAI time scale) (1xN)
// flag   : (optional) "n" -> NUT,F; "p" NUT,F and NUTP,Fp; "s" -> F, Fp, Fpp and NUT,NUTp,NUTPP (default depends on number of output variables)
// NUT    : longitude, obliquity NUT=[DPSI;DEPS] [rad] (2xN)
// F      : fundamentals arguments of nutation [rad] (5xN)
// NUTP : (optional with FP) first time derivative of NUT [rad/s] (2xN)
// NUTPP : (optional with FPP) second time derivative of NUT [rad/s^2] (5xN)
// FP : (optional with NUTP) first time derivative of F [rad/s] (2xN)
// FPP : (optional with NUTPP) second time derivative of F [rad/s^2] (5xN)
//
// Authors
// CNES - DCT/SB
//
// Bibliography
// 1) IERS Conventions (1996), Dennis D. McCarthy
//
// See also
// CL_mod_nutationArg
//
// Examples
// jj_tai = [19500:1:20500];
// ss_tai = 43200*ones(jj_tai);
// [NUT,F,NUTP,FP,NUTPP,FPP]=CL_mod_nutationAngles(jj_tai,ss_tai,"s");
// plot(jj_tai'+ss_tai'/86400,NUT');
//

// Declarations:


// Code:

CL__warnDeprecated(); // deprecated function

[lhs,rhs]=argn(0);

if (~exists('flag','local'))
  select lhs
    case 2
      flag='n';
    case 4
      flag='p';
    case 6
      flag='s';
    else
      CL__error("check number of output arguments");
  end
end

if (rhs>3)|(rhs==1) CL__error("check number of input arguments"); end,

if (flag~="s" & flag~="p" & flag~="n") then
  CL__error("value of parameter ''flag'' = "+flag+" unknown");
end

J2000 = 36525. / 2.;
TT_MOINS_TAI = 32.184;
t = (jj_tai - J2000 + (TT_MOINS_TAI+sec_tai)/86400) / 36525;
t = t(:)';

//coefficients de la nutation
A=[0.,0.,0.,0.,1.,-171996.,-174.2,92025.,8.9
0.,0.,0.,0.,2.,2062.,0.2,-895.,0.5
-2.,0.,2.,0.,1.,46.,0.0,-24.,0.0
2.,0.,-2.,0.,0.,11.,0.0,0.,0.0
-2.,0.,2.,0.,2.,-3.,0.0,1.,0.0
1.,-1.,0.,-1.,0.,-3.,0.0,0.,0.0
0.,-2.,2.,-2.,1.,-2.,0.0,1.,0.0
2.,0.,-2.,0.,1.,1.,0.0,0.,0.0
0.,0.,2.,-2.,2.,-13187.,-1.6,5736.,-3.1
0.,1.,0.,0.,0.,1426.,-3.4,54.,-0.1
0.,1.,2.,-2.,2.,-517.,1.2,224.,-0.6
0.,-1.,2.,-2.,2.,217.,-0.5,-95.,0.3
0.,0.,2.,-2.,1.,129.,0.1,-70.,0.0
2.,0.,0.,-2.,0.,48.,0.0,1.,0.0
0.,0.,2.,-2.,0.,-22.,0.0,0.,0.0
0.,2.,0.,0.,0.,17.,-0.1,0.,0.0
0.,1.,0.,0.,1.,-15.,0.0,9.,0.0
0.,2.,2.,-2.,2.,-16.,0.1,7.,0.0
0.,-1.,0.,0.,1.,-12.,0.0,6.,0.0
-2.,0.,0.,2.,1.,-6.,0.0,3.,0.0
0.,-1.,2.,-2.,1.,-5.,0.0,3.,0.0
2.,0.,0.,-2.,1.,4.,0.0,-2.,0.0
0.,1.,2.,-2.,1.,4.,0.0,-2.,0.0
1.,0.,0.,-1.,0.,-4.,0.0,0.,0.0
2.,1.,0.,-2.,0.,1.,0.0,0.,0.0
0.,0.,-2.,2.,1.,1.,0.0,0.,0.0
0.,1.,-2.,2.,0.,-1.,0.0,0.,0.0
0.,1.,0.,0.,2.,1.,0.0,0.,0.0
-1.,0.,0.,1.,1.,1.,0.0,0.,0.0
0.,1.,2.,-2.,0.,-1.,0.0,0.,0.0
0.,0.,2.,0.,2.,-2274.,-0.2,977.,-0.5
1.,0.,0.,0.,0.,712.,0.1,-7.,0.0
0.,0.,2.,0.,1.,-386.,-0.4,200.,0.0
1.,0.,2.,0.,2.,-301.,0.0,129.,-0.1
1.,0.,0.,-2.,0.,-158.,0.0,-1.,0.0
-1.,0.,2.,0.,2.,123.,0.0,-53.,0.0
0.,0.,0.,2.,0.,63.,0.0,-2.,0.0
1.,0.,0.,0.,1.,63.,0.1,-33.,0.0
-1.,0.,0.,0.,1.,-58.,-0.1,32.,0.0
-1.,0.,2.,2.,2.,-59.,0.0,26.,0.0
1.,0.,2.,0.,1.,-51.,0.0,27.,0.0
0.,0.,2.,2.,2.,-38.,0.0,16.,0.0
2.,0.,0.,0.,0.,29.,0.0,-1.,0.0
1.,0.,2.,-2.,2.,29.,0.0,-12.,0.0
2.,0.,2.,0.,2.,-31.,0.0,13.,0.0
0.,0.,2.,0.,0.,26.,0.0,-1.,0.0
-1.,0.,2.,0.,1.,21.,0.0,-10.,0.0
-1.,0.,0.,2.,1.,16.,0.0,-8.,0.0
1.,0.,0.,-2.,1.,-13.,0.0,7.,0.0
-1.,0.,2.,2.,1.,-10.,0.0,5.,0.0
1.,1.,0.,-2.,0.,-7.,0.0,0.,0.0
0.,1.,2.,0.,2.,7.,0.0,-3.,0.0
0.,-1.,2.,0.,2.,-7.,0.0,3.,0.0
1.,0.,2.,2.,2.,-8.,0.0,3.,0.0
1.,0.,0.,2.,0.,6.,0.0,0.,0.0
2.,0.,2.,-2.,2.,6.,0.0,-3.,0.0
0.,0.,0.,2.,1.,-6.,0.0,3.,0.0
0.,0.,2.,2.,1.,-7.,0.0,3.,0.0
1.,0.,2.,-2.,1.,6.,0.0,-3.,0.0
0.,0.,0.,-2.,1.,-5.,0.0,3.,0.0
1.,-1.,0.,0.,0.,5.,0.0,0.,0.0
2.,0.,2.,0.,1.,-5.,0.0,3.,0.0
0.,1.,0.,-2.,0.,-4.,0.0,0.,0.0
1.,0.,-2.,0.,0.,4.,0.0,0.,0.0
0.,0.,0.,1.,0.,-4.,0.0,0.,0.0
1.,1.,0.,0.,0.,-3.,0.0,0.,0.0
1.,0.,2.,0.,0.,3.,0.0,0.,0.0
1.,-1.,2.,0.,2.,-3.,0.0,1.,0.0
-1.,-1.,2.,2.,2.,-3.,0.0,1.,0.0
-2.,0.,0.,0.,1.,-2.,0.0,1.,0.0
3.,0.,2.,0.,2.,-3.,0.0,1.,0.0
0.,-1.,2.,2.,2.,-3.,0.0,1.,0.0
1.,1.,2.,0.,2.,2.,0.0,-1.,0.0
-1.,0.,2.,-2.,1.,-2.,0.0,1.,0.0
2.,0.,0.,0.,1.,2.,0.0,-1.,0.0
1.,0.,0.,0.,2.,-2.,0.0,1.,0.0
3.,0.,0.,0.,0.,2.,0.0,0.,0.0
0.,0.,2.,1.,2.,2.,0.0,-1.,0.0
-1.,0.,0.,0.,2.,1.,0.0,-1.,0.0
1.,0.,0.,-4.,0.,-1.,0.0,0.,0.0
-2.,0.,2.,2.,2.,1.,0.0,-1.,0.0
-1.,0.,2.,4.,2.,-2.,0.0,1.,0.0
2.,0.,0.,-4.,0.,-1.,0.0,0.,0.0
1.,1.,2.,-2.,2.,1.,0.0,-1.,0.0
1.,0.,2.,2.,1.,-1.,0.0,1.,0.0
-2.,0.,2.,4.,2.,-1.,0.0,1.,0.0
-1.,0.,4.,0.,2.,1.,0.0,0.,0.0
1.,-1.,0.,-2.,0.,1.,0.0,0.,0.0
2.,0.,2.,-2.,1.,1.,0.0,-1.,0.0
2.,0.,2.,2.,2.,-1.,0.0,0.,0.0
1.,0.,0.,2.,1.,-1.,0.0,0.,0.0
0.,0.,4.,-2.,2.,1.,0.0,0.,0.0
3.,0.,2.,-2.,2.,1.,0.0,0.,0.0
1.,0.,2.,-2.,0.,-1.,0.0,0.,0.0
0.,1.,2.,0.,1.,1.,0.0,0.,0.0
-1.,-1.,0.,2.,1.,1.,0.0,0.,0.0
0.,0.,-2.,0.,1.,-1.,0.0,0.,0.0
0.,0.,2.,-1.,2.,-1.,0.0,0.,0.0
0.,1.,0.,2.,0.,-1.,0.0,0.,0.0
1.,0.,-2.,-2.,0.,-1.,0.0,0.,0.0
0.,-1.,2.,0.,1.,-1.,0.0,0.,0.0
1.,1.,0.,-2.,1.,-1.,0.0,0.,0.0
1.,0.,-2.,2.,0.,-1.,0.0,0.,0.0
2.,0.,0.,2.,0.,1.,0.0,0.,0.0
0.,0.,2.,4.,2.,-1.,0.0,0.,0.0
0.,1.,0.,1.,0.,1.,0.0,0.,0.0];

//argnut: calcul des arguments astronomique pour chaque armonique
//attention... FP en rad par second , FPP en rad per second^2
siecle2sec=36525*86400;
arcsec2rad=(1/3600)*(%pi/180);
DPSIP=[];DEPSP=[];FP=[];DPSIPP=[];DEPSPP=[];FPP=[];

[F,FP,FPP]=CL_mod_nutationArg(t,flag);
DPSI   = sum((A(:,6)*ones(t) + A(:,7)*t).*sin(A(:,1:5)*F),'r')*arcsec2rad*1d-4;
DEPS   = sum((A(:,8)*ones(t) + A(:,9)*t).*cos(A(:,1:5)*F),'r')*arcsec2rad*1d-4;
NUT=[DPSI; DEPS];

NUTP=[];
if (flag=="p" | flag=="s") then
  FP=FP*siecle2sec;
  DPSIP  = sum((A(:,7)*ones(t)).*sin(A(:,1:5)*F) + (A(:,6)*ones(t) + ...
    A(:,7)*t).*( A(:,1:5)*FP).*cos(A(:,1:5)*F),'r')*arcsec2rad*1d-4;
  DPSIP  = DPSIP/(siecle2sec);
  DEPSP  = sum((A(:,9)*ones(t)).*cos(A(:,1:5)*F) + (A(:,8)*ones(t) + ...
    A(:,9)*t).*(-A(:,1:5)*FP).*sin(A(:,1:5)*F),'r')*arcsec2rad*1d-4;
  DEPSP  = DEPSP/(siecle2sec);
  NUTP=[DPSIP; DEPSP];

end

NUTPP=[];
if flag=="s" then
  FPP=FPP*siecle2sec*siecle2sec;
  DPSIPP = sum(2*(A(:,7)*ones(t)).*(A(:,1:5)*FP).*cos(A(:,1:5)*F) + ...
       (A(:,6)*ones(t) + A(:,7)*t).*(A(:,1:5)*FPP).*cos(A(:,1:5)*F) - ...
       (A(:,6)*ones(t) + A(:,7)*t).*(A(:,1:5)*FP).*(A(:,1:5)*FP).*sin(A(:,1:5)*F) ,'r')*arcsec2rad*1d-4;
  DPSIPP=DPSIPP/(siecle2sec^2);
  DEPSPP = sum(-2*(A(:,9)*ones(t)).*(A(:,1:5)*FP).*sin(A(:,1:5)*F) + ...
       (A(:,8)*ones(t) + A(:,9)*t).*(-A(:,1:5)*FPP).*sin(A(:,1:5)*F) + ...
       (A(:,8)*ones(t) + A(:,9)*t).*(-A(:,1:5)*FP).*(A(:,1:5)*FP).*cos(A(:,1:5)*F) ,'r')*arcsec2rad*1d-4;
  DEPSPP=DEPSPP/(siecle2sec^2);
  NUTPP=[DPSIPP; DEPSPP];
end
FP=FP/siecle2sec;
FPP=FPP/(siecle2sec*siecle2sec);

endfunction