//  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 [pos_moon_EOD,vel_moon_EOD] = CL__eph_moon_Meeus(jd,nbcoef, cvel)
// Moon position and velocity in EOD frame (Meeus algorithm)
//
// Calling Sequence
// [pos_moon_EOD,vel_moon_EOD] = CL__eph_moon_Meeus(jd,nbcoef [,cvel])
//
// Description
// <itemizedlist><listitem>
// <p>Computes Moon position and velocity relative to the ecliptic of date frame (EOD) using Meeus algorithm.</p>
// <p>The frame is centered on Earth. (pos_moon_EOD is the Earth-Moon vector).</p>
// <p>The precision can optionaly be modified using nbcoef:</p>
// <p> - nbcoef(1) is the number of sinusoidal terms used to compute longitude. (max=62)</p>
// <p> - nbcoef(2) is the number of sinusoidal terms used to compute latitude. (max=66)</p>
// <p> - nbcoef(3) is the number of sinusoidal terms used to compute distance. (max=46)</p>
// <p></p>
// <p>Note: if cvel is set to %t but only one output argument is required, the velocity is not computed.</p>
// </listitem>
// </itemizedlist>
//
// Parameters
// jd: Julian date (TT) (1xN)
// nbcoef: (integer) Number of sinusoidal terms to be used. Maximum is: [62,66,46] (1x3)
// cvel: (optional, boolean) %t if velocity should be computed. Default is %t. (1x1)
// pos_moon_EOD : Position relative to EOD frame (m) (3xN)
// vel_moon_EOD : Velocity relative to EOD frame (m/s) (3xN)
//
// Bibliography
// Astronomical Algorithms J.Meeus - Second edition 1998, Chapter 47
//
// Authors
// CNES - DCT/SB
//
// Examples
// // Julian date, in TT time scale
// jd = 21512210.3;
//
// // Highest precision
// [pos_moon_EOD,vel_moon_EOD] = CL__eph_moon_Meeus(jd, [62,66,46])
//
// // Lower precision
// [pos_moon_EOD,vel_moon_EOD] = CL__eph_moon_Meeus(jd , [9,4,3]);

// Declarations:
DEG2RAD = %pi / 180; 

// Code:
if (~exists("cvel","local")); cvel = %t; end;
if (argn(2) <> 2 & argn(2) <> 3)
  CL__error("Invalid number of input arguments");
end

// Handle [] case
if (jd == [])
  pos_moon_EOD = [];
  vel_moon_EOD = [];
  return; // <-- RETURN
end

// Compute velocity only if needed
if (argn(1) == 1)
  cvel = %f;
end

// Control on number of coefs : 
nlon = min(nbcoef(1), 62);
nlat = min(nbcoef(2), 66);
ndist = min(nbcoef(3), 46);

// Time in centuries from epoch 2000 12h (TT time scale)
TT   = (jd - 2451545.0) / 36525.0;
TT2   = TT .* TT ;
TT3   = TT2 .* TT ;
TT4   = TT2 .* TT2 ;


// Moon mean longitude (Lp), referred to the mean equinox of the date, 
// and including the constant term of the effect of light time (eq. 47.1)
Lp = 218.3164477 + 481267.88123421 * TT - 0.0015786 * TT2 + TT3 / 538841 - TT4 / 65194000 ; // deg
Lp = pmodulo(Lp * DEG2RAD, 2*%pi);

// Mean elongation of the Moon (eq. 47.2)
D = 297.8501921 + 445267.1114034 * TT - 0.0018819 * TT2 + TT3 / 545868 - TT4 / 113065000 ; // deg
D = pmodulo(D * DEG2RAD, 2*%pi);

// Sun mean anomaly (eq. 47.3)
M   = 357.5291092 + 35999.0502909 * TT - 0.0001536 * TT2 + TT3 / 24490000; // deg
M   = pmodulo(M * DEG2RAD, 2*%pi);

