Cannot plot the full interval

function [] = FuelJetPlot

H_climb = 0:500:35000;         % Altitude [ft]
for k1 = 1:length(H_climb)
% CLIMB PHASE

% Thrust calculation

C_Tc_1 = .75979E+06;
C_Tc_2 = .52423E+05;
C_Tc_3 = .40968E-10;



Thr_jet_climb_ISA (k1)= C_Tc_1* (1-(H_climb(k1)/C_Tc_2)+C_Tc_3*H_climb(k1)^2);     %  Maximum climb and take-off thrust for jet engine [N]
Thr (k1)= Thr_jet_climb_ISA(k1)/1000;       % [N] -> [kN]

% Vtas calculation

Vcl_1 = 335;             % Standard calibrated airspeed [kt]
Vcl_2 = 172.3;           % Standard calibrated airspeed [kt] -> [m/s] (To find the  Mach transition altitude)
Vcl_2_in_knots = 335;    % Standard calibrated airspeed [kt] (To find the result in knots, if altitude is between 10,000 ft and Mach transition altitude)
M_cl  = 0.86;            % Standard calibrated airspeed [kt]

K = 1.4;                 % Adiabatic index of air
R = 287.05287;           % Real gas constant for air [m2/(K·s2)]
Bt = - 0.0065;           % ISA temperature gradient with altitude below the tropopause [K/m]
T0 = 288.15;             % Standard atmospheric temperature at MSL [K]
g0 = 9.80665;            % Gravitational acceleration [m/s2]
a0= 340.294;             % Speed of Sound [m/s]

CV_min = 1.3;            % Minimum speed coefficient 
Vstall_TO = .14200E+03;  % Stall speed at take-off [KCAS]

Vd_CL1 = 5;              % Climb speed increment below 1500 ft (jet)
Vd_CL2 = 10;             % Climb speed increment below 3000 ft (jet)
Vd_CL3 = 30;             % Climb speed increment below 4000 ft (jet)
Vd_CL4 = 60;             % Climb speed increment below 5000 ft (jet)
Vd_CL5 = 80;             % Climb speed increment below 6000 ft (jet)


% CLIMB SPEEDS

% 1) Jet aircraft

CAS_climb  = Vcl_2;   
Mach_climb = M_cl;
delta_trans = (((1+((K-1)/2)*(CAS_climb/a0)^2)^(K/(K-1)))-1)/(((1+(K-1)/2*Mach_climb^2)^(K/(K-1)))-1);    % Pressure ratio at the transition altitude

teta_trans = delta_trans ^ (-Bt*R/g0);      % Temperature ratio at the transition altitude

H_p_trans_climb = (1000/0.348/6.5)*(T0*(1-teta_trans));   %  Transition altitude for climb [ft]

if (0<=H_climb(k1)&&H_climb(k1)<=1499)
Vnom_climb_jet(k1) = CV_min * Vstall_TO + Vd_CL1;

elseif (1500<=H_climb(k1)&&H_climb(k1)<=2999)
        Vnom_climb_jet(k1) = CV_min * Vstall_TO + Vd_CL2;
        
elseif (3000<=H_climb(k1)&&H_climb(k1)<=3999)
        Vnom_climb_jet(k1) = CV_min * Vstall_TO + Vd_CL3; 
        
elseif (4000<=H_climb(k1)&&H_climb(k1)<=4999)
        Vnom_climb_jet(k1) = CV_min * Vstall_TO + Vd_CL4;
        
elseif (5000<=H_climb(k1)&&H_climb(k1)<=5999)
        Vnom_climb_jet (k1)= CV_min * Vstall_TO + Vd_CL5;
        
elseif (6000<=H_climb(k1)&&H_climb(k1)<=9999)
        Vnom_climb_jet(k1) = min(Vcl_1,250);
        
elseif (10000<=H_climb(k1)&&H_climb(k1)<=H_p_trans_climb)
        Vnom_climb_jet(k1) = Vcl_2_in_knots;
        
elseif (H_p_trans_climb<H_climb(k1))
        Vnom_climb_jet (k1)= M_cl;
end

Vcas(k1) = Vnom_climb_jet(k1)* 0.514;               % [kn] -> [m/s]

