function [pos,vel] = CL_eph_planet(planet,cjd, tt_tref,model)
// Planet position and velocity vectors in EOD frame
//
// Calling Sequence
// [pos,vel] = CL_eph_planet(planet, cjd, [,tt_tref, model])
//
// Description
// <itemizedlist><listitem>
// <p>Computes planets position and velocity vectors relative to the ecliptic of date (EOD)</p>
// <p></p>
// <p>The origin of the frame is the Sun.</p>
// <p>Available planets are "Mercury", "Venus", "Earth", "Mars", "Jupiter", "Saturn", "Uranus" and "Neptune"</p>
// <p></p>
// <p>The model used can optionaly be modified:</p>
// <p> - model = "med": VSOP87 "med" model is used.</p>
// <p> - model = "high": VSOP87 "high" model is used.</p>
// <p> - model = "full":  VSOP87 "full" model is used.</p>
// <p></p></listitem>
// <listitem>
// <p>Notes:</p>
// <p> - Chosen model greatly affects the time and memory consumed by the function</p>
// <p> - Planet names are case sensitive</p>
// <p></p></listitem>
// <listitem>
// <p>See <link linkend="Reference frames">Reference frames</link> for more details on the EOD frame,
// and <link linkend="Ephemerides">Ephemerides</link> for more details on the models.</p>
// <p></p></listitem>
// </itemizedlist>
//
// Parameters
// planet: (string) Name of planet ("Mercury","Venus","Earth","Mars","Jupiter","Saturn","Uranus" or "Neptune") (1x1)
// cjd: Modified (1950.0) julian day (Time scale: TREF) (1xN)
// tt_tref: (optional) TT-TREF [seconds]. Default is %CL_TT_TREF. (1xN or 1x1)
// model : (optional, string) Name of the model. "med", "high" or "full". Default is "med". (1x1)
// pos : Planet position vector in output frame [m] (3xN)
// vel : Planet velocity vector in output frame [m/s] (3xN)
//
// Bibliography
// Planetary Theories in rectangular and spherical variables: VSOP87 solution. Bretagnon P., Francou G.
// Astronomical Algorithms J.Meeus, Appendix II
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_eph_de405
//
// Examples
// cjd = CL_dat_cal2cjd(2010,02,03,05,35,25);
//
// // Default = standard precision
// pos_mars = CL_eph_planet("Mars", cjd)
//
// // Higher precision
// [pos_uranus_EOD,vel_uranus_EOD] = CL_eph_planet("Uranus", cjd, model = "high");
//
// // Maximum precision
// [pos_saturn_EOD,vel_saturn_EOD] = CL_eph_planet("Saturn", cjd, model = "full");