// Moon mean anomaly (eq. 47.4)
Mp = 134.9633964 + 477198.8675055 * TT + 0.0087414 * TT2 + TT3 / 69699 - TT4 / 14712000 ; // deg
Mp = pmodulo(Mp * DEG2RAD, 2*%pi);

// Moon argument of latitude (eq. 47.5)
F = 93.2720950 + 483202.0175233 * TT - 0.0036539 * TT2 - TT3 / 3526000 + TT4 / 863310000 ; // deg
F = pmodulo(F * DEG2RAD, 2*%pi);

// A1 : impact of Venus
A1 = 119.75 + 131.849 * TT ; // deg
A1 = pmodulo(A1 * DEG2RAD, 2*%pi);

// A2 : impact of Jupiter
A2 = 53.09 + 479264.290 * TT ; // deg
A2 = pmodulo(A2 * DEG2RAD, 2*%pi);

// A3 
A3 = 313.45 + 481266.484 * TT; // deg
A3 = pmodulo(A3 * DEG2RAD, 2*%pi);

// Term used to correct harmonics that contain (M) or (2M) (eq. 47.5)
// (eccentricity of Earth's orbit around the Sun is currently decreasing with time)
E = 1.0 - 0.002516 * TT - 0.0000074 * TT2 ;
E2 = E.^2;



// ---------------------
// Longitude computation
// ---------------------

// Magnitude of each sine (rad)
loncoefs = DEG2RAD * 1e-6 * ..
            [6288774;1274027;658314;213618;-185116; ..
              -114332; 58793; 57066; 53322;  45758; ..
              -40923;-34720;-30383; 15327; -12528; ..
               10980; 10675; 10034; 8548; -7888; .. // 20
              -6766;-5163; 4987; 4036;  3994; .. 
               3958;3861; 3665;-2689; -2602; ..
               2390;-2348; 2236;-2120; -2069; ..
               2048; 1962;-1773;-1595;  1215; .. // 40
              -1110;-892;-810; 759; -713; ..
              -700; 691; 596; 549;  537; ..
               520;-487;-399;-381;  351; ..
              -340; 330; 327;-323;  318; ..  // 60
               299; 294];

