//  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_ELP(jd, cvel)
// Moon position and velocity in EOD frame (ELP/MPP02 theory)
//
// Calling Sequence
// [pos_moon_EOD,vel_moon_EOD] = CL__eph_moon_ELP(jd [,cvel])
//
// Description
// <itemizedlist><listitem>
// <p>Computes Moon position and velocity relative to the ecliptic of date frame (EOD) using ELP/MPP02 theory.</p>
// <p>The frame is centered on Earth. (pos_moon_EOD is the Earth-Moon vector).</p>
// <p></p>
// <p>Notes:</p>
// <p>- The time scale should be TDB in theory. However, the difference between TDB and TT is not higher than 2ms, 
// inducing a negligible error.</p>
// <p>- The implemented version of ELP/MPP02 is the one fitted to DE405 JPL ephemeris, 
// including the secular terms which ensures that the solution ELP/MPP02 approaches closely the 
// JPL Ephemeris DE406 on a long range (a few seconds over 6 millennia)</p>
// <p>- 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)
// 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
// LUNAR SOLUTION ELP version ELP/MPP02 - Jean CHAPRONT and Gerard FRANCOU Observatoire de Paris - SYRTE department - UMR 8630/CNRS - October 2002
//
// Authors
// CNES - DCT/SB
//
// Examples
// // Julian date, in TT time scale
// jd = 21512210.3;
//
// [pos_moon_EOD,vel_moon_EOD] = CL__eph_moon_ELP(jd)


  function [del, Fpl, cst] = fundArgs(TT,cvel)
  
  del = struct();
  Fpl = struct();
  cst = struct();
  
  // -------------------------------------
  // Fundamental arguments (Moon and EMB)
  // -------------------------------------
  
  // Secular developments of the Moon and Earth-Moon arguments 
  // before corrections (paragraph 3.1, table 1)

  // W1 = mean mean longitude of the Moon [rad]
  W10 = (218 +18/60 + 59.95571/3600) * DEG2RAD;
  W11 = 1732559343.73604  * SEC2RAD; // W11 = nu
  W12 = -6.8084 * SEC2RAD;
  W13 = 0.66040e-2 * SEC2RAD;
  W14 = -0.31690e-4 * SEC2RAD;

  // W2 = mean longitude of the lunar perigee [rad]
  W20 = (83 + 21/60 + 11.67475/3600) * DEG2RAD;
  W21 = 14643420.3171 * SEC2RAD;
  W22 = -38.2631 * SEC2RAD;
  W23 = -0.45047e-1 * SEC2RAD;
  W24 = 0.21301e-3 * SEC2RAD;

  // W3 = mean longitude of the lunar ascending node [rad]
  W30 = (125 + 2/60 + 40.39816/3600) * DEG2RAD;
  W31 = -6967919.5383 * SEC2RAD;
  W32 = 6.3590 * SEC2RAD;
  W33 = 0.76250e-2 * SEC2RAD;
  W34 = -0.35860e-4 * SEC2RAD;

  // T = heliocentric mean longitude of the Earth-Moon barycenter [rad]
  T0 = (100 + 27/60 + 59.13885/3600) * DEG2RAD;
  T1 = 129597742.29300 * SEC2RAD; // T1 = n' = np
  T2 = -0.020200d0 * SEC2RAD;
  T3 = 0.90000e-5 * SEC2RAD;
  T4 = 0.15000e-6 * SEC2RAD;

  // op = mean longitude of the perihelion of the Earth-Moon barycenter [rad]
  op0 = (102 + 56/60 + 14.45766/3600) * DEG2RAD;
  op1 = 1161.24342 * SEC2RAD;
  op2 = 0.529265d0 * SEC2RAD;
  op3 = -0.11814e-3 * SEC2RAD;
  op4 = 0.11379e-4 * SEC2RAD;

  // Corrections to the nominal values (paragraph 4.2, table 3) 
  // and corrections to secular terms to fit DE406 over a long period (paragraph 4.3.3, table 6)
  // [rad/cy^n]
  DW10 = -0.07008 * SEC2RAD;
  DW11 = -0.35106 * SEC2RAD;  // = Dnu
  DW12 = -0.03743 * SEC2RAD;
  DW13 = -0.00018865 * SEC2RAD;
  DW14 = -0.00001024 * SEC2RAD;

  DW20 = 0.20794 * SEC2RAD;
  DW21 = 0.08017 * SEC2RAD;
  DW22 = 0.00470602 * SEC2RAD;
  DW23 = -0.00025213 * SEC2RAD;

  DW30 = -0.07215 * SEC2RAD;
  DW31 = -0.04317 * SEC2RAD;
  DW32 = -0.00261070 * SEC2RAD;
  DW33 = -0.00010712 * SEC2RAD;

  DT0 = -0.00033 * SEC2RAD;
  DT1 = 0.00732 * SEC2RAD; // = Dn' = Dnp

  Dop0 = -0.00749 * SEC2RAD;

  DG = 0.00085 * SEC2RAD;
  DE = -0.00006 * SEC2RAD;
  Dep = 0.00224 * SEC2RAD;

  // ------------
  // Constants
  // ------------
  cst.m = 0.074801329;  // [dimensionless]
  cst.alpha = 0.002571881; // [dimensionless]
  
  // Table 4
  cst.dnu = 0.55604 * SEC2RAD + DW11;
  cst.dG = -0.08066 * SEC2RAD + DG;
  cst.dE = 0.01789 * SEC2RAD + DE;
  cst.dep = -0.12879 * SEC2RAD + Dep;
  cst.dnp = -0.0642 * SEC2RAD + DT1; 
  // NB: in elmpp02.for: they have -0.06424 for cst.dnp.
  
  // Additional corrections (paragraph 4.3.2)
  Bp2 = [0.311079095, -0.004482398, -0.001102485, 0.001056062, 0.000050928];
  Bp3 = [-0.103837907, 0.000668287, -0.001298072, -0.000178028, -0.000037342];

  dW21 = (W21/W11- cst.m*(Bp2(1) + 2/3*cst.alpha/cst.m * Bp2(5)))* DW11 + ...
        (Bp2(1) + 2/3*cst.alpha/cst.m * Bp2(5)) * DT1 + ...
        W11 * (Bp2(2)*DG + Bp2(3)*DE + Bp2(4)*Dep) ;

  dW31 = (W31/W11- cst.m*(Bp3(1) + 2/3*cst.alpha/cst.m * Bp3(5)))* DW11 + ...
        (Bp3(1) + 2/3*cst.alpha/cst.m * Bp3(5)) * DT1 + ...
        W11 * (Bp3(2)*DG + Bp3(3)*DE + Bp3(4)*Dep) ;

  // W1 = mean mean longitude of the Moon [rad]
  W1 = (W10 + DW10) + (W11 + DW11) * TT(2,:) + (W12 + DW12) * TT(3,:) + ... 
       (W13 + DW13) * TT(4,:) + (W14 + DW14) * TT(5,:);

  // W2 = mean longitude of the lunar perigee [rad]
  W2 = (W20 + DW20) + (W21 + DW21 + dW21) * TT(2,:) + (W22 + DW22) * TT(3,:) + ... 
       (W23 + DW23) * TT(4,:) + W24 * TT(5,:);

  // W3 = mean longitude of the lunar ascending node [rad]
  W3 = (W30 + DW30) + (W31 + DW31 + dW31) * TT(2,:) + (W32 + DW32) * TT(3,:) + ... 
       (W33 + DW33) * TT(4,:) + W34 * TT(5,:);

  // T = heliocentric mean longitude of the Earth-Moon barycenter [rad]
  T = (T0 + DT0) + (T1 + DT1) * TT(2,:) + T2 * TT(3,:) + T3 * TT(4,:) + T4 * TT(5,:);

  // op = the mean longitude of the perihelion of the Earth-Moon barycenter [rad]
  op = (op0 + Dop0) + op1 * TT(2,:) + op2 * TT(3,:) + op3 * TT(4,:) + op4 * TT(5,:);

  // nu will be needed when computing the series to apply the corrections to A
  cst.nu = (W11 + DW11);
  
  // ---------------------
  // Delaunay arguments [rad]
  // ---------------------
  del.D = W1 - T + %pi;
  del.F = W1 -W3;
  del.l = W1 -W2;
  del.lp = T - op;
  
  // NB: W1 is not a Delaunay argument but stored in the structure for convenience
  //     (will be needed to later to compute longitude)
  del.W1 = W1;
  
  // -------------------------------
  // Planetary arguments (VSOP2000) [rad]
  // -------------------------------
  // Mercury, Venus, Earth, Mars, Jupiter, Saturn, Uranus, Neptune
  Fpl.Me = (252 + 15/60 +  3.216919/3600) * DEG2RAD + 538101628.66888 * SEC2RAD * TT(2,:);
  Fpl.V = (181 + 58/60 + 44.758419/3600) * DEG2RAD + 210664136.45777 * SEC2RAD * TT(2,:);
  Fpl.T = T ; 
  // NB : In elpmpp02.for they use the value below. In the paper (paragraph 2.2.2) they explicitely say to use T.
  // Fpl.T = (100 + 27/60 + 59.138850/3600) * DEG2RAD + 129597742.29300 * SEC2RAD * TT(2,:); 
  Fpl.Ma = (355 + 26/60 +  3.642778/3600) * DEG2RAD + 68905077.65936 * SEC2RAD * TT(2,:);
  Fpl.J = ( 34 + 21/60 +  5.379392/3600) * DEG2RAD + 10925660.57335 * SEC2RAD * TT(2,:); 
  Fpl.S = ( 50 +  4/60 + 38.902495/3600) * DEG2RAD + 4399609.33632 * SEC2RAD * TT(2,:);
  Fpl.U = (314 +  3/60 +  4.354234/3600) * DEG2RAD + 1542482.57845 * SEC2RAD * TT(2,:);
  Fpl.N = (304 + 20/60 + 56.808371/3600) * DEG2RAD + 786547.89700 * SEC2RAD * TT(2,:);
  
  // Zeta : Mean longitude W1 + Rate of the precession
  cst.dprec = -0.29965 * SEC2RAD;
  Dprec = 5029.0966 * SEC2RAD + cst.dprec ;
  Fpl.zeta = (W10 + DW10) + (W11 + DW11 + Dprec) * TT(2,:) + (W12 + DW12) * TT(3,:) + ... 
            (W13 + DW13) * TT(4,:) + (W14 + DW14) * TT(5,:);
     
  // Time derivatives (if needed)
  if (cvel)
    // rad/century
    W1dot = (W11 + DW11) + (W12 + DW12) * 2*TT(2,:) + ... 
            (W13 + DW13) * 3*TT(3,:) + (W14 + DW14) * 4*TT(4,:);
    W2dot = (W21 + DW21 + dW21) + (W22 + DW22) * 2*TT(2,:) + ... 
            (W23 + DW23) * 3*TT(3,:) + W24 * 4*TT(4,:);
    W3dot = (W31 + DW31 + dW31) + (W32 + DW32) * 2*TT(2,:) + ... 
            (W33 + DW33) * 3*TT(3,:) + W34 * 4*TT(4,:);
    Tdot = (T1 + DT1) + T2 * 2*TT(2,:) + T3 * 3*TT(3,:) + T4 * 4*TT(4,:);
    opdot = op1 + op2 * 2*TT(2,:) + op3 * 3*TT(3,:) + op4 * 4*TT(4,:);
    
    // rad/s
    del.Ddot = (W1dot - Tdot) * CEN2SEC;
    del.Fdot = (W1dot -W3dot) * CEN2SEC;
    del.ldot = (W1dot -W2dot) * CEN2SEC;
    del.lpdot = (Tdot - opdot) * CEN2SEC;
    del.W1dot = W1dot * CEN2SEC;
    
    // rad/s
    Fpl.Medot = 538101628.66888 * SEC2RAD * TT(1,:) * CEN2SEC;
    Fpl.Vdot = 210664136.45777 * SEC2RAD * TT(1,:) * CEN2SEC;
    Fpl.Tdot = Tdot * CEN2SEC; 
    // Fpl.Tdot = 129597742.29300 * SEC2RAD * CEN2SEC; // In elpmpp02.for they use this value. Not justified!
    Fpl.Madot = 68905077.65936 * SEC2RAD * TT(1,:) * CEN2SEC;
    Fpl.Jdot = 10925660.57335 * SEC2RAD * TT(1,:) * CEN2SEC; 
    Fpl.Sdot = 4399609.33632 * SEC2RAD * TT(1,:) * CEN2SEC;
    Fpl.Udot = 1542482.57845 * SEC2RAD * TT(1,:) * CEN2SEC;
    Fpl.Ndot = 786547.89700 * SEC2RAD * TT(1,:) * CEN2SEC;
    Fpl.zetadot = ((W11 + DW11 + Dprec) + (W12 + DW12) * 2*TT(2,:) + ...
                  (W13 + DW13) * 3*TT(3,:) + (W14 + DW14) * 4*TT(4,:)) * CEN2SEC;
  end
  
  endfunction


  function L = decoupe(P,Pmax)
  // Utility function that returns a list of index from 1 to P divided into sub-groups
  // Each member of the list (i.e each sub-group) contains Pmax index (at max)
  // The last member of the list contains from 1 to Pmax index.
  // NB: if P is [], the returned list is empty
  // NB2: P and Pmax must be (1x1) (not checked)
  
    nb = ceil(P/Pmax);
    L = list();
    
    for k = 1 : nb
      L(k) = [1+(k-1)*Pmax : min(P,k*Pmax)]
    end
  
  endfunction
  
  
  function [pos_sum, vel_sum] = series_main_lon_lat(s,del,cst,cvel)
  // Computes series for the main problem (longitude or latitude)
  //
  // s: structure containing fields A,i1,i2,i3,i4,B1,B2,B3,B4,B5,B6
  //    A expressed in arcsec or km. i1,i2,i3,i4 are dimensionless integers (int8 Scilab variable)
  //    each field is of size (Px1)
  //
  // del: Structure containing Delaunay arguments (D,F,l,lp) (1xN)
  // cst: Structure containing constants
  // pos_sum: sum( Acor * sin(i1*D + i2*F + i3*l + i4*lp)) (1xN)
  // vel_sum: derivative of pos_sum with respect to T (1xN)
  
  P = size(s.A,1);
  N = size(del.D,2);
  pos_sum = [];
  vel_sum = [];
  
  // Chapter 4.3.1: computation of Acor
  // NB: Acor is (Px1) (same as s.A)
  Acor = s.A + (s.B1 + 2/3 * cst.alpha/cst.m * s.B5) * (cst.dnp/cst.nu - cst.m * cst.dnu/cst.nu) +  ...
               (s.B2*cst.dG + s.B3*cst.dE + s.B4*cst.dep);
  
  // Semi-vectorized version (Pmax rows are treated at once)
  // NB: Read additional comments about semi-vectorized versions in function series_pert.
  Pmax = 50;
  for p = decoupe(P,Pmax)
    phi = double(s.i1(p))*del.D + double(s.i2(p))*del.F + double(s.i3(p))*del.l + double(s.i4(p))*del.lp;
    pos_sum  = pos_sum + sum( (Acor(p)*ones(1,N)) .* sin(phi), "r");
    
    if (cvel)
      phidot = double(s.i1(p))*del.Ddot + double(s.i2(p))*del.Fdot + double(s.i3(p))*del.ldot + double(s.i4(p))*del.lpdot;
      vel_sum = vel_sum + sum( (Acor(p)*ones(1,N)) .* phidot .* cos(phi), "r");             
    end
  end
  
  endfunction
  
  function [pos_sum, vel_sum] = series_main_dist(s,del,cst,cvel)
  // Computes series for the main problem (distance)
  //
  // s: structure containing fields A,i1,i2,i3,i4,B1,B2,B3,B4,B5,B6
  //    A expressed in arcsec or km. i1,i2,i3,i4 are dimensionless integers (int8 Scilab variable)
  //    each field is of size (Px1)
  //
  // del: Structure containing Delaunay arguments (D,F,l,lp) (1xN)
  // cst: Structure containing constants
  // pos_sum: sum( Acor * cos(i1*D + i2*F + i3*l + i4*lp)) (1xN)
  // vel_sum: derivative of pos_sum with respect to T (1xN)
  
  P = size(s.A,1);
  N = size(del.D,2);
  pos_sum = [];
  vel_sum = [];
  
  // Chapter 4.3.1: computation of Acor:
  // NB: Acor is (Px1) (same as s.A)
  Acor = s.A + (s.B1 + 2/3 * cst.alpha/cst.m * s.B5 + 2/3 * s.A/cst.m) * (cst.dnp/cst.nu - cst.m * cst.dnu/cst.nu) +  ...
               (s.B2*cst.dG + s.B3*cst.dE + s.B4*cst.dep);
  
  // Semi-vectorized version (Pmax rows are treated at once)
  // NB: Read additional comments about semi-vectorized versions in function series_pert.
  Pmax = 50;
  for p = decoupe(P,Pmax)
    phi = double(s.i1(p))*del.D + double(s.i2(p))*del.F + double(s.i3(p))*del.l + double(s.i4(p))*del.lp;
    pos_sum  = pos_sum + sum( (Acor(p)*ones(1,N)) .* cos(phi),"r");
    
    if (cvel)
      phidot = double(s.i1(p))*del.Ddot + double(s.i2(p))*del.Fdot + double(s.i3(p))*del.ldot + double(s.i4(p))*del.lpdot;
      vel_sum = vel_sum - sum( (Acor(p)*ones(1,N)) .* phidot .* sin(phi), "r");             
    end
  end
  
  endfunction

  function [pos_sum, vel_sum] = series_pert(s,del,cst,TT,Fpl,cvel)
  // Computes series for the perturbation problem 
  // 
  // s: list containing the series
  //    First element of the list corresponds to the series of T^0
  //    Second element of the list corresponds to the series of T^1
  //    ...
  //    Fourth element of the list corresponds to the series of T^3
  //
  // Each element of the list is a structure containing the fields S,C,i1,i2,i3,i4,i5,i6,i7,i8,i9,i10,i11,i12,i13
  //    S,C are expressed in arcsec or km. ix are dimensionless integers (int8 Scilab variable)
  //    each field is of size (Px1)
  //
  // del: Structure containing Delaunay arguments (D,F,l,lp) (1xN)
  // cst: Structure containing constants.
  // TT: T^0; T^1; ... T^5 (5xN) with T = time in centuries from epoch 2000 12h (TT time scale)
  // Fpl: Planetary arguments (8xN)
  // pos_sum: sum( s * sin(phi) + C * cos(phi))*t^i (1xN)
  // vel_sum: derivative of pos_sum with respect to T (1xN)
  
  pos_sum = [];
  vel_sum = [];
  N = size(del.D,2);
  
  // NB : The fully vectorized version is avoided because of stacksize errors.
  //      (Even when setting stacksize("max"), errors appears for a number of dates (N) as low as 1000)
  //
  // NB2 : Performance wise, the computation is slower when computing for only a few dates (0.2 sec instead of 0.12 for one date)
  //       But the computation is faster when computing for dozens or hundreds of dates (5.4 sec instead of 8 sec for 800 dates)
  //       And the function can be used for 50000 dates and higher without any stacksize errors (although it takes 450 sec!)
  
  // Semi-vectorized version (Pmax rows are treated at once)
  // NB: Setting Pmax to a high value (1e8 for example) will actually correspond to the fully vectorized version.
  Pmax = 50;
  for k = 1 : 4
    P = size(s(k).S,1);
    for p = decoupe(P,Pmax)
      phi = double(s(k).i1(p))*del.D + double(s(k).i2(p))*del.F + double(s(k).i3(p))*del.l + double(s(k).i4(p))*del.lp + double(s(k).i5(p))*Fpl.Me + ...
            double(s(k).i6(p))*Fpl.V + double(s(k).i7(p))*Fpl.T + double(s(k).i8(p))*Fpl.Ma + double(s(k).i9(p))*Fpl.J + ...
            double(s(k).i10(p))*Fpl.S + double(s(k).i11(p))*Fpl.U + double(s(k).i12(p))*Fpl.N + double(s(k).i13(p))*Fpl.zeta;
      
      pos_sum = pos_sum + sum( s(k).S(p)*TT(k,:) .* sin(phi) + s(k).C(p)*TT(k,:) .* cos(phi), "r");
    
      if (cvel)
        phidot = double(s(k).i1(p))*del.Ddot + double(s(k).i2(p))*del.Fdot + double(s(k).i3(p))*del.ldot + double(s(k).i4(p))*del.lpdot + double(s(k).i5(p))*Fpl.Medot + ...
                 double(s(k).i6(p))*Fpl.Vdot + double(s(k).i7(p))*Fpl.Tdot + double(s(k).i8(p))*Fpl.Madot + double(s(k).i9(p))*Fpl.Jdot + ...
                 double(s(k).i10(p))*Fpl.Sdot + double(s(k).i11(p))*Fpl.Udot + double(s(k).i12(p))*Fpl.Ndot + double(s(k).i13(p))*Fpl.zetadot;
        vel_sum = vel_sum + sum( s(k).S(p)*TT(k,:) .* phidot .* cos(phi) - s(k).C(p)*TT(k,:) .* phidot .* sin(phi), "r");  
        if (k >= 2)
          vel_sum = vel_sum + sum( (s(k).S(p)*(k-1)*TT(k-1,:)*CEN2SEC) .* sin(phi) + (s(k).C(p)*(k-1)*TT(k-1,:)*CEN2SEC) .* cos(phi), "r");  
        end
      end
    end
  end
  
  
  endfunction


  function [pos_moon_EOD,vel_moon_EOD] = convFrame(res, del, TT, cvel);
  // Conversion to the mean ecliptic of date frame (EOD)
  // Conversion to rectangular coordinates
  // See chapter 5.1
  
  // pA is the accumulated precession between J2000 and the date derived from (Laskar, 1986) [rad]
  pA = (5029.0966 * TT(2,:) + 1.1120* TT(3,:) + 0.000077* TT(4,:) - 0.00002353* TT(5,:)) * SEC2RAD; 
  
  // Dp is a correction to the precession constant from (Herring et al., 2002) [rad/century]
  Dp = -0.29965 * SEC2RAD;
  
  // Longitude in EOD frame [rad]
  Vd = (res.lon_main + res.lon_pert) * SEC2RAD + del.W1 + pA + Dp * TT(2,:);
  
  // Latitude in EOD frame [rad]
  Ud = (res.lat_main + res.lat_pert) * SEC2RAD;
  
  // Distance [m]
  ra0 = 384747.961370173 / 384747.980674318; // a0(DE405)/a0(ELP)
  r = (res.dist_main + res.dist_pert)*1000 * ra0;
  
  if (~cvel)
    pos_moon_EOD = CL_co_sph2car([Vd;Ud;r]);
    vel_moon_EOD = []; 
    
  else
    // rad/century
    pAdot = (5029.0966 + 1.1120* 2*TT(2,:) + 0.000077* 3*TT(3,:) - 0.00002353* 4*TT(4,:)) * SEC2RAD;
    // rad/s
    Vddot = (res.londot_main + res.londot_pert) * SEC2RAD + (pAdot + Dp) * CEN2SEC + del.W1dot;
    Uddot = (res.latdot_main + res.latdot_pert) * SEC2RAD;
    // m/s
    rdot = (res.distdot_main + res.distdot_pert)*1000 * ra0;
    [pos_moon_EOD,vel_moon_EOD] = CL_co_sph2car([Vd;Ud;r], [Vddot;Uddot;rdot]);
  end

  endfunction

  