// Declarations:
  function [nbcoef] = nbcoef_model(model,planet)
  
    // Full model (all the terms are used)
    if (model == "full")
      nbcoef = %inf * ones(3,6);
      
    // High model, number of terms found in J.Meeus book
    elseif (model == "high")
      if (planet == "Mercury")
        nbcoef = [ 38, 16, 10, 8 , 6 , 1 ; .. 
                   14, 11, 9 , 7 , 2 , 0 ; ..
                   13, 8 , 7 , 5 , 0 , 0];
      elseif (planet == "Venus")
        nbcoef = [ 24, 12, 8 , 3 , 3 , 1 ; .. 
                   9 , 4 , 4 , 4 , 1 , 0 ; ..
                   12, 3 , 3 , 1 , 1 , 0];
      elseif (planet == "Earth")
        nbcoef = [ 64, 34, 20, 7 , 3 , 1 ; .. 
                   5 , 2 , 0 , 0 , 0 , 0 ; ..
                   40, 10, 6 , 2 , 1 , 0];
      elseif (planet == "Mars")
        nbcoef = [ 69, 46, 33, 12, 8 , 2 ; .. 
                   16 , 9, 7 , 4 , 3 , 0 ; ..
                   45, 27, 11, 6 , 4 , 0];
      elseif (planet == "Jupiter")
        nbcoef = [ 64, 61, 57, 39, 19, 5 ; .. 
                   26, 22, 14, 9 , 6 , 1 ; ..
                   46, 43, 36, 28, 15, 7];
      elseif (planet == "Saturn")
        nbcoef = [ 90, 79, 63, 48, 27, 12 ; .. 
                   34, 32, 29, 21, 12, 2 ; ..
                   44, 38, 32, 28, 23, 18];
      elseif (planet == "Uranus")
        nbcoef = [ 91, 57, 35, 18, 4 , 0 ; .. 
                   28, 20, 11, 4 , 1 , 0 ; ..
                   59, 35, 18, 10, 2 , 0];
      elseif (planet == "Neptune")
        nbcoef = [ 38, 18, 7 , 4 , 1 , 0 ; .. 
                   17, 13, 6 , 4 , 1 , 0 ; ..
                   32, 15, 5 , 1 , 0 , 0];
      else
        CL__error("Unrecognized planet name")
      end
    
    // Medium model, resuling in 40'' in [2000,2100]
    elseif (model == "med")
      if (planet == "Mercury")
        nbcoef = [ 7, 4, 0, 0, 0, 0 ; ..
                   7, 2, 0, 0, 0, 0 ; ..
                   6, 2, 0, 0, 0, 0  ];
      elseif (planet == "Venus")
        nbcoef = [ 3, 2, 0, 0, 0, 0 ; ..
                   3, 1, 0, 0, 0, 0 ; ..
                   2, 0, 0, 0, 0, 0 ];
      elseif (planet == "Earth")
        nbcoef = [ 3, 2, 0, 0, 0, 0 ; ..
                   0, 0, 0, 0, 0, 0 ; ..
                   3, 0, 0, 0, 0, 0];
      elseif (planet == "Mars")
        nbcoef = [ 13, 3, 0, 0, 0, 0 ; ..
                    4, 1, 0, 0, 0, 0 ; ..
                    4, 1, 0, 0, 0, 0];
      elseif (planet == "Jupiter")
        nbcoef = [ 22, 3, 0, 0, 0, 0 ; ..
                    6, 1, 0, 0, 0, 0 ; ..
                   13, 1, 0, 0, 0, 0];
      elseif (planet == "Saturn")
        nbcoef = [ 26, 5, 0, 0, 0, 0 ; ..
                   10, 1, 0, 0, 0, 0 ; ..
                   15, 1, 0, 0, 0, 0];
      elseif (planet == "Uranus")
        nbcoef = [ 25, 2, 0, 0, 0, 0 ; ..
                    5, 1, 0, 0, 0, 0 ; ..
                   15, 0, 0, 0, 0, 0];
      elseif (planet == "Neptune")
        nbcoef = [ 11, 1, 0, 0, 0, 0 ; ..
                    5, 1, 0, 0, 0, 0 ; ..
                    9, 0, 0, 0, 0, 0];
      else
        CL__error("Unrecognized planet name (case sensitive!)")
      end
      
    else
      CL__error("Unrecognized model")
    end
  endfunction


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

if (~exists("tt_tref", "local")); tt_tref = CL__dataGetEnv("TT_TREF"); end
if (~exists("model","local")); model = "med"; end;

if (typeof(model) <> "string"); CL__error("Invalid input argument ''model''"); end
if (typeof(planet) <> "string"); CL__error("Invalid input argument ''planet''"); end

// Compute velocity only if needed
cvel = %f;
if (argn(1) == 2)
  cvel = %t;
end
  
// Number of coefficients to be used (depends on the planet and the model)
// This function also checks the validity of "model" and "planet"
nbcoef = nbcoef_model(model,planet);

// Conversion from cjd (TREF) to jd (TT)
jd = CL_dat_convert("cjd","jd",cjd+tt_tref/86400);

// Compute position (and velocity if needed) in EOD frame
[pos,vel] = CL__eph_planet_VSOP87(planet, jd, nbcoef, cvel);

endfunction