// Arguments of each sine (linear combination of angles : D, M, Mp, F, A1, A2, Lp)
lonargs = [ 0, 0, 1, 0, 0, 0, 0 ;          2, 0,-1, 0, 0, 0, 0 ; ..
            2, 0, 0, 0, 0, 0, 0 ;          0, 0, 2, 0, 0, 0, 0 ; ..
            0, 1, 0, 0, 0, 0, 0 ;          0, 0, 0, 2, 0, 0, 0 ; ..
            2, 0,-2, 0, 0, 0, 0 ;          2,-1,-1, 0, 0, 0, 0 ; ..
            2, 0, 1, 0, 0, 0, 0 ;          2,-1, 0, 0, 0, 0, 0 ; .. // 10
            0, 1,-1, 0, 0, 0, 0 ;          1, 0, 0, 0, 0, 0, 0 ; ..
            0, 1, 1, 0, 0, 0, 0 ;          2, 0, 0,-2, 0, 0, 0 ; ..
            0, 0, 1, 2, 0, 0, 0 ;          0, 0, 1,-2, 0, 0, 0 ; ..
            4, 0,-1, 0, 0, 0, 0 ;          0, 0, 3, 0, 0, 0, 0 ; ..
            4, 0,-2, 0, 0, 0, 0 ;          2, 1,-1, 0, 0, 0, 0 ; .. // 20
            2, 1, 0, 0, 0, 0, 0 ;          1, 0,-1, 0, 0, 0, 0 ; ..
            1, 1, 0, 0, 0, 0, 0 ;          2,-1, 1, 0, 0, 0, 0 ; ..
            2, 0, 2, 0, 0, 0, 0 ;          0, 0, 0, 0, 1, 0, 0 ; ..
            4, 0, 0, 0, 0, 0, 0 ;          2, 0,-3, 0, 0, 0, 0 ;           ..
            0, 1,-2, 0, 0, 0, 0 ;          2, 0,-1, 2, 0, 0, 0 ; .. // 30
            2,-1,-2, 0, 0, 0, 0 ;          1, 0, 1, 0, 0, 0, 0 ; ..
            2,-2, 0, 0, 0, 0, 0 ;          0, 1, 2, 0, 0, 0, 0 ; ..
            0, 2, 0, 0, 0, 0, 0 ;          2,-2,-1, 0, 0, 0, 0 ; ..
            0, 0, 0,-1, 0, 0, 1 ;          2, 0, 1,-2, 0, 0, 0 ; ..
            2, 0, 0, 2, 0, 0, 0 ;          4,-1,-1, 0, 0, 0, 0 ; .. // 40
            0, 0, 2, 2, 0, 0, 0 ;          3, 0,-1, 0, 0, 0, 0 ; ..
            2, 1, 1, 0, 0, 0, 0 ;          4,-1,-2, 0, 0, 0, 0 ; ..
            0, 2,-1, 0, 0, 0, 0 ;          2, 2,-1, 0, 0, 0, 0 ; ..
            2, 1,-2, 0, 0, 0, 0 ;          2,-1, 0,-2, 0, 0, 0 ; ..
            4, 0, 1, 0, 0, 0, 0 ;          0, 0, 4, 0, 0, 0, 0 ; .. // 50
            4,-1, 0, 0, 0, 0, 0 ;          1, 0,-2, 0, 0, 0, 0 ; ..
            2, 1, 0,-2, 0, 0, 0 ;          0, 0, 2,-2, 0, 0, 0 ; ..
            1, 1, 1, 0, 0, 0, 0 ;          3, 0,-2, 0, 0, 0, 0 ; ..
            4, 0,-3, 0, 0, 0, 0 ;          2,-1, 2, 0, 0, 0, 0 ; ..
            0, 2, 1, 0, 0, 0, 0 ;          0, 0, 0, 0, 0, 1, 0 ; .. // 60
            1, 1,-1, 0, 0, 0, 0 ;          2, 0, 3, 0, 0, 0, 0 ] ;
            
// NB: Reasons to use a loop instead of a vectorized version:
// - less memory used
// - faster for lower number of harmonics (26/13/13)
// - BUT slower for higher number of harmonics (62/66/46)
// --> Maybe a loop dealing with blocks of 5 or 10 harmonics would be faster
            
// Loop over number of harmonics 
lon = zeros(jd);
for k = 1 : nlon

  // Arguments of sine (linear combination of angles : D, M, Mp, F, A1, A2, Lp)
  angle = lonargs(k,1) * D + lonargs(k,2) * M + lonargs(k,3) * Mp + ..
          lonargs(k,4) * F + lonargs(k,5) * A1 + lonargs(k,6) * A2 + lonargs(k,7) * Lp ;
  res = loncoefs(k) * sin(angle);

  // If M (or -M) is in the arguments, correct the coefficient by E
  if (abs(lonargs(k,2)) == 1) 
    res = res .* E;
  end

  // If 2M (or -2M) is in the arguments, correct the coefficient by E^2
  if (abs(lonargs(k,2)) == 2) 
    res = res .* E2;
  end

  lon = lon + res;
end
lon = lon + Lp;

   
   
// ---------------------
// Latitude computation
// ---------------------                                                  