H_climb (k1)= H_climb(k1) * 0.3048;         % [feet] -> [m]
K = 1.4;                                    % Adiabatic index of air
R = 287.05287;                              % Real gas constant for air [m2/(K·s2)]
Bt = - 0.0065;                              % ISA temperature gradient with altitude below the tropopause [K/m]
deltaT = 0;                                 % Value of the real temperature T in ISA conditions [K]
T0 = 288.15;                                % Standard atmospheric temperature at MSL [K]
P0 = 101325;                                % Standard atmospheric pressure at MSL [Pa]
g0 = 9.80665;                               % Gravitational acceleration [m/s2]
p0 = 1.225;                                 % Standard atmospheric density at MSL [kg/m3]
visc = (K-1)./K;

T(k1) = T0 + deltaT + Bt * H_climb(k1);         % Temperature [K]
P (k1)= P0*((T(k1)-deltaT)/T0).^((-g0)/(Bt*R));     % Pressure [Pa]
p (k1)= P(k1) ./ (R*T(k1));                             % Density [kg/m^3]
Vtas(k1) = (2*P(k1)/visc/p(k1)*((1 + P0/P(k1)*((1 + visc*p0*Vcas(k1)*Vcas(k1)/2/P0).^(1/visc)-1)).^(visc)-1)).^(1/2); % True Air Speed [m/s]
Vtas (k1)= Vtas(k1) * 1.94;                         % [m/s] -> [knots]
% Fuel consumption calculation

C_f1 = .58259E+00;
C_f2 = .12839E+04;

n_jet(k1) = C_f1*(1+Vtas(k1)/C_f2);                     %  Thrust specific fuel consumption for jet engine  [kg/(min·kN)]
f_nom_jet(k1) = n_jet(k1) * Thr(k1);                    %  Nominal fuel flow for jet engine [kg/min], can be used in climb phase
end

H_cruise = 0:500:35000;         % Altitude[ft]
for k2 = 1:length(H_cruise)
    
% CRUISE PHASE

% Speed Calculation

% CRUISE SPEEDS

Vcr_1 = 320;             % Standard calibrated airspeed [kt]
Vcr_2 = 164.62;          % Standard calibrated airspeed [kt] -> [m/s] (To find the  Mach transition altitude)
Vcr_2_in_knots = 320;    % Standard calibrated airspeed [kt] (To find the result in knots, if altitude is between 10,000 ft and Mach transition altitude)
M_cr  = 0.84;            % Standard calibrated airspeed [kt]

K = 1.4;                 % Adiabatic index of air
R = 287.05287;           % Real gas constant for air [m2/(K·s2)]
Bt = - 0.0065;           % ISA temperature gradient with altitude below the tropopause [K/m]
%deltaT = 0;             % Value of the real temperature T in ISA conditions [K]
T0 = 288.15;             % Standard atmospheric temperature at MSL [K]
g0 = 9.80665;            % Gravitational acceleration [m/s2]
a0= 340.294;             % Speed of Sound [m/s]

% 1) Jet Aircraft

CAS_cruise  = Vcr_2;   
Mach_cruise = M_cr;
delta_trans = (((1+((K-1)/2)*(CAS_cruise/a0)^2)^(K/(K-1)))-1)/(((1+(K-1)/2*Mach_cruise^2)^(K/(K-1)))-1);    % Pressure ratio at the transition altitude

teta_trans = delta_trans ^ (-Bt*R/g0);      % Temperature ratio at the transition altitude

H_p_trans_cruise = (1000/0.348/6.5)*(T0*(1-teta_trans));   %  Transition altitude for cruise [ft]

if (0<=H_cruise(k2)&&H_cruise(k2)<=2999)
Vnom_cruise_jet(k2) = min(Vcr_1,170);

elseif (3000<=H_cruise(k2)&&H_cruise(k2)<=5999)
        Vnom_cruise_jet(k2) = min(Vcr_1,220);
        
elseif (6000<=H_cruise(k2)&&H_cruise(k2)<=13999)
        Vnom_cruise_jet(k2) = min(Vcr_1,250);
        
elseif (14000<=H_cruise(k2)&&H_cruise(k2)<=H_p_trans_cruise)
        Vnom_cruise_jet(k2) = Vcr_2_in_knots;
 
elseif (H_p_trans_cruise < H_cruise(k2))
        Vnom_cruise_jet(k2) = M_cr;
end

Vcas (k2)= Vnom_cruise_jet(k2) * 0.514;                     %[knots] -> [m/s]

% Thrust Calculation

