//  Copyright (c) CNES  2010
//
//  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 [angles, angles_der] = CL_mod_IAUBodyAngles(cjd, body ,tt_tref)
// Rotation angles defining the IAU body-fixed rotating frame
//
// Calling Sequence
// [angles, angles_der] = CL_mod_IAUBodyAngles(cjd, body [,tt_tref])
//
// Description
// <itemizedlist><listitem>
// <p>This function computes the three rotation angles that define the IAU body-fixed rotating frame.  </p>
// <p>The body fixed frame axes are obtained from the ICRS axes by applying 3 successive rotations of respective angles 
// [alpha0+pi/2, pi/2-delta0, W] around the respectives axes [3,1,3], where:  </p>
// <p>- alpha0 and delta0: right ascension and declination of the North pole of the body in the ICRS frame, </p>
// <p>- W: angle defining the prime meridian of the body. </p>
// <p></p> 
// <p>The function also optionally computes the time derivatives of these angles: angles_der = [alpha0dot;delta0dot;Wdot]. </p>
// <p></p>
// <p>Available bodies are: "Mercury", "Venus", "Earth", "Mars", "Jupiter", "Saturn", "Uranus", "Neptune", 
// "Sun" and "Moon". </p>
// <p></p></listitem>
// <listitem>
// <p>Notes:</p>
// <p>- The date should be given in TDB time scale. In practice using TT time scale leads to a negligible 
// error (about 2 ms).</p>
// <p>- Body names are case sensitive</p>
// <p>- This implementation of the function is compliant with the 2009 IAU report. (See bibliography)</p>
// </listitem>
// </itemizedlist>
//
// Parameters
// cjd: Modified (1950.0) julian day (Time scale: TREF) (1xN)
// body: (string) Name of the body. ("Mercury","Venus","Earth","Mars","Jupiter","Saturn","Uranus", "Neptune", "Sun" or "Moon") (1x1)
// tt_tref: (optional) TT-TREF [seconds]. Default is %CL_TT_TREF. (1xN or 1x1)
// angles: Angles [alpha0;delta0;W] defining the orientation of the body [rad]. (1xN)
// angles_der: Time derivatives of Angles = [alpha0dot;delta0dot;Wdot] [rad/s]. (1xN)
//
// Authors
// CNES - DCT/SB
//
// Bibliography
// Report of the IAU/IAG working group on cartographic coordinates and rotational elements: 2009
//
// See also
// CL_fr_bodyConvertMat
// CL_fr_bodyConvert
//
// Examples
// cjd = CL_dat_cal2cjd(2010,02,03,05,35,25); // TREF time scale
//
// [angles,angles_der] = CL_mod_IAUBodyAngles(cjd, "Mars")


// Declarations:
DEG2RAD = %pi / 180;

// Code:
if (argn(2) <> 2 & argn(2) <> 3) 
  CL__error("Invalid number of input arguments"); 
end

if (~exists("tt_tref", "local")); tt_tref = CL__dataGetEnv("TT_TREF"); end
if (typeof(body) <> "string"); CL__error("Invalid input argument ''body''"); end
if (size(body,"*") <> 1); CL__error("Invalid input argument ''body''"); end

if (cjd == [])
  angles = [];
  angles_der = [];
  return; // <-- RETURN
end


// t = Number of days since reference epoch (J2000.0 = 1 January 2000 12:00:00 TDB)
// NB : To avoid loosing numerical precision: avoid conversion to julian days
t = (cjd - 18262.5) + tt_tref/86400 ;

// T = Number of centuries since reference epoch
T = t / 36525;

// RA,DEC and W are expressed in degrees.
// RAdot and DECdot are expressed in degrees/century
// Wdot is expressed in degrees/days


if (body == "Mercury")
  // Mercury: The 20 degrees meridian is defined by the crater Hun Kal
  
  // M1 to M5 in radians
  M1 = (174.791086 + 4.092335*t) * DEG2RAD;
  M2 = (349.582171 + 8.184670*t) * DEG2RAD;
  M3 = (164.373257 + 12.277005*t) * DEG2RAD;
  M4 = (339.164343 + 16.369340*t) * DEG2RAD;
  M5 = (153.955429 + 20.461675*t) * DEG2RAD;
  
  // M1d to M5d in radians/days
  M1d = 4.092335 * DEG2RAD;
  M2d = 8.184670 * DEG2RAD;
  M3d = 12.277005 * DEG2RAD;
  M4d = 16.369340 * DEG2RAD;
  M5d = 20.461675 * DEG2RAD;

  RA = 281.0097 - 0.0328*T; 
  DEC = 61.4143 - 0.0049*T;
  W = 329.5469 + 6.1385025*t + 0.00993822 * sin(M1) - ..
      0.00104581 * sin(M2) - 0.00010280 * sin(M3) - ..
      0.00002364 * sin(M4) - 0.00000532 * sin(M5);

  RAdot = -0.0328*ones(t); 
  DECdot = -0.0049*ones(t); 
  Wdot = 6.1385025 + 0.00993822 * M1d * cos(M1) - ..
      0.00104581 * M2d * cos(M2) - 0.00010280 * M3d * cos(M3) - ..
      0.00002364 * M4d * cos(M4) - 0.00000532 * M5d * cos(M5);
  