// Declarations:
DEG2RAD = %pi / 180;

// Code:
if (~exists("cvel","local")); cvel = %t; end;
if (argn(2) <> 1 & argn(2) <> 2)
  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

// Arcsecond to radian
SEC2RAD = 1/3600 * DEG2RAD;
// Centuries to seconds
CEN2SEC = 1 / (36525*86400);

// Time in centuries from epoch 2000 12h (TT time scale)
TT = (jd - 2451545.0) / 36525.0;
TT = [ ones(TT); TT ; TT.^2 ; TT.^3 ; TT.^4] ;

// Structure to store constants
res = struct();

// Note : 
// The corrections used are the ones to fit DE405 
// (including additive corrections to secular terms to fit DE406 over a long period)
// (See paragraph 4)

// -------------------------------
// Fundamental arguments (Delaunay arguments, Planetary arguments)
// + constants
// -------------------------------
[del, Fpl, cst] = fundArgs(TT,cvel);
     
// -------------------------------
// Load series
// -------------------------------   
fpath = fullfile(CL_home(), "data", "ephem", "elp_mpp02.dat");
if (~isfile(fpath)); CL__error("File " + fpath + "does not exist"); end;
// Loads a variable named "s"
load(fpath, "s");
     
// -------------------------------
// Computation of main problem and perturbations
// in the natural coordinate system of the lunar theories ELP
// (inertial mean ecliptic of date and departure point gamma' 2000.
// [arcsec, km and arcsec/s, km/s]
// -------------------------------
res = struct();

[res.lon_main, res.londot_main] = series_main_lon_lat(s.main.lon,del,cst,cvel);
[res.lat_main, res.latdot_main] = series_main_lon_lat(s.main.lat,del,cst,cvel);
[res.dist_main, res.distdot_main] = series_main_dist(s.main.dist,del,cst,cvel);

[res.lon_pert, res.londot_pert] = series_pert(s.pert.lon,del,cst,TT,Fpl,cvel);
[res.lat_pert, res.latdot_pert] = series_pert(s.pert.lat,del,cst,TT,Fpl,cvel);
[res.dist_pert, res.distdot_pert] = series_pert(s.pert.dist,del,cst,TT,Fpl,cvel);

// -------------------------------
// Conversion to the mean ecliptic of date frame (EOD)
// Conversion to rectangular coordinates
// -------------------------------
[pos_moon_EOD,vel_moon_EOD] = convFrame(res, del, TT, cvel);


endfunction