L = .44225E+06;                             % .44225E+03 tons = .44225E+06 kg and W = L.
S = .51097E+03;                             % Surface Area [m^2]
H_cruise(k2) = H_cruise(k2) * 0.3048;       % [feet] -> [m]
K = 1.4;                                    % Adiabatic index of air
R = 287.05287;                              % Real gas constant for air [m2/(K·s2)]
Bt = - 0.0065;                              % ISA temperature gradient with altitude below the tropopause [K/m]
deltaT = 0;                                 % Value of the real temperature T in ISA conditions [K]
T0 = 288.15;                                % Standard atmospheric temperature at MSL [K]
P0 = 101325;                                % Standard atmospheric pressure at MSL [Pa]
g0 = 9.80665;                               % Gravitational acceleration [m/s2]
p0 = 1.225;                                 % Standard atmospheric density at MSL [kg/m3]
visc = (K-1)./K;

T (k2)= T0 + deltaT + Bt * H_cruise(k2);        % Temperature [K]
P (k2)= P0*((T(k2)-deltaT)/T0).^((-g0)/(Bt*R));     % Pressure [Pa]
p (k2)= P(k2) ./ (R*T(k2));                             % Density [kg/m^3]

Vtas(k2) = (2*P(k2)/visc/p(k2)*((1 + P0/P(k2)*((1 + visc*p0*Vcas(k2)*Vcas(k2)/2/P0).^(1/visc)-1)).^(visc)-1)).^(1/2); % True Air Speed [m/s]

Cl(k2) = 2*L / (p(k2)*Vtas(k2)*Vtas(k2)*S);                 % Lift coefficient
C_D0cr = .25669E-01;
C_D2cr = .39082E-01;
C_D(k2) = C_D0cr + C_D2cr * (Cl(k2)).^2;            % Drag Coefficient
D (k2)= 1/2 * p(k2) * Vtas(k2)* Vtas(k2) * S * C_D(k2);         % Drag Force [N]

Thr(k2) = D(k2)/1000;                               % Cruise drag (Thrust) [N] -> [kN]
Vtas(k2) = Vtas(k2) * 1.94;                         % [m/s] -> [knots]

% Fuel Consumption Calculation

C_fcr = .89625E+00;
C_f1 = .58259E+00;
C_f2 = .12839E+04;

n_jet(k2) = C_f1*(1+Vtas(k2)/C_f2);                 %  Thrust specific fuel consumption for jet engine  [kg/(min·kN)]
f_cr_jet(k2) = n_jet(k2) * Thr(k2) * C_fcr;         %  Cruise fuel flow [kg/min]
end

H_descent = 0:500:35000;         % Altitude [ft]
for k3 = 1:length(H_descent)

% DESCENT PHASE

% The minimum fuel flow, fmin [kg/min], corresponds to idle thrust descent conditions for both jet and turboprop engines
% so Vtas and Thrust calculations are unnecessary.

% Fuel Consumption Calculation

C_f3 = .45380E+02;
C_f4 = .72929E+05;

f_min_jet(k3) = C_f3 * (1-H_descent(k3)/C_f4);   % Minimum fuel flow [kg/min]
end

figure(1)
plot(H_climb,f_nom_jet);
xlabel('Altitude [ft]');
ylabel('Nominal fuel flow for jet engine [kg/min]');
title('H vs f for climb phase');

figure(2)
plot(H_cruise,f_cr_jet);
xlabel('Altitude [ft]');
ylabel('Cruise fuel flow for jet engine [kg/min]');
title('H vs f for cruise phase');

figure(3)
plot(H_descent,f_min_jet);
xlabel('Altitude [ft]');
ylabel('Minimum fuel flow for jet engine [kg/min]');
title('H vs f for descent phase');

end

Hi, the code plots fuel consumption for climb, cruise and descent phases of a jet engine aircraft with respect to altitude. Whenever I run the ccode, Fig1 and Fig2's altitude stops at H=x=10668 ft but I want them to end at H=x=35,000 ft. How to fix this problem?

And if you see any mistakes in cruise phase part of the code, please inform me because its plot looks wrong.

Thanks.

Answer 1

[Simon Chan]

You converted H_climb & H_cruise from [feet] to [m] by multiplying the value by 0.3048.

However, there is no such conversion for variable H_descent.

You can just remove the conversion and get your required figures.

%H_climb (k1)= H_climb(k1) * 0.3048;         % [feet] -> [m]
%
%H_cruise(k2) = H_cruise(k2) * 0.3048;       % [feet] -> [m]

[Turgut Ataseven]

Thanks for pointing the mistake out!

Those conversions were necessary for climb and cruise phases, so I redefined altitudes at the last row of climb and cruise phases as

H_climb = 0:500:35000;          % Redefining altitude[ft] because it was converted to [m] above.
H_cruise = 0:500:35000;         % Redefining altitude[ft] because it was converted to [m] above.

Hope it won't change the calculations.