function [] = Analisi_sensibilita(i,j)
addpath("funzioni\")
addpath("Input\")

tic
%analisi di sensibilità
%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);

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

Wpay=rho_Wpay*passengers;
Vpay=Wpay/rho_Vpay;
Vaf=0;


if j==1
    Rt=9000;
end

if j==13
    delta=0.05;
end

flag=0;

% Ciclo di Analisi di sensibilità

z=20;

if j==7 || j==14
    z=10;
end

% definisco il vettore sens
if j==1 ||j==2 || j==4 || j==5 || j==7 || j==8 || j==11 || j==13 || j==14
    sens=zeros(z,4);
else
    sens=zeros(z,3);
end




for r=1:z  
    Lb=guess_inputs(1);
    Dbe=guess_inputs(2);
    Fr=Lb/Dbe;
    if r==1
    Vtot=guess_inputs(3);
    Wgtot=guess_inputs(4);
    end
    
    WS=carico_alare(sLFL,AR,Mach_land,lambda_mezzi,alpha_land,i);

    %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
    q=1000;
    % definisco i vettori per velocizzare il ciclo
    check_breguet=zeros(q,5);
    matrix=zeros(q,34);
    matrix_converted=zeros(q,34);
for k=1:q

    %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>1
        if abs(matrix_converted(k,34)-matrix_converted(k-1,34))<=5 %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,:);

    if j==1
        sens(r,:)=[Rt,Wfuel*0.4536, E, Wgtot*0.4536];
        Rt=Rt+500;
    end

    if j==2
        sens(r,:)=[passengers,Wgtot*0.4536,Vtot*0.0283, Wfuel*0.4536];
        passengers=passengers+5;
        Wpay=rho_Wpay*passengers;
    end
    
    if j==3
        sens(r,:)=[WS*4.8825,Wgtot*0.4536,Vtot*0.0283];
        WS=WS-2;
    end

    if j==4
        sens(r,:)=[rho_a/0.06242691,Wgtot*0.4536,Vtot*0.0283, Lb*0.3048];
        rho_a=rho_a-5*0.06242691;
    end
    
    if j==5
        sens(r,:)=[Mach,Wfuel*0.4536, E, Wgtot*0.4536];
        Mach=Mach+0.1;
    end

    if j==6
        sens(r,:)=[sLFL,Wgtot*0.4536,Vtot*0.0283];
        sLFL=sLFL+100;
    end

    if j==7
        sens(r,:)=[ISP,Wfuel*0.4536, E, Wgtot*0.4536];
        ISP=ISP-100;
    end

    if j==8
        sens(r,:)=[Altitude,Wfuel*0.4536, E, Wgtot*0.4536];
        Altitude=Altitude-0.5;
    end
    
    if j==9
        sens(r,:)=[guess_inputs(3)*0.0283,Wgtot*0.4536,Vtot*0.0283];
        guess_inputs(3)=guess_inputs(3)+100/0.0283;
        Vtot=guess_inputs(3);
    end

    if j==10
        sens(r,:)=[guess_inputs(4)*0.4536,Wgtot*0.4536,Vtot*0.0283];
        guess_inputs(4)=guess_inputs(4)+5000/0.4536;
        Wgtot=guess_inputs(4);
    end

    if j==11
        sens(r,:)=[AR,Wgtot*0.4536,WS*4.8825,Ww*0.4536];
        if i==5 || i==6
        AR=AR+0.04;
        else
        AR=AR-0.04;
        end
    end

    if j==12
        sens(r,:)=[Ttott*0.4536,Wgtot*0.4536,Vtot*0.0283];
        Ttott=Ttott+10^4/0.4536;
    end
        
    if j==13
        sens(r,:)=[delta,Wgtot*0.4536,Vtot*0.0283, Wfuel*0.4536];
        delta=delta+0.05;
    end

    if j==14
        sens(r,:)=[mf,Wgtot*0.4536,Ww*0.4536,Wb*0.4536];
        if i==5 || i==6
        mf=mf+0.02;
        else
        mf=mf-0.02;
        end
    end
    
    % reimposto parametri di guess 
    if j~=9
        Vtot=guess_inputs(3);
    end
    
    if j~=10
        Wgtot=guess_inputs(4);
    end

end
    
    if flag==1
        "l'analisi diverge"
    else
    
    grafici_sensibilita(sens,i,j)
    end
    toc
end



