function [TW_LO,TW_climb,TW_subcrBR,TW_subcrBE] = matching_chart(Wgtot,WS,Altitude,Mach)
    
    Wgtot=Wgtot*0.4536; % kg
    WS=WS*4.8828; % kg/m^2

    %TW di Lift Off
    sigma_0=1;
    [T_0,rho_0] = atmosphere(0);
    l_LO=4000; % m 
    CL_LO=0.18;
        
    TW_LO=@(WS) WS/(rho_0*sigma_0*l_LO*CL_LO);

    %TW di Climb
    [T_climb,rho_climb] = atmosphere(Altitude/4);
    a_climb=sqrt(287*T_climb*1.4); % m/s
    Mach_climb=0.50; 
    q_climb=0.5*(Mach_climb*a_climb)^2*rho_climb; % Pa
    Gsubclimb=0.03;
    Cd0=0.02;
    sigma_climb=rho_climb/rho_0;
    manetta_climb=1;
        
    TW_climb=@(WS) (q_climb*Cd0./(WS*9.81)+Gsubclimb)/manetta_climb/sigma_climb;

    %TW di cruise
    [T_subcr,rho_subcr] = atmosphere(Altitude/2);
    a_subcr=sqrt(287*T_subcr*1.4); % m/s
    Mach_subcr=0.7;
    q_subcr=0.5*(Mach_subcr*a_subcr)^2*rho_subcr;
    sigma_subcr=rho_subcr/rho_0;
    manetta_subcr=1;

    TW_subcrBR=@(WS) (q_subcr*(4/3)*Cd0./(WS*9.81))/manetta_subcr/sigma_subcr;
    TW_subcrBE=@(WS)(q_subcr*2*Cd0./(WS*9.81))/manetta_subcr/sigma_subcr;


    x1=50;
    x2=800;
    x=linspace(x1,x2,10000);
    
    figure
    plot(x,TW_LO(x),'-','Color','black')
    hold on
    plot(x,TW_climb(x),'-','Color','blue')
    hold on
    plot(x,TW_subcrBE(x),'-','Color','green')
    hold on
    plot(x,TW_subcrBR(x),'-','Color','magenta')
    hold on
    plot([WS,WS,WS,WS,WS,WS],[0,TW_LO(WS),TW_climb(WS),TW_subcrBR(WS),TW_subcrBE(WS),2],'-o','Color','red')
    title('Matching Chart, regione subsonica')
    xlabel('W/S')
    ylabel('T/W')
    legend('TW LO','TW climb','TW subcrBE','TW subcrBR','WS')
    
    %TW di cruise
    [T_supcr,rho_supcr] = atmosphere(Altitude);
    a_supcr=sqrt(287*T_supcr*1.4); % m/s
    q_supcr=0.5*(Mach*a_supcr)^2*rho_supcr;
    sigma_supcr=rho_supcr/rho_0;
    manetta_supcr=1;
    Cd0=0.004;

    TW_supcrBR=@(WS) (q_supcr*(4/3)*Cd0./(WS*9.81))/manetta_supcr;
    TW_supcrBE=@(WS)(q_supcr*2*Cd0./(WS*9.81))/manetta_supcr;

    figure
    plot(x,TW_supcrBE(x),'-','Color','green')
    hold on
    plot(x,TW_supcrBR(x),'-','Color','magenta')
    hold on
    plot([WS,WS,WS,WS],[0,TW_supcrBR(WS),TW_supcrBE(WS),1],'-o','Color','red')
    title('Matching Chart, regione supersonica')
    xlabel('W/S')
    ylabel('T/W')
    legend('TW supcrBE','TW supcrBR','WS')
    
    TW_LO=TW_LO(WS);
    TW_climb=TW_climb(WS);
    TW_subcrBE=TW_subcrBE(WS);
    TW_subcrBR=TW_subcrBR(WS);
end