%% codice per effettuare più ottimizzazioni di seguito

clc
clear
close all
load_kernels

t0 = cspice_str2et('2026-jan-1 00:00:00');                                 % tempo iniziale di riferimento
errH = 5;                                                                  % max errore di quota
errI = 0.5;                                                                % max errore di inclinazione
h = [100,2e5,2.5e5];                                                       % quota dell'orbita bassa lunare [km]
N = 4;                                                                   % numero di impulsi
pop_size = [50,200];                                                      % dimensione popolazione
gen_max = 25;                                                              % massimo numero di generazioni
plt = 2;                                                                   % flag dei plot 
dV_max = 0.1;                                                              % modulo massimo velocità
t_max = [30,20];                                                           % tempo massimo della trasferta
ref = 2;                                                                   % riferimento luna

INC = [20,30];

for i = 1:length(INC)
    inc = INC(i);
    [X_GA_i, X_sqp_i, Inc_i, Quota_i, DV_tot_i, DV_i, Dt_i] = ...
        ottimizzazione(t0,@fitness_NRHO_LLO,@vincoli_LLO,h,inc,N,errH,errI, ...
        pop_size,gen_max,dV_max,t_max,plt,ref);
    
    X_GA{i}   = X_GA_i;
    Xsqp{i}  = X_sqp_i;
    Inc{i}    = Inc_i;
    Quota{i}  = Quota_i;
    DV_tot{i} = DV_tot_i;
    DV{i}     = DV_i;
    Dt{i}     = Dt_i;
end

%%
dV_sqp = cellfun(@(x) x(2), DV_tot);
figure; bar(dV_sqp); hold on; grid on;
xticklabels(INC);
xlabel('Inclinazione [°]'); ylabel('\DeltaV [km/s]');
title('DV totale per inclinazione');