// Magnitude of each sine (rad)
latcoefs = DEG2RAD * 1e-6 * ..
              [5128122; 280602;  277693;  173237;  55413; ..
               46271; 32573;  17198;  9266;  8822; ..
               8216; 4324;  4200; -3359;  2463; ..
              -2235; 2211;  2065; -1870;  1828; .. // 20
              -1794;-1749; -1565; -1491; -1475; ..
              -1410;-1344; -1335;  1107;  1021; ..
               833; 777;  671;  607;  596; ..
               491;-451;  439;  422;  421; .. // 40
               382;-366; -351;  331;  315; ..
               302;-283; -229;  223;  223; ..
              -220;-220; -185;  181; -177; ..
               176; 175;  175;  166; -164; .. // 60
               132; 127; -119; -115;  115; ..
               107];
                           
// Arguments of each sine (linear combination of angles : D, M, Mp, F, A1, A3, Lp)
latargs = [ 0, 0, 0, 1, 0, 0, 0 ;       0, 0, 1, 1, 0, 0, 0 ; ..
            0, 0, 1,-1, 0, 0, 0 ;       2, 0, 0,-1, 0, 0, 0 ; ..
            2, 0,-1, 1, 0, 0, 0 ;       2, 0,-1,-1, 0, 0, 0 ; ..
            2, 0, 0, 1, 0, 0, 0 ;       0, 0, 2, 1, 0, 0, 0 ; ..
            2, 0, 1,-1, 0, 0, 0 ;       0, 0, 2,-1, 0, 0, 0 ; .. // 10
            2,-1, 0,-1, 0, 0, 0 ;       2, 0,-2,-1, 0, 0, 0 ; ..
            2, 0, 1, 1, 0, 0, 0 ;       2, 1, 0,-1, 0, 0, 0 ; ..
            2,-1,-1, 1, 0, 0, 0 ;       0, 0, 0, 0, 0, 0, 1 ; ..
            2,-1, 0, 1, 0, 0, 0 ;       2,-1,-1,-1, 0, 0, 0 ; ..
            0, 1,-1,-1, 0, 0, 0 ;       4, 0,-1,-1, 0, 0, 0 ; .. // 20
            0, 1, 0, 1, 0, 0, 0 ;       0, 0, 0, 3, 0, 0, 0 ; ..
            0, 1,-1, 1, 0, 0, 0 ;       1, 0, 0, 1, 0, 0, 0 ; ..
            0, 1, 1, 1, 0, 0, 0 ;       0, 1, 1,-1, 0, 0, 0 ; ..
            0, 1, 0,-1, 0, 0, 0 ;       1, 0, 0,-1, 0, 0, 0 ; ..
            0, 0, 3, 1, 0, 0, 0 ;       4, 0, 0,-1, 0, 0, 0 ; .. // 30
            4, 0,-1, 1, 0, 0, 0 ;       0, 0, 1,-3, 0, 0, 0 ; ..
            4, 0,-2, 1, 0, 0, 0 ;       2, 0, 0,-3, 0, 0, 0 ; ..
            2, 0, 2,-1, 0, 0, 0 ;       2,-1, 1,-1, 0, 0, 0 ; ..
            2, 0,-2, 1, 0, 0, 0 ;       0, 0, 3,-1, 0, 0, 0 ; ..
            2, 0, 2, 1, 0, 0, 0 ;       2, 0,-3,-1, 0, 0, 0 ; .. // 40
            0, 0, 0, 0, 0, 1, 0 ;       2, 1,-1, 1, 0, 0, 0 ; ..
            2, 1, 0, 1, 0, 0, 0 ;       4, 0, 0, 1, 0, 0, 0 ; ..
            2,-1, 1, 1, 0, 0, 0 ;       2,-2, 0,-1, 0, 0, 0 ; ..
            0, 0, 1, 3, 0, 0, 0 ;       2, 1, 1,-1, 0, 0, 0 ; ..
            1, 1, 0,-1, 0, 0, 0 ;       1, 1, 0, 1, 0, 0, 0 ; .. // 50
            0, 1,-2,-1, 0, 0, 0 ;       2, 1,-1,-1, 0, 0, 0 ; ..
            1, 0, 1, 1, 0, 0, 0 ;       2,-1,-2,-1, 0, 0, 0 ; ..
            0, 1, 2, 1, 0, 0, 0 ;       4, 0,-2,-1, 0, 0, 0 ; ..
            0, 0, 0,-1, 1, 0, 0 ;       0, 0, 0, 1, 1, 0, 0 ; ..
            4,-1,-1,-1, 0, 0, 0 ;       1, 0, 1,-1, 0, 0, 0 ; .. // 60
            4, 0, 1,-1, 0, 0, 0 ;       0, 0,-1, 0, 0, 0, 1 ; ..
            1, 0,-1,-1, 0, 0, 0 ;       0, 0, 1, 0, 0, 0, 1 ; ..
            4,-1, 0,-1, 0, 0, 0 ;       2,-2, 0, 1, 0, 0, 0 ] ;

