%% funzione per calcolare il dV per circolarizzare un'orbita
function dV = circolarizzazione(Xfin, parametro, t,inc)
% funzione per calcolare il dV per circolarizzare un'orbita
% parametro = 1 se terra, 2 se luna

% condizione di default
if nargin < 2
    parametro = 1;
    t = 0;
    inc = [];
end
if nargin < 4
    inc = [];
end

% calcolo parametri orbitali
kep1 = state2kep(Xfin,parametro,t);

% impongo dV di circolarizzazione
if parametro == 2
    s_moon = cspice_spkpos('moon',t,'j2000','none','earth');               % posizione luna
    r_rel = Xfin(1:3)-s_moon;                                            % posizione relativa
    r = norm(r_rel);
else
    r_rel = Xfin(1:3);
    r = norm(r_rel);
end
kep_targ = kep1;
kep_targ(1:2) = [r;0];                                                     % impongo eccentricità nulla
state_fin = kep2state(kep_targ,parametro,t);                               % calcolo stato finale dell'orbita circolare

% dV di circolarizzazione
dV = state_fin(4:6) - Xfin(4:6);      

if ~isempty(inc)
    kep_fin = state2kep(state_fin);
    inc1 = kep_fin(3)*180/pi;
    sin_lat = state_fin(3)/norm(state_fin(1:3));
    lat = asin(sin_lat)*180/pi;
    inc2 = inc1;
    dteta = 0;
    i = 0;
    if ~(inc >= abs(lat) && inc <= 180 - abs(lat))
        % manovra impossibile
        aa = [abs(lat),180-abs(lat)];
        [~,ix] = min(abs(aa-inc));
        inc = aa(ix);
    end
    V2 = state_fin(4:6);
    while abs(inc-inc2) > 0.1
        dtet = abs(inc-inc2)/2;
        if inc > inc2
            dteta = dteta-dtet;
        else
            dteta = dteta+dtet;
        end
        V2 = ruota_vettore(state_fin(4:6),state_fin(1:3),dteta);
        state2 = [state_fin(1:3);V2];
        kep2 = state2kep(state2);
        inc2 = kep2(3)*180/pi;
        i = i +1;
        if i >1e5
            break
        end
    end
    dV = V2-Xfin(4:6);


end