%% funzione di propagazione della traiettoria con problema N corpi a N impulsi partenza NRHO
function [Yt,t_t,Xf_t,if_t,ie] = propagazione_traj(X,t0)

% X = [t1, dt(i) dV(i)], X0 stato iniziale, t0 tempo iniziale
% --> Yt contiene lo stato istante per istante di tutte le traiettorie (N-1)
% --> t_t gli istanti di tutte le traiettorie,
% --> Xf_t continene gli stati finali delle traiettorie prima del dV 
% --> if_t contiene gli indici degli istanti finali delle traiettorie 

% estraggo numero di impulsi e divido tempi e dV
N = length(X)/4+3/4;
t_v = X(1:N)*60*60*24;
dv = X(N+1:end);
dV = reshape(dv,3,[]);                                                     % matrice che ha i dV(i) nelle colonne                                                     

% condizioni iniziali
t_t = [];                                                                  % set output tempi
Yt = [];                                                                   % set output stati
ie = [];                                                                   % nessuna collisione iniziale
tnew = t_v(1)+t0;
Xnew = cspice_spkezr('-60000',tnew,'j2000','none','earth');                % partenza da NRHO

% integrazione delle traiettorie
for i = 1:N-1                                                              % N-1 traiettorie
    % set parametri da integrare
    Xin = Xnew + [0; 0; 0; dV(:,i)];                                       % impongo la manovra
    tin = tnew;                                                            % aggiorno t iniziale
    tfin = tin + t_v(i+1);                                                 % estraggo t finale

    % INTEGRAZIONE
    [t,Y,ie_0] = int_eqmoto(Xin,tin,tfin,1,4);
    Xfin = Y(end,:)';                                                      % estraggo stato finale

    % aggiorno collisione
    if ~isempty(ie_0)
        ie = 1;
    end

    % aggiorno output
    t_t = [t_t; t];
    Yt = [Yt; Y];
    Xf_t(i,:) = Xfin';
    if_t(i) = length(t_t);

    % aggiorno condizioni di partenza nex traj.
    Xnew = Xfin;
    tnew = tfin;
end

