%% equazione del moto problema 4 corpi sole, luna, terra, sat
function dsdt = eqmoto_nbody(t,s,n_body,t_ref)                             % n_body indica il numero di corpi che considero
                                                                           % 2 = terra e sat, 3 = aggiungo luna, 4 = aggiungo anche sole
% costanti gravitazionali
mu_e = cspice_bodvrd('399','GM',1);
mu_m = cspice_bodvrd('301','GM',1);
mu_s = cspice_bodvrd('10','GM',1);

% estrazione stato
r = s(1:3);
v = s(4:6);

% accelerazioni luna e sole fissati a 0 per problema dei 2 corpi
dvdt_moon = 0;
dvdt_sun = 0;

% modifica accelerazioni se problema a più corpi
if n_body > 2
    r_moon = cspice_spkpos('moon', t+t_ref, 'J2000', 'NONE', 'earth');     % posizione luna efemeridi
    r_rel_moon = r_moon-r;                                                 % posizione relativa
    dvdt_moon = mu_m*(r_rel_moon/norm(r_rel_moon)^3-r_moon/norm(r_moon)^3);% contributo luna

    if n_body > 3
        r_sun  = cspice_spkpos('sun', t+t_ref, 'J2000', 'NONE', 'earth');  % posizione sole efemeridi
        r_rel_sun = r_sun-r;                                               % posizione relativa
        dvdt_sun = mu_s*(r_rel_sun/norm(r_rel_sun)^3-r_sun/norm(r_sun)^3); % contributo sole
    end
end


dvdt_earth = -r*mu_e/norm(r)^3;                                            % equazione del moto problema due corpi accelerazioni
dvdt = dvdt_earth+dvdt_moon+dvdt_sun;                                      % aggiungo contributi
dsdt = [v; dvdt];                                                          % vettore completo da integrare
end