elseif (body == "Venus")
  // Venus: The 0 degrees meridian is defined by the central peak in the crater Ariadne
  RA = 272.76*ones(t);
  DEC = 67.16*ones(t);
  W = 160.20 - 1.4813688*t;
  
  RAdot = zeros(T);
  DECdot = zeros(T);
  Wdot = -1.4813688*ones(t);
  
elseif (body == "Earth")
  RA = -0.641*T;
  DEC = 90.00 - 0.557*T;
  W = 190.147 + 360.9856235*t;
  
  RAdot = -0.641*ones(t);
  DECdot = -0.557*ones(t);
  Wdot = 360.9856235*ones(t);
  
elseif (body == "Mars")
  // Mars: The 0 degrees meridian is defined by the crater Airy-0
  RA = 317.68143 - 0.1061*T;
  DEC = 52.88650 - 0.0609*T;
  W = 176.630 + 350.89198226*t;
  
  RAdot = -0.1061*ones(t);
  DECdot = -0.0609*ones(t);
  Wdot = 350.89198226*ones(t);
  
elseif (body == "Jupiter")
  // The equations for W for Jupiter, Saturn, Uranus and Neptune refer to the rotation of their magnetic fields

  // Ja to Je in radians
  Ja = (99.360714 + 4850.4046*T) * DEG2RAD;
  Jb = (175.895369 + 1191.9605*T) * DEG2RAD;
  Jc = (300.323162 + 262.5475*T) * DEG2RAD;
  Jd = (114.012305 + 6070.2476*T) * DEG2RAD;
  Je = (49.511251 + 64.3000*T) * DEG2RAD;
  
  // Jad to Jed in radians/century
  Jad = 4850.4046 * DEG2RAD;
  Jbd = 1191.9605 * DEG2RAD;
  Jcd = 262.5475 * DEG2RAD;
  Jdd = 6070.2476 * DEG2RAD;
  Jed = 64.3000 * DEG2RAD;

  RA = 268.056595 - 0.006499*T + 0.000117*sin(Ja) + 0.000938*sin(Jb) + 0.001432*sin(Jc) + 0.000030*sin(Jd) + 0.002150*sin(Je);
  DEC = 64.495303 + 0.002413*T + 0.000050*cos(Ja) + 0.000404*cos(Jb) + 0.000617*cos(Jc) - 0.000013*cos(Jd) + 0.000926*cos(Je);
  W = 284.95 + 870.536*t;
  
  RAdot = -0.006499 + 0.000117*Jad*cos(Ja) + 0.000938*Jbd*cos(Jb) + 0.001432*Jcd*cos(Jc) + 0.000030*Jdd*cos(Jd) + 0.002150*Jed*cos(Je);
  DECdot = 0.002413 - 0.000050*Jad*sin(Ja) - 0.000404*Jbd*sin(Jb) - 0.000617*Jcd*sin(Jc) + 0.000013*Jdd*sin(Jd) - 0.000926*Jed*sin(Je);
  Wdot = 870.536*ones(t);
  
elseif (body == "Saturn")
  // The equations for W for Jupiter, Saturn, Uranus and Neptune refer to the rotation of their magnetic fields
  RA = 40.589 - 0.036*T;
  DEC = 83.537 - 0.004*T;
  W = 38.90 + 810.7939024*t;
  
  RAdot = -0.036*ones(t);
  DECdot = -0.004*ones(t);
  Wdot = 810.7939024*ones(t);
  
elseif (body == "Uranus")
  // The equations for W for Jupiter, Saturn, Uranus and Neptune refer to the rotation of their magnetic fields
  RA = 257.311*ones(t) ;
  DEC = -15.175*ones(t);
  W = 203.81 - 501.1600928*t;
  
  RAdot = zeros(T);
  DECdot = zeros(T);
  Wdot = -501.1600928*ones(t);
  
elseif (body == "Neptune")
  // The equations for W for Jupiter, Saturn, Uranus and Neptune refer to the rotation of their magnetic fields
  N = (357.85 + 52.316*T) * DEG2RAD; // rad
  Ndot = 52.316*ones(t) * DEG2RAD; // rad/century
  
  RA = 299.36 + 0.70 * sin(N);
  DEC = 43.46 - 0.51 * cos(N);
  W = 253.18 + 536.3128492 * t - 0.48 * sin(N);
  
  RAdot = 0.70 * Ndot .* cos(N);    // deg/century
  DECdot = -0.51 * Ndot .* -sin(N); // deg/century 
  Wdot = 536.3128492 - 0.48 * (Ndot/36525) .* cos(N); // deg/day
  
elseif (body == "Sun")
  RA = 286.13*ones(t);
  DEC = 63.87*ones(t);
  W = 84.176 + 14.1844000*t;
  RAdot = zeros(T);
  DECdot = zeros(T);
  Wdot = 14.1844000*ones(t);
  