// Loop over number of harmonics
lat = zeros(jd);
for k = 1 : nlat

  // Arguments of sine (linear combination of angles : D, M, Mp, F, A1, A2, Lp)
  angle = latargs(k,1) * D + latargs(k,2) * M + latargs(k,3) * Mp + ..
          latargs(k,4) * F + latargs(k,5) * A1 + latargs(k,6) * A3 + latargs(k,7) * Lp ;
  res = latcoefs(k) * sin(angle);

  // If M (or -M) is in the arguments, correct the coefficient by E
  if (abs(latargs(k,2)) == 1);  
    res = res .* E;
  end

  // If 2M (or -2M) is in the arguments, correct the coefficient by E^2
  if (abs(latargs(k,2)) == 2);  
    res = res .* E2;
  end

  lat = lat + res;
end


// ---------------------
// Distance computation
// ---------------------

// Magnitude of each cosine (m)
distcoefs = [-20905355;-3699111;-2955968;  -569925; 246158; ..
                 -204586; -170733; -152138;  -129620; 108743; ..
                 104755;    79661;   48888;   -34782;  30824; ..
                  24208;   -23210;  -21636;   -16675;  14403; .. // 20
                 -12831;   -11650;  -10445;    10321;  10056; ..
                  -9884;     8752;   -8379;    -7003;   6322; ..
                   5751;    -4950;   -4421;     4130;  -3958; ..
                   3258;    -3149;    2616;     2354;  -2117; .. // 40
                  -1897;    -1739;   -1571;    -1423;   1165; ..
                  -1117];

// Arguments of each cosine (linear combination of angles : D, M, Mp, F) 
distargs = [  0, 0, 1, 0 ;                       2, 0, -1, 0; .. 
              2, 0, 0, 0 ;                       0, 0, 2, 0; .. 
              2, 0,-2, 0 ;                       2,-1, 0, 0; .. 
              2, 0, 1, 0 ;                       2,-1,-1, 0; .. 
              0, 1,-1, 0 ;                       1, 0, 0, 0; .. // 10
              0, 1, 1, 0 ;                       0, 0, 1,-2; .. 
              0, 1, 0, 0 ;                       4, 0,-1, 0; .. 
              2, 1, 0, 0 ;                       2, 1,-1, 0; .. 
              0, 0, 3, 0 ;                       4, 0,-2, 0; .. 
              1, 1, 0, 0 ;                       2, 0,-3, 0; .. // 20
              2,-1, 1, 0 ;                       4, 0, 0, 0; .. 
              2, 0, 2, 0 ;                       2, 0, 0,-2; .. 
              2,-1,-2, 0 ;                       2,-2, 0, 0; .. 
              2, 0,-1,-2 ;                       1, 0,-1, 0; .. 
              0, 1,-2, 0 ;                       1, 0, 1, 0; .. // 30
              0, 1, 2, 0 ;                       2,-2,-1, 0; .. 
              0, 0, 2,-2 ;                       2, 0, 1,-2; .. 
              4,-1,-1, 0 ;                       3, 0,-1, 0; .. 
              0, 0, 0, 2 ;                       2, 1, 1, 0; .. 
              2, 2,-1, 0 ;                       0, 2,-1, 0; .. // 40
              4,-1,-2, 0 ;                       1, 0,-2, 0; .. 
              4,-1, 0, 0 ;                       4, 0, 1, 0; .. 
              0, 2, 1, 0 ;                       0, 0, 4, 0];

