//  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 [mean_oe, cjd] = CL_tle_evalAvgElem(tle, type_oe, frame, ut1_tref, tt_tref, mu)
// Evaluates TLE mean orbital elements at epoch
//
// Calling Sequence
// [mean_oe, cjd] = CL_tle_evalAvgElem(tle, type_oe [, frame, ut1_tref, tt_tref])
//
// Description
// <itemizedlist>
// <listitem>
// <p>Evaluates the mean orbital elements at epoch numerically. The processing consists in removing 
// short period terms from the osculating elements (in the specified frame) over a short period 
// of time (a few orbits) The mean elements obtained should be comparable with the mean elements 
// defined by the "eckhech" or "lydlp" models.  </p>
// <p></p>
// <p>The orbital elements can be returned in various forms: Keplerian elements, position and 
// velocity, etc... </p>
// <p></p>
// <p>In a standard situation, the orbital elements are converted to the specified reference 
// <b>frame</b>. The epoch (internally in UTC) is converted to the TREF time scale. Note that if 
// the frame is not a true equator frame, the transformation may not be as accurate as should be
// as the frame transformation arguments that can be provided are limited to ut1_tref and tt_tref. </p>
// <p></p>
// <p>It is also possible to give <b>frame</b> a special value: "native". It means that the 
// output reference frame AND time scale are those natively used for TLEs: TEME and UTC 
// respectively. No conversion is performed. </p> 
// </listitem>
// </itemizedlist>
// 
// Parameters
// tle: TLE structure (size N)
// type_oe: (string) Type of output orbital elements: "kep", "cir", "cireq", "equin" or "pv" (1x1)
// frame: (string, optional) Output frame, see above for details. Default is "ECI". 
// ut1_tref: (optional) UT1-TREF [seconds]. Default is %CL_UT1_TREF (1xN or 1x1)
// tt_tref: (optional) TT-TREF [seconds]. Default is %CL_TT_TREF (1xN or 1x1)
// mean_oe: Mean orbital elements (of type given by "type_oe") relative to the specified frame [m, m/s, rad] (6xN)
// cjd: TLE epoch (modified julian date from 1950.0, TREF time scale except if frame == "native") (1xN)
//
// Authors
// CNES - DCT/SB/MS
//
// See also
// CL_tle_genEphem
//
// Examples
// fnames = CL_path("tle_examples.txt", CL_home(), "all");
// tle = CL_tle_load(fnames); 
// // all "SPOT" satellites
// I = grep(tle.desc, "SPOT"); 
// tle.desc(I)
// // cjd in TREF, mean elements relative to ECI
// [mean_cir, cjd] = CL_tle_evalAvgElem(tle(I), "cir") 


// Declarations:

// eval mean orbital elements for one TLE
// mean_kep: relative to TEME 
function [mean_kep] = tle_evalAvgElem(tle, type_oe, frame, ut1_tref, tt_tref, mu)

  // mean motion, eccentricity and period
  n = tle.n; 
  ecc = tle.ecc; 
  T = 2*%pi ./ n; 

  tle = CL_tle_set(tle, "bstar", 0); 
  cjd0 = tle.epoch_cjd; 

  nborb = 4; 
  nbpts_orb = 20; 
  nb = nbpts_orb * nborb + 3; 

  // mean anomaly (approx)
  anm = linspace(-nborb*%pi, nborb*%pi, nb);
  cjd = cjd0 + (T / 86400) * (anm / (2*%pi));  

  [p, v] = CL_tle_genEphem(tle, cjd, "native"); 

  equin = CL__oe_car2equin(p, v, mu); 
  equin(6,:) = CL_unMod(equin(6,:), 2*%pi); 

  x = anm; 
  A = [ones(x); x; cos(x); sin(x); cos(2*x); sin(2*x); cos(3*x); sin(3*x)]'; 
  X = lsq(A, equin'); 
  
  if (%f)
  // DEBUG
  for (k = 1:6)
  scf(); 
  plot(x, equin(k,:)); 
  plot(x, X(1,k) + X(2,k)*x, "r"); 
  plot(0, X(1,k), "ro"); 

  xtitle(string(k)); 
  end
  end
  
  equin0 = X(1,:)'; 
  mean_kep = CL__oe_equin2kep(equin0); 
  
endfunction


// Code
if (~exists("frame", "local")); frame="ECI"; end
if (~exists("ut1_tref", "local")); ut1_tref = CL__dataGetEnv("UT1_TREF"); end 
if (~exists("tt_tref", "local")); tt_tref = CL__dataGetEnv("TT_TREF"); end 

// Orbital elements types
Types_oe = [ "kep", "cir", "cireq", "equin", "pv" ]; 

// check validity of TLE struct
CL__tle_checkValid(tle); 

if (typeof(frame) <> "string" | size(frame, "*") <> 1)
  CL__error("Invalid argument type or size for frame"); 
end

if (typeof(type_oe) <> "string" | size(type_oe, "*") <> 1)
  CL__error("Invalid argument type or size for type_oe"); 
end

if (find(Types_oe == type_oe) == [])
  CL__error("Invalid orbital element type"); 
end

cst = CL_tle_getConst(); 
N = size(tle); 

mean_kep = %nan * ones(6,N); 
cjd = %nan * ones(1,N); 

// filters valid TLEs
// generates Keplerian elements
I = find(tle.n >= 0 & tle.ecc >= 0 & tle.status == 0); 
for k = I
  [mean_kep(:,k)] = tle_evalAvgElem(tle(k), "kep", frame, ut1_tref, tt_tref, cst.mu); 
end

// change time scale and frame
// done here to improve efficiency (only one call to CL_fr_convert) 
if (frame == "native")
  cjd = tle.epoch_cjd; 
  mean_oe = CL_oe_convert("kep", type_oe, mean_kep, cst.mu); 
else
  cjd = CL_dat_scaleConvert("UTC", "TREF", tle.epoch_cjd, ut1_tref, tt_tref); 
  pv = CL_oe_convert("kep", "pv", mean_kep, cst.mu); 
  [pv(1:3,:), pv(4:6,:)] = CL_fr_convert("TEME", frame, cjd, pv(1:3,:), pv(4:6,:), ut1_tref=ut1_tref, tt_tref=tt_tref); 
  mean_oe = CL_oe_convert("pv", type_oe, pv, cst.mu); 
end
     

endfunction