function [] = Ciclo_dati_modificati(i,ii,rho_a0,mf0,ISP0)
tic
addpath("funzioni\")
addpath("Input\")

%STEP 1 definizione degli input

[Wgtot,Vtot,Lb,Dbe,Fr,mf,ULF,Qmax,AR,lambda,tsuc,lambda_mezzi,Ahfp,Avfp,Wins,...
 Wa,Ttott,Nengtj,Nengsj,Nengrj,Nengtr,Htsjm,rho_f,rho_th,rho_a,eta_vol,theta_f,theta_r,...
 Alorb,ISP,rho_Wpay,rho_Vpay,delta,Mach,Altitude,passengers,Rt,sLFL,alpha_land,...
 Mach_land,x_lesuLb,croot,air_intakesuWspan] = read_input(i);


Wpay=rho_Wpay*passengers;
Vpay=Wpay/rho_Vpay;
Vaf=0;
WS=carico_alare(sLFL,AR,Mach_land,lambda_mezzi,alpha_land,i);

%impongo alcuni dati
if ii==1
rho_a=rho_a0*0.06242691;
ISP=ISP0;
mf=mf0;
end

%STEP 2 calibrazione

[kn,kc,kb] = calibrazione(Dbe,theta_f,Alorb,theta_r,Vtot,Lb);

k_vector=[kn,kc,kb];

%STEP 3
%imposto il ciclo iterativo per Vtot e Wgtot
flag=0;
q=1000;
% definisco i vettori per velocizzare il ciclo
check_breguet=zeros(q,5);
matrix=zeros(q,34);
matrix_converted=zeros(q,34);
input_iterati=zeros(q,4);

for k=1:q

    input_kmeno1=[Wgtot,Vtot,Lb,Dbe];

    %STEP 3.1
    %definizione geometria
    [Lb,Sbtot,Dbe,Fr] = geometria(kb,Fr,Vtot,eta_vol,kc,kn);
    Sref=Wgtot/WS;
    

    %STEP 3.2.1
    %Calcolo frazione di Fuel
   
    tau=Vtot/Sref^1.5;
    [WH2,E] = breguet(Rt,Altitude,ISP,Mach,tau);

    check_breguet(k,:)=[tau,Sref,Vtot,E,WH2];

    %STEP 3.2.2
    %definizione dei pesi
    %calcolo pesi strutturali
    
    [Wtps,Wstr,Wfuel,Wb,Ww,Wfinh,Wfinv,Wgear,Wthrua,Swfh,Swfv]...
    = pesi_strutturali(Lb,ULF,Dbe,Qmax,Sbtot,Wgtot,WH2,mf,Sref,AR,...
    lambda,tsuc,lambda_mezzi,Ahfp,Avfp,Ttott,Wins,rho_f,rho_th,delta);
    W_tail=Wfinh+Wfinv;
    Tail_Area=Swfh+Swfv;


    %calcolo pesi motori
    %ATTENZIONE: CONTROLLARE I RAMJET
   [Weng,Wpros,Wtnk,Wttj,Wttr,Wtrj,Wtsj] = pesi_propulsione(Nengtj,Nengtr,Nengsj,Nengrj,Wa,Htsjm,Wfuel...
    ,rho_f,rho_th,Ttott);
    
    %calcolo pesi sottosistemi
    Wspan=(AR*Sref)^0.5;
    [Wsub,Wtaves,Whydr,Welect,Wequip] = pesi_sottosistemi(Sref,Swfv,Swfh,Qmax,Lb,Wspan,Wgtot);
    
   
    %STEP 3.3
    %ricalcolo il peso totale ed il volume
    Wgtot=Wfuel+Wstr+Wpay+Wpros+Wsub;
    Vtot=(Wgtot-delta*Wfuel-Wpay-delta*Wtnk-Wtps)/rho_a+delta*(Wfuel/rho_f)+Vpay+Vaf;

    %STEP 3.4
    % rinnovo la calibrazione
    [kn,kc,kb] = ricalibrazione(Dbe,theta_f,Alorb,theta_r,Vtot,Lb,Wgtot,...
    delta,Wfuel,Wtnk,Wpay,rho_a,rho_f,Vpay,Vaf);
    
    k_vector(k+1,:)=[kn,kc,kb];

     matrix(k,:)=[Lb Dbe Fr Sref Wspan AR WS Swfh Swfv Tail_Area... 
        Sbtot Vtot Vpay Wpay Wtnk Wttj Wttr Wtrj Wtsj Wpros...
        Wb Ww W_tail Wtps Wgear Wthrua Wstr...
        Wfuel Wtaves Whydr Welect Wequip Wsub Wgtot...
        ];
    matrix_converted(k,:)=[Lb*0.3048 Dbe*0.3048 Fr Sref*0.0929 Wspan*0.3048...
        AR WS*4.8825 Swfh*0.0929 Swfv*0.0929 Tail_Area*0.0929 Sbtot*0.0929 Vtot*0.0283 Vpay*0.0283...
        Wpay*0.4536 Wtnk*0.4536 Wttj*0.4536 Wttr*0.4536 Wtrj*0.4536 Wtsj*0.4536 Wpros*0.4536...
        Wb*0.4536 Ww*0.4536 W_tail*0.4536 Wtps*0.4536 Wgear*0.4536 Wthrua*0.4536 Wstr*0.4536...
        Wfuel*0.4536 Wtaves*0.4536 Whydr*0.4536 Welect*0.4536 Wequip*0.4536 Wsub*0.4536 Wgtot*0.4536...
        ];

    input_k=[Wgtot,Vtot,Lb,Dbe];
    input_iterati(k,:)=(input_k-input_kmeno1)./input_kmeno1;
    
    %CONTROLLO DI DIVERGENZA
    if tau<=0.005 || E>100
        flag=1;
        break
    end
    
    if k>10
        if abs(matrix_converted(k,34)-matrix_converted(k-1,34))<=1 %tolleranza sul Wgtot in kg
            break
        end
    end
    
end

% riduco i vettori
check_breguet=check_breguet(1:k,:);
matrix=matrix(1:k,:);
matrix_converted=matrix_converted(1:k,:);
input_iterati=input_iterati(1:k,:);

% RISULTATI

[TW_LO,TW_climb,TW_subcrBR,TW_subcrBE] = matching_chart(Wgtot,WS,Altitude,Mach);
write_output(check_breguet,matrix,matrix_converted,k,input_iterati,TW_LO,Weng);

Sketch2D(Lb,theta_f,Dbe,Mach,x_lesuLb,lambda_mezzi,Wspan,...
    lambda,croot,theta_r,air_intakesuWspan,i)

if flag==1
    "Il ciclo iterativo diverge"
else
    "Il ciclo iterarivo converge"
end
toc
end

