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

// This function extrapolates an initial state vector using STELA
// Parameters:
// - cjd0: initial epoch, cjd date, TREF time scale (double)
// - x0: initial orbital elements (CelestLab type: "cireq"), Frame = MOD (6x1 double)
// - delta_t: times from cjd0 where the results are computed. Should be multiple of extrapolation step (s) (1xN double)
// - params: parameters using for computing (struct)
// - res: if equals 1, the transition matrix is computed (double)
// - ut1_tref: UT1 minus TREF (double)
// - tt_tref: TT minus TREF (double)
// Returns:
// - x: mean orbital elements (6xN double)
// - xosc: osculating orbital elements (6xN double)
// - dxd0: transition matrix = d(x)/d(x0) (6x8xN double)
// - info: informations about the computing (struct):
//    - status: 0 if the computing was completed successfully, -1 if error occurred (integer)
//    - status_msg: "OK" if the computing was completed successfully, "Error" if error occured (string)
//    - err_msg: error message (string)
//    - sim_xml: simulation parameters in STELA XML format (string)
//==================================================================================================
function [x, xosc, dxd0, info] = CLx_stela_extrap(cjd0, x0, delta_t, params, res, ut1_tref, tt_tref)
//==================================================================================================
    x = [];
    xosc = [];
    dxd0 = [];
    info = struct("status", 0, "status_msg", "OK", "sim_xml", "", "err_msg", "");

    // Check STELA is available
    if (~CLx__stelaAvailable())
      CLx__error("STELA is not available");
    end
    
    // check sizes: cjd=1xN and x=6xN
    if (size(cjd0,"*") <> 1 | size(x0,1) <> 6 | size(x0,2) <> 1)
        CLx__error("Invalid arguments sizes (initial time or state)");
    end

    // check size of delta_t
    if (size(delta_t, "*") == 0)
        return; 
        
    elseif (size(delta_t,1) <> 1)
        CLx__error("Invalid arguments sizes (final times)");
    end

    // Pass the extrapolation parameters
    jparams = CLx__createStelaParams(params);

    // Get the Extrapolation instance
    jimport fr.cnes.celestlab.celestlabx.stela.Extrapolation;
    jextrapolation = jnewInstance(Extrapolation);

    // Run extrapolation
    jinformations = jinvoke(jextrapolation, "extrapolate", cjd0, x0, delta_t, jparams, res, ut1_tref, tt_tref);

    // Get status
    status = jinvoke(jinformations, "getStatus");
    if (status == 0)
        statusMsg = "OK";
    elseif (status == 1) 
        statusMsg = "Stop";
    else
        statusMsg = "Error";
    end

    errorMessage = jinvoke(jinformations, "getErrorMessage");
    simulationXML = jinvoke(jinformations, "getSimulationXML");

    info = struct("status", double(status), ..
                  "status_msg", statusMsg, ..
                  "sim_xml", simulationXML', ..
                  "err_msg", errorMessage);

    // If extrapolation succeeded
    if (info.status >= 0) 
        // Get the params result matrixes
        jextrapolationResults = jinvoke(jextrapolation, "getExtrapolationResults");

        x = jinvoke(jextrapolationResults, "getMeanParamMatrix");
        xosc = []; // not computed any more

        // Get transition matrix result
        // Convert manually the returned jarray to hyper matrix
        if (res == 1) 
            dxd0_matrix = jinvoke(jextrapolationResults, "getTransitionMatrix");
            n = size(delta_t, '*');
            dxd0 = matrix(dxd0_matrix, [6, 8, n])
        end
        
        // Clear memory - should not be necessary !! 
        jremove(jextrapolationResults); 
    end
    
    // Clear memory - should not be necessary
    jremove(Extrapolation, jextrapolation, jinformations); 
    
endfunction