// Loop over number of harmonics
dist = zeros(jd); 
for k = 1 : ndist

  // Arguments of cosine (linear combination of angles : D, M, Mp, F, A1, A2, Lp)
  angle = distargs(k,1) * D + distargs(k,2) * M + distargs(k,3) * Mp + distargs(k,4) * F ;
  res = distcoefs(k) .* cos(angle);

  // If M (or -M) is in the arguments, correct the corresponding coefficient by E
  if (abs(distargs(k,2)) == 1);  
    res = res .* E;
  end

  // If 2M (or -2M) is in the arguments, correct the corresponding coefficient by E^2
  if (abs(distargs(k,2)) == 2);  
    res = res .* E2;
  end

  dist = dist + res;
end
dist = dist + 385000.56e3;

// ---------------------  
// Velocity computation
// ---------------------
if (cvel)
  // Derivative of Moon mean longitude (rad/century)
  Lp_d = (481267.88123421  - 0.0015786 * 2*TT + 3*TT2 / 538841.0 - 4*TT3 / 65194000.0) * DEG2RAD ;

  // Derivate of mean elongation of the Moon (eq. 45.2) (rad/century)
  D_d = (445267.1114034  - 0.0018819 * 2*TT + 3*TT2 / 545868 - 4*TT3 / 113065000) * DEG2RAD ;

  // Derivate of Sun mean anomaly (eq. 45.3) (rad/century)
  M_d = (35999.0502909  - 0.0001536 * 2*TT + 3*TT2 / 24490000) * DEG2RAD ;

  // Derivate of Moon mean anomaly (eq. 45.4) (rad/century)
  Mp_d = (477198.8675055  + 0.0087414 * 2*TT + 3*TT2 / 69699 - 4*TT3 / 14712000) * DEG2RAD ;

  // Derivate of Moon argument of latitude (rad/century)
  F_d = (483202.0175233  - 0.0036539 * 2*TT - 3*TT2 / 3526000 + 4*TT3 / 863310000) * DEG2RAD ;

  // Derivate of A1 : impact of Venus (rad/century)
  A1_d = 131.849 * ones(TT) * DEG2RAD ;

  // Derivate of A2 : impact of Jupiter (rad/century)
  A2_d = 479264.290 * ones(TT) * DEG2RAD ;

  // Derivate of A3 (rad/century) 
  A3_d = 481266.484 * ones(TT) * DEG2RAD ;

  // Derivative of term used to correct harmonics that contain (M) or (2M)
  // (eccentricity of Earth's orbit around the Sun is currently decreasing with time)
  E_d = -0.002516  - 0.0000074 * 2*TT ;
  E2_d = 2*E.*E_d;
  
  
  
  // -----------
  // Longitude :
  // -----------
  
  // Loop over number of harmonics:
  lon_d = zeros(jd);
  for k = 1 : nlon

    // Arguments of sine (linear combination of angles : D, M, Mp, F, A1, A2, Lp)
    angle = lonargs(k,1) * D + lonargs(k,2) * M + lonargs(k,3) * Mp + ..
            lonargs(k,4) * F + lonargs(k,5) * A1 + lonargs(k,6) * A2 + lonargs(k,7) * Lp ;
    angle_d = lonargs(k,1) * D_d + lonargs(k,2) * M_d + lonargs(k,3) * Mp_d + ..
            lonargs(k,4) * F_d + lonargs(k,5) * A1_d + lonargs(k,6) * A2_d + lonargs(k,7) * Lp_d ;
    // (K*sin(a)') = K*a'*cos(a)
    res = loncoefs(k) * angle_d .* cos(angle);

    // If M (or -M) is in the arguments, correct the coefficient by E
    // (K*E*sin(a))' = K*E*a'*cos(a) + K*E'*sin(a)
    if (abs(lonargs(k,2)) == 1)
      res = res .* E + loncoefs(k) * E_d .* sin(angle);
    end
    
    // If 2M (or -2M) is in the arguments, correct the coefficient by E^2
    if (abs(lonargs(k,2)) == 2);  
      res = res .* E2 + loncoefs(k) * E2_d .* sin(angle);
    end
    
    lon_d = lon_d + res;
  end
  lon_d = lon_d + Lp_d;
  
  
  
  // -----------
  // Latitude
  // -----------
  
  // Loop over number of harmonics:
  lat_d = zeros(jd);
  for k = 1 : nlat

    // Arguments of sine (linear combination of angles : D, M, Mp, F, A1, A2, Lp)
    angle = latargs(k,1) * D + latargs(k,2) * M + latargs(k,3) * Mp + ..
            latargs(k,4) * F + latargs(k,5) * A1 + latargs(k,6) * A3 + latargs(k,7) * Lp ;
    angle_d = latargs(k,1) * D_d + latargs(k,2) * M_d + latargs(k,3) * Mp_d + ..
            latargs(k,4) * F_d + latargs(k,5) * A1_d + latargs(k,6) * A3_d + latargs(k,7) * Lp_d ;
    // (K*sin(a)') = K*a'*cos(a)
    res = latcoefs(k) * angle_d .* cos(angle);

    // If M (or -M) is in the arguments, correct the coefficient by E
    // (K*E*sin(a))' = K*E*a'*cos(a) + K*E'*sin(a)
    if (abs(latargs(k,2)) == 1)
      res = res .* E + latcoefs(k) * E_d .* sin(angle);
    end

    // If 2M (or -2M) is in the arguments, correct the coefficient by E^2
    if (abs(latargs(k,2)) == 2);  
      res = res .* E2 + latcoefs(k) * E2_d .* sin(angle);
    end

    lat_d = lat_d + res;
  end
  
  
  
  // -----------
  // Distance
  // -----------
  
  // Loop over number of harmonics:
  dist_d = zeros(jd);
  for k = 1 : ndist

    // Arguments of cosine (linear combination of angles : D, M, Mp, F, A1, A2, Lp)
    angle = distargs(k,1) * D + distargs(k,2) * M + distargs(k,3) * Mp + distargs(k,4) * F ;
    angle_d = distargs(k,1) * D_d + distargs(k,2) * M_d + distargs(k,3) * Mp_d + distargs(k,4) * F_d ;
    // (K*cos(a)') = -K*a'*sin(a)
    res = - distcoefs(k) * angle_d .* sin(angle);

    // If M (or -M) is in the arguments, correct the coefficient by E
    // (K*E*cos(a))' = -K*E*a'*sin(a) + K*E'*cos(a)
    if (abs(distargs(k,2)) == 1)
      res = res .* E + distcoefs(k) * E_d .* cos(angle);
    end
    
    // If 2M (or -2M) is in the arguments, correct the coefficient by E^2
    if (abs(distargs(k,2)) == 2)
      res = res .* E2 + distcoefs(k) * E2_d .* cos(angle);
    end
    
    dist_d = dist_d + res ;
  end
  
  
end
   
   
   
// ------------------   
// Output
// ------------------ 
  
// Position (in EOD frame)
if (~cvel) 
  pos_moon_EOD = CL_co_sph2car([lon;lat;dist]);
  vel_moon_EOD = [];
  
// Position and velocity (in EOD frame)
else
  [pos_moon_EOD,vel_moon_EOD] = CL_co_sph2car([lon;lat;dist],[lon_d;lat_d;dist_d]/(36525*86400));
end
  



endfunction