//  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 computes the derivatives of a satellite state at specified times using STELA
// Parameters:
// - cjd: vector of dates at which the derivatives are computed (Calendar Julian Date), TREF time scale (1xN double)
// - x: vector of satellite state for dates specified in cjd (1xN double)
// - params: parameters using for computing (struct)
// - res: ignored for moment
// - ut1_tref: UT1 minus TREF (double)
// - tt_tref: TT minus TREF (double)
// Returns:
// - dxdt: derivatives of satellite state (6xN 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 [dxdt, info] = CLx_stela_deriv(cjd, x, params, res, ut1_tref, tt_tref)
//=============================================================================
    
    // Default outputs
    dxdt = [];
    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 size: cjd=1xN and x=6xN
    if (size(cjd,1) <> 1 | size(x,1) <> 6 | size(cjd,2) <> size(x,2))
        CLx__error("Invalid arguments sizes (dates or state vectors).");
    end
    
    // initializes vector of derivatives (same size as x)
    // Only columns of x with no %nan are considered    
    dxdt = %nan * ones(x); 
    
    Inotnan = find(~isnan(sum(x, "r"))); 
    cjd = cjd(Inotnan); 
    x = x(:,Inotnan); 
    
    if (cjd == [])
      return; 
    end

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

    // Get a Derivation instance
    jimport fr.cnes.celestlab.celestlabx.stela.Derivation;
    jderivation = jnewInstance(Derivation);

    // Compute derivatives
    jinformations = jinvoke(jderivation, "derivate", cjd, x, 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 computation succeeded, get the results
    if (info.status >= 0) 
        res = jinvoke(jderivation, "getDerivativeMatrix");
        dxdt(:,Inotnan) = res; 
    end

    // Clear memory - should not be necessary
    jremove(Derivation, jderivation, jinformations); 
    
endfunction