//  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 [pars] = CL_mod_planetsG50ecl(cjd,planet_name)
// Planet positions in Gamma50 ecliptic reference frame - DEPRECATED
//
// Calling Sequence
// pars = CL_mod_planetsG50ecl(cjd,planet_name)
//
// Description
// <itemizedlist><listitem>
// <p>This function is deprecated. </p>
// <p>Replacement function: <link linkend="CL_eph_planet">CL_eph_planet</link></p>
// <p></p></listitem>
// <listitem>
// <p>Computes the Keplerian orbit parameters of planets in the Gamma 50 ecliptic 
// reference frame, using a simple model. </p>
// </listitem>
// </itemizedlist>
//
// Parameters
// cjd: Modified (1950.0) julian day (1xN)
// planet_name: (string) Name of planet (Mercury,Venus,Earth,Mars,Jupiter,Saturn,Uranus,Neptune or Pluto)
// pars: Planet's keplerian parameters at time cjd [sma,ecc,inc,pom,gom,M] [6xN]
//
// See also
// CL_mod_moonSunG50
//
// Authors
// CNES - DCT/SB
//
// Examples
// cjd = CL_dat_cal2cjd(2016,1,1,0,0,0);
// dates = cjd : 1 : cjd+365*5;
// kep_Earth = CL_mod_planetsG50ecl(dates,'Earth');
// kep_Mars = CL_mod_planetsG50ecl(dates,'Mars');


// 17/10/2008 ThM Correction bug sur le cas Terre + aide


// Declarations:

// Code:
CL__warnDeprecated(); // deprecated function

ua = CL__dataGetEnv("au", internal=%t);

Ncjd = size(cjd,2);

rdda = 57.295779513082;
ua = ua/1000*ones(1,Ncjd); // conversion en km
tu1 = cjd  +  18262.5
te = tu1
T = tu1 / 36525
t2 = T.* T
t3 = t2.* T
tp = cjd./ 365250
preces = (50267.525  +  111.16.* tp).* tp./ 3600

select planet_name

