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

// -------------------------------
// Computation of "mean" elements (generic, internal function)
// => Determines "moy" such that: fct_extrap(moy) = par
//
// Interface: fct_extrap(moy, I, args), 
// I: indices being analysed, args: any (constant) variable
//
// typ: "kep", "cir" or "equin" only
// moy and par: type = typ
//
// Algorithm: 
// Iterative process: 
// Loop: 
//   p = fct(moy);
//   Apply corrections to moy: moy = moy - (p - par)
//   after conversion to equinoxial elements. 
// Stop when (p - par) small enough (relative difference)
// Returns %nan if no convergence. 
//
// Author: AL
// -------------------------------

function [moy] = CL__ex_inverse(typ, par, fct_extrap, args)

  // Conversion: equin -> kep or cir or equin (if cir: equin->kep->cir)
  // NB: mu not used (not included in function call)
  function [oe] = equin2oe(equin, typ)
     if (typ == "equin") 
       oe = equin;
     elseif (typ == "kep")
       oe = CL__oe_equin2kep(equin); 
     else // if (typ == "cir")
       oe = CL__oe_equin2cir(equin); 
     end
  endfunction

  // kep or cir or equin -> equin (if cir: cir->kep->equin)
  // NB: mu not used (not included in function call)
  function [equin] = oe2equin(oe, typ)
    if (typ == "equin") 
      equin = oe; 
    elseif (typ == "kep")
      equin = CL__oe_kep2equin(oe); 
    else // if (typ == "cir")
      equin = CL__oe_cir2equin(oe); 
    end
  endfunction
  
  // Check type to avoid wrong computations
  if (typ <> "cir" & typ <> "kep" & typ <> "equin")
    CL_error("Invalid argument: ''typ''"); 
  end
  
  moy = par; 
  par_equin = oe2equin(par, typ); 
  moy_equin = par_equin; 

  // (relative) tolerance ("equin" type)  
  tol = 10 * %eps * (1 + abs(par_equin)); 

  // K = indices to examine (initially: all)
  K = 1:size(par,2); 

  // relative error (ratio to tol) for convergence checking
  errel = %inf * ones(K); 

  // max number of iterations and current number
  maxiter = 50; 
  iter = 1; 

  while (iter <= maxiter & K <> [])
    // note: "K" in names means "vector has K columns" 

    oeK = fct_extrap(moy(:,K), K, args); 
    oeK_equin = oe2equin(oeK, typ);  

    deltaK = oeK_equin - par_equin(:,K); 
    deltaK(6,:) = CL_rMod(deltaK(6,:), -%pi, %pi); 

    errel(K) = max(CL_dMult(abs(deltaK), 1.0 ./ tol(:,K)), "r"); 

    moy_equin(:,K) = moy_equin(:,K) - deltaK; 
    moy(:,K) = equin2oe(moy_equin(:,K), typ); 

    K = find(errel > 1); 
    iter = iter+1; 
  end

  // if K <> []  (no convergence)  => %nan 
  moy(:,K) = %nan; 

  // adjusts "modulo" 
  // value of angles made close to values in "par"
  if (typ == "kep")
    moy(4,:) = CL_rMod(moy(4,:), par(4,:)-%pi, par(4,:)+%pi);
  end
  if (typ == "kep" | typ == "cir")
    moy(5,:) = CL_rMod(moy(5,:), par(5,:)-%pi, par(5,:)+%pi);
  end
  moy(6,:) = CL_rMod(moy(6,:), par(6,:)-%pi, par(6,:)+%pi);

endfunction