elseif(body == "Moon")
  // NB: These formulae are precise to approximately 150m
  
  // E1 to E13 : in radians
  E1 = (125.045 - 0.0529921*t) * DEG2RAD;
  E2 = (250.089 - 0.1059842*t) * DEG2RAD;
  E3 = (260.008 + 13.0120009*t) * DEG2RAD;
  E4 = (176.625 + 13.3407154*t) * DEG2RAD;
  E5 = (357.529 + 0.9856003*t) * DEG2RAD;
  E6 = (311.589 + 26.4057084*t) * DEG2RAD;
  E7 = (134.963 + 13.0649930*t) * DEG2RAD;
  E8 = (276.617 + 0.3287146*t) * DEG2RAD;
  E9 = (34.226 + 1.7484877*t) * DEG2RAD;
  E10 = (15.134 - 0.1589763*t) * DEG2RAD;
  E11 = (119.743 + 0.0036096*t) * DEG2RAD;
  E12 = (239.961 + 0.1643573*t) * DEG2RAD;
  E13 = (25.053 + 12.9590088*t) * DEG2RAD;
  
  // E1dot to E13dot : in radians/days
  E1dot = -0.0529921 * DEG2RAD;
  E2dot = -0.1059842 * DEG2RAD;
  E3dot = 13.0120009 * DEG2RAD;
  E4dot = 13.3407154 * DEG2RAD;
  E5dot = 0.9856003 * DEG2RAD;
  E6dot = 26.4057084 * DEG2RAD;
  E7dot = 13.0649930 * DEG2RAD;
  E8dot = 0.3287146 * DEG2RAD;
  E9dot = 1.7484877 * DEG2RAD;
  E10dot = -0.1589763 * DEG2RAD;
  E11dot = 0.0036096 * DEG2RAD;
  E12dot = 0.1643573 * DEG2RAD;
  E13dot = 12.9590088 * DEG2RAD;

  RA = 269.9949 + 0.0031*T - 3.8787*sin(E1) - 0.1204*sin(E2) + 0.0700*sin(E3) + ..
              -0.0172*sin(E4) + 0.0072*sin(E6) - 0.0052*sin(E10) + 0.0043*sin(E13);
              
  DEC = 66.5392 + 0.0130*T + 1.5419*cos(E1) + 0.0239*cos(E2) + ..
              -0.0278*cos(E3) + 0.0068*cos(E4) - 0.0029*cos(E6) + ..
               0.0009*cos(E7) + 0.0008*cos(E10) - 0.0009*cos(E13);
               
  W = 38.3213 + 13.17635815*t - 1.4e-12*t.^2 + 3.5610*sin(E1) + 0.1208*sin(E2) + ..
            -0.0642*sin(E3) + 0.0158*sin(E4) + 0.0252*sin(E5) - 0.0066*sin(E6) + ..
            -0.0047*sin(E7) - 0.0046*sin(E8) + 0.0028*sin(E9) + 0.0052*sin(E10) + ..
             0.0040*sin(E11) + 0.0019*sin(E12) - 0.0044*sin(E13);

  // deg/century
  RAdot = 0.0031 + 36525 * ( ..
       -3.8787*E1dot*cos(E1) - 0.1204*E2dot*cos(E2) + 0.0700*E3dot*cos(E3) + ..
       -0.0172*E4dot*cos(E4) + 0.0072*E6dot*cos(E6) -0.0052*E10dot*cos(E10)+0.0043*E13dot*cos(E13));
  
  // deg/century
  DECdot = 0.0130 + 36525 * ( ..
      -1.5419*E1dot*sin(E1) - 0.0239*E2dot*sin(E2) + ..
       0.0278*E3dot*sin(E3) - 0.0068*E4dot*sin(E4) + 0.0029*E6dot*sin(E6) + ..
       -0.0009*E7dot*sin(E7) - 0.0008*E10dot*sin(E10) + 0.0009*E13dot*sin(E13));
  
  // deg/day
  Wdot = 13.17635815 - 1.4e-12*2*t + 3.5610*E1dot*cos(E1) + 0.1208*E2dot*cos(E2) + ..
            -0.0642*E3dot*cos(E3) + 0.0158*E4dot*cos(E4) + 0.0252*E5dot*cos(E5) - 0.0066*E6dot*cos(E6) + ..
            -0.0047*E7dot*cos(E7) - 0.0046*E8dot*cos(E8) + 0.0028*E9dot*cos(E9) + 0.0052*E10dot*cos(E10) + ..
             0.0040*E11dot*cos(E11) + 0.0019*E12dot*cos(E12) - 0.0044*E13dot*cos(E13);          
else
  CL__error("Invalid body name (case sensitive!)");
end


// Output angles are in radians
W = pmodulo(W,360);
angles = [RA; DEC; W] * DEG2RAD;

// Output time derivatives of angles are in rad/s
angles_der = [RAdot/36525; DECdot/36525; Wdot] * DEG2RAD / 86400;

endfunction