//  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'.

// Converts mean motion to sma or sma to mean motion
// mean motion to sma: "unkozai"  (cmd = "n2a")
// sma to mean motion: "kozai"    (cmd = "a2n") 
// x, ecc, inc: same size (1xN) 

// TEST: 
// whichconst = "wgs72";
// sma = 7200.e3; 
// ecc = 0.1;
// inc = 1;
// n = CL__tle_sma("a2n", sma, ecc, inc, whichconst)
// sma2 = CL__tle_sma("n2a", n, ecc, inc, whichconst)

function [y] = CL__tle_sma(cmd, x, ecc, inc, whichconst)

// Unkozai mean motion and compute semi major axis 
// => n and sma as computed follow "Brouwer" convention
// Reference: 
// A HISTORY OF ANALYTICAL ORBIT MODELING IN THE UNITED STATES SPACE SURVEILLANCE SYSTEM
// Felix R. Hoots, Paul W. Schumacher, Jr. and Robert A. Glover
// 
// Note: The results is such that: dM/dt = n (approx)

// interface for fsolveb 
function [dsma] = tle_F(n, ind, args)
  dsma = CLx_tle_unkozai(n, args.ecc(ind), args.inc(ind), args.whichconst) - args.sma(ind); 
endfunction

// check extension module is loaded
CL__checkExtensionModule(); 

if (cmd == "n2a")
  // x is the mean motion 
  // => compute sma
  y = CLx_tle_unkozai(x, ecc, inc, whichconst); 

elseif (cmd == "a2n")
  // x is the semi major axis 
  // => compute n

  // arguments for fsolveb
  args = struct();  
  args.ecc = ecc;
  args.inc = inc;
  args.whichconst = whichconst; 
  args.sma = x; // objective sma

  // computes n (consistent with TLE convention)
  // n1 and n2: bounds
  // ytol: tolerance on sma (meters) 
  // => relative tolerance ~ 3.e-15
  // (tested for orbits of altitude between 0 and 40000km)
  cst = CLx_tle_getConst(whichconst);
  n1 = 0.9*sqrt(cst.mu ./ args.sma.^3);
  n2 = 1.1*sqrt(cst.mu ./ args.sma.^3);
  // Relative tolerance on Earth surface: 3.e-15
  // (tested for orbits of altitude between 0 and 40000km)
  ytol = 2.e-8; 
  
  y = CL_fsolveb(tle_F, n1, n2, args, ytol=ytol, meth="ds");

else
  CL__error("Invalid value for argument cmd");
  
end

endfunction