//MERCURY
  case 'Mercury'
    xl = 178 + 10/60 + (44.68 + 538106654.80 .* T + 1.084 .* t2) ./ 3600
    xlper = 75 + 53/60 + (58.91 + 5599.76 .* T + 1.061 .* t2) ./ 3600
    gom = 47 + 8/60 + (45.40 + 4266.75 .* T + 0.626 .* t2) ./ 3600
    ecc = 0.20561421 + 2.046e-5 .* T - 3e-8 .* t2
    inc = 7 + (10.37 + 6.699 .* T - 0.066 .* t2) ./ 3600
    sma = ua .* 0.38709860

 //VENUS
  case 'Venus'
    xl = 342 + 46 ./ 60 + (1.39 + 210669162.88 .* T + 1.1148 .* t2) ./ 3600
    xlper = 130 + 9 ./ 60 + (49.8 + 5068.93 .* T - 3.515 .* t2) ./ 3600
    gom = 75 + 46 ./ 60 + (46.73 + 3239.46 .* T + 1.476 .* t2) ./ 3600
    ecc = 0.00682069 - 4.774e-5 .* T + 9.1e-8 .* t2
    inc = 3 + 23 ./ 60 + (37.07 + 3.621 .* T - 0.0035 .* t2) ./ 3600
    sma = ua .* 0.72333162

 //EARTH
  case 'Earth'
    xl = 99 + 41 ./ 60 + (48.04 + 129602768.13 .* T + 1.089 .* t2) ./ 3600
    xlper = 101 + 13 ./ 60 + (15.0 + 6189.03 .* T + 1.63 .* t2 + 0.012 .* t3) ./ 3600
    gom = 0.* T
    ecc = 0.01675104 - 4.18e-5 .* T - 1.26e-7 .* t2
    inc = 0.* T
    sma = ua .* 1.00000023

 //MARS
  case 'Mars'
    xl = 293 + 44 ./ 60 + (51.46 + 68910117.33 .* T + 1.1184 .* t2) ./ 3600
    xlper = 334 + 13 ./ 60 + (5.53 + 6626.73 .* T + 0.4675 .* t2 - 0.0043 .* t3) ./ 3600
    gom = 48 + 47 ./ 60 + (11.19 + 2775.57 .* T - 0.005 .* t2 - 0.0192 .* t3) ./ 3600
    ecc = 0.09331290 + 9.2064e-5 .* T - 7.7e-8 .* t2
    inc = 1 + 51 ./ 60 + (1.20 - 2.430 .* T + 0.0454 .* t2) ./ 3600
    sma = ua .* 1.52368840

 //JUPITER
  case 'Jupiter'
    xl = 238 + 2 ./ 60 + (57.32 + 10930687.148 .* T + 1.20486 .* t2 - 0.005936 .* t3) ./ 3600
    xlper = 12 + 43 ./ 60 + (15.34 + 5795.862 .* T + 3.80258 .* t2 - 0.01236 .* t3) ./ 3600
    gom = 99 + 26 ./ 60 + (36.19 + 3637.908 .* T + 1.2680 .* t2 - 0.03064 .* t3) ./ 3600
    ecc = 0.04833475 + 1.64180e-4 .* T - 4.676e-7 .* t2 - 1.7e-9 .* t3
    inc = 1 + 18 ./ 60 + (31.45 - 20.506 .* T + 0.014 .* t2) ./ 3600
    sma = ua .* 5.202561

 //SATURN
  case 'Saturn'
    xl = 266 + 33 ./ 60 + (51.76 + 4404635.581 .* T + 1.16835 .* t2 - 0.0021 .* t3) ./ 3600
    xlper = 91 + 5 ./ 60 + (53.38 + 7050.297 .* T + 2.9749 .* t2 + 0.0166 .* t3) ./ 3600
    gom = 112 + 47 ./ 60 + (25.40 + 3143.5025 .* T - 0.54785 .* t2 - 0.0191 .* t3) ./ 3600
    ecc = 0.05589232 - 3.4550e-4 .* T - 7.28e-7 .* t2 + 7.4e-10 .* t3
    inc = 2 + 29 ./ 60 + (33.07 - 14.108 .* T - 0.05576 .* t2 + 0.00016 .* t3) ./ 3600
    sma = ua .* 9.554747

 //URANUS
  case 'Uranus'
    xl = 244 + 11 ./ 60 + (50.89 + 1547508.765 .* T + 1.13774 .* t2 - 0.002176 .* t3) ./ 3600
    xlper = 171 + 32 ./ 60 + (55.14 + 5343.958 .* T + 0.8539 .* t2 - 0.00218 .* t3) ./ 3600
    gom = 73 + 28 ./ 60 + (37.55 + 1795.204 .* T + 4.722 .* t2) ./ 3600
    ecc = 0.0463444 - 2.658e-5 .* T + 7.7e-8 .* t2
    inc = 46 ./ 60 + (20.87 + 2.251 .* T + 0.1422 .* t2) ./ 3600
    sma = ua .* 19.21814

 //NEPTUNE
  case 'Neptune'
    xl = 84 + 27 ./ 60 + (28.78 + 791589.291 .* T + 1.15374 .* t2 - 0.002176 .* t3) ./ 3600
    xlper = 46 + 43 ./ 60. + (38.37 + 5128.468 .* T + 1.40694 .* t2 - 0.002176 .* t3) ./ 3600
    gom = 130 + 40 ./ 60 + (52.89 + 3956.166 .* T + 0.89952 .* t2 - 0.016984 .* t3) ./ 3600
    ecc = 0.00899704 + 6.330e-6 .* T - 2e-9 .* t2
    inc = 1 + 46 ./ 60 + (45.27 - 34.357 .* T - 0.0328 .* t2) ./ 3600
    sma = ua .* 30.10957

 //Pluto
  case 'Pluto'
    sma = ua .* 39.438712
    ecc = 0.250236
    inc = 17.17047
    gom = 109.73994
    pom = 114.27032
    aM = 317.48609 + 0.0039794275 .* (te - 37200.5)

  else
    CL__error('unknown planet');

end

 //FINAL CALCULATION
xl = CL_rMod(xl,360.0)
xlper = CL_rMod(xlper,360.0)
pom = CL_rMod(xlper - gom,360.0)
aM = xl - xlper

aM = CL_rMod(aM,360.0)
gom = gom - preces
gom = CL_rMod(gom,360.0)

pars = [sma*1000; ecc ; CL_deg2rad([inc;pom;gom;aM])]; // conversion en m

endfunction