Is this how I find the area under my graph?

I have this code which takes too long to run so I would save you the trouble and just attach a snippet of how I plotted the graph and called the trapz for finding the area. However, I did this but my answer seems quite far off from my caluclations so I am wondering if this was the right method to find the area under the graph. I'd like to know if I did something wrong or this answer is correct.

So this is the graph I generated and I'll attach the code below.

%% Maximum ROC vs V by taking the maximum ROC at each altitude 

% ------ Functions to solve for FSOLVE (steady climbing flight) ------ %
% func1 = v' = 1/m *(T-D-mgsin(gamma)) = 0 
% func2 = gamma' = 1/mu * (L-mgcos(gamma)) = 0
% func3 = Cm = 0 

% ------ FSOLVE variables ------ %
% x(1) = velocity
% x(2) = deltae
% x(3) = gamma

% ------ Altitude range from flight envelope ------ %
hlist_2_5 = linspace(0,9000,300);
hlist_3 = linspace(0,17000,300);

% ------ Array for storing ------ %
% 1: V   2:deltae   3:gamma   4:aoa   5:roc   6. altitude   7. thrust   8. 1/ROC 
final_list_2_5 = zeros(length(hlist_2_5),8);
final_list_3 = zeros(length(hlist_3),8);


%% PS = 2.5
%% --REMOVED-- the calculations
for j = 1:length(hlist_2_5)
    [max_roc, index] = max(store_v_de_gam_aoa_roc(:,5));
    final_list_2_5(j,5) = max_roc;  % ROC
    final_list_2_5(j,1) = store_v_de_gam_aoa_roc(index,1); % V
    final_list_2_5(j,2) = store_v_de_gam_aoa_roc(index,2); % deltae
    final_list_2_5(j,3) = store_v_de_gam_aoa_roc(index,3); % gamma
    final_list_2_5(j,4) = store_v_de_gam_aoa_roc(index,4) % aoa
    final_list_2_5(j,6) = hlist_2_5(j) % h
    final_list_2_5(j,7) = store_v_de_gam_aoa_roc(index,6) % thrust
    final_list_2_5(j,8) = 1/max_roc; % 1/ROC
end 

%% PS = 3
%% --REMOVED-- the calculations
for j = 1:length(hlist_3)
    [max_roc, index] = max(store_v_de_gam_aoa_roc(:,5));
    final_list_3(j,5) = max_roc;  % ROC
    final_list_3(j,1) = store_v_de_gam_aoa_roc(index,1); % V
    final_list_3(j,2) = store_v_de_gam_aoa_roc(index,2); % deltae
    final_list_3(j,3) = store_v_de_gam_aoa_roc(index,3); % gamma
    final_list_3(j,4) = store_v_de_gam_aoa_roc(index,4) % aoa
    final_list_3(j,6) = hlist_3(j) % h
    final_list_3(j,7) = store_v_de_gam_aoa_roc(index,6) % thrust
    final_list_3(j,8) = 1/max_roc; % 1/ROC
end

% ------ Finding time to climb to each altitude ------ %
time_2_5_climb = trapz(final_list_2_5(:,8),final_list_2_5(:,6))
% 92475.4221575 s = 25.7 hours <-- this is the result I got 
time_3_climb = trapz(final_list_3(:,8),final_list_3(:,6))
% 164693.805554 s = 45.7 hours <-- this is the result I got

% ------ Plotting graph (X vs Y)------ %
figure     % 1/ROC vs Altitude
plot(final_list_2_5(:,6),final_list_2_5(:,8))
hold on
plot(final_list_3(:,6),final_list_3(:,8))
xlabel('Altitude [m]')
ylabel('1/Rate of climb [s/m]')
title('1/Rate of climb vs Altitude')
legend('PS = 2.5','PS = 3')
grid on

[Rachel Ong]

Incase the snippet of script didnt make sense, here is the full code. However, it takes almost 600s to run.

%% Maximum ROC vs V by taking the maximum ROC at each altitude 

clc
clear all

% ------ Constants ------ %
m = 1248.5;
g = 9.81;
W = m*g;
S = 17.1;

% ------ Functions to solve for FSOLVE (steady climbing flight) ------ %
% func1 = v' = 1/m *(T-D-mgsin(gamma)) = 0 
% func2 = gamma' = 1/mu * (L-mgcos(gamma)) = 0
% func3 = Cm = 0 

% ------ FSOLVE variables ------ %
% x(1) = velocity
% x(2) = deltae
% x(3) = gamma

% ------ Altitude range from flight envelope ------ %
hlist_2_5 = linspace(0,9000,300);
hlist_3 = linspace(0,17000,300);

% ------ Array for storing ------ %
% 1: V   2:deltae   3:gamma   4:aoa   5:roc   6. altitude   7. thrust   8. 1/ROC 
final_list_2_5 = zeros(length(hlist_2_5),8);
final_list_3 = zeros(length(hlist_3),8);

% ------ Initial guess ------ %
x = [40 0 0];

%% PS = 2.5
for j = 1:length(hlist_2_5)
    PS = 2.5 
    h = hlist_2_5(j)
    [rho, speedsound] = atmos(h)
    [aoa_start,aoa_end] = aoa_guess(PS)
    aoalist = linspace(aoa_start,aoa_end,150);
    store_v_de_gam_aoa_roc = zeros(length(aoalist),6); 
    % 1: V   2:deltae   3:gamma   4:aoa   5:roc   6. thrust
    
    for i = 1:length(aoalist)
        qbar = @(x) 0.5*rho*x(1)^2*S;
        CL = @(x) 6.44*aoalist(i) + 0.355*x(2);
        L = @(x) CL(x)*qbar(x);
        CD = @(x) 0.03 + 0.05*(CL(x))^2;
        D = @(x) CD(x)*qbar(x);
        Cm = @(x) 0.05 - 0.683*aoalist(i) - 0.923*x(2);
        T = @(x) PS *( (7+x(1)/speedsound)*200/3 + h/1000*(2*x(1)/speedsound -11));

        f1 = @(x) (1/m)*(T(x)-D(x)-m*g*sin(x(3)));
        f2 = @(x) (1/(m*x(1)))*(L(x)-m*g*cos(x(3)));
        f3 = @(x) Cm(x);

        func = @(x) [f1(x);f2(x);f3(x)];

        x = fsolve(@(x) func(x), x);
        
        store_v_de_gam_aoa_roc(i,1) = x(1); % V
        store_v_de_gam_aoa_roc(i,2) = x(2); % delta e
        store_v_de_gam_aoa_roc(i,3) = x(3); % gamma
        store_v_de_gam_aoa_roc(i,4) = aoalist(i); % aoa
        
        thrust = PS * ( 200/3*(7 + x(1)/speedsound) + h/1000*(2*x(1)/speedsound - 11)  );
        liftcoeff = 6.44*aoalist(i) + 0.355*x(2);
        drag = 0.5*rho*x(1)^2*S*(0.03 + 0.05*liftcoeff^2);
        ROC = (thrust-drag)*x(1)/W;
        store_v_de_gam_aoa_roc(i,5)= ROC; % ROC
        store_v_de_gam_aoa_roc(i,6)= thrust; % ROC
    end    
     
    [max_roc, index] = max(store_v_de_gam_aoa_roc(:,5));
    final_list_2_5(j,5) = max_roc;  % ROC
    final_list_2_5(j,1) = store_v_de_gam_aoa_roc(index,1); % V
    final_list_2_5(j,2) = store_v_de_gam_aoa_roc(index,2); % deltae
    final_list_2_5(j,3) = store_v_de_gam_aoa_roc(index,3); % gamma
    final_list_2_5(j,4) = store_v_de_gam_aoa_roc(index,4) % aoa
    final_list_2_5(j,6) = hlist_2_5(j) % h
    final_list_2_5(j,7) = store_v_de_gam_aoa_roc(index,6) % thrust
    final_list_2_5(j,8) = 1/max_roc; % 1/ROC
    
end

% ------ Initial guess ------ %
x = [40 0 0];

%% PS = 3
for j = 1:length(hlist_3)
    PS = 3; 
    h = hlist_3(j);
    [rho, speedsound] = atmos(h);
    [aoa_start,aoa_end] = aoa_guess(PS);
    aoalist = linspace(aoa_start,aoa_end,150);
    store_v_de_gam_aoa_roc = zeros(length(aoalist),5); 
    % 1: V   2:deltae   3:gamma   4:aoa   5:roc   6. thrust
    
    for i = 1:length(aoalist)
        qbar = @(x) 0.5*rho*x(1)^2*S;
        CL = @(x) 6.44*aoalist(i) + 0.355*x(2);
        L = @(x) CL(x)*qbar(x);
        CD = @(x) 0.03 + 0.05*(CL(x))^2;
        D = @(x) CD(x)*qbar(x);
        Cm = @(x) 0.05 - 0.683*aoalist(i) - 0.923*x(2);
        T = @(x) PS *( (7+x(1)/speedsound)*200/3 + h/1000*(2*x(1)/speedsound -11));

        f1 = @(x) (1/m)*(T(x)-D(x)-m*g*sin(x(3)));
        f2 = @(x) (1/(m*x(1)))*(L(x)-m*g*cos(x(3)));
        f3 = @(x) Cm(x);

        func = @(x) [f1(x);f2(x);f3(x)];

        x = fsolve(@(x) func(x), x);
        
        store_v_de_gam_aoa_roc(i,1) = x(1); % V
        store_v_de_gam_aoa_roc(i,2) = x(2); % delta e
        store_v_de_gam_aoa_roc(i,3) = x(3); % gamma
        store_v_de_gam_aoa_roc(i,4) = aoalist(i); % aoa
        
        thrust = PS * ( 200/3*(7 + x(1)/speedsound) + h/1000*(2*x(1)/speedsound - 11)  );
        liftcoeff = 6.44*aoalist(i) + 0.355*x(2);
        drag = 0.5*rho*x(1)^2*S*(0.03 + 0.05*liftcoeff^2);
        ROC = (thrust-drag)*x(1)/W;
        store_v_de_gam_aoa_roc(i,5)= ROC; % ROC
        store_v_de_gam_aoa_roc(i,6)= thrust; % thrust
    end     
    
    [max_roc, index] = max(store_v_de_gam_aoa_roc(:,5));
    final_list_3(j,5) = max_roc;  % ROC
    final_list_3(j,1) = store_v_de_gam_aoa_roc(index,1); % V
    final_list_3(j,2) = store_v_de_gam_aoa_roc(index,2); % deltae
    final_list_3(j,3) = store_v_de_gam_aoa_roc(index,3); % gamma
    final_list_3(j,4) = store_v_de_gam_aoa_roc(index,4) % aoa
    final_list_3(j,6) = hlist_3(j) % h
    final_list_3(j,7) = store_v_de_gam_aoa_roc(index,6) % thrust
    final_list_3(j,8) = 1/max_roc; % 1/ROC
end

% ------ Finding time to climb to each altitude ------ %
time_2_5_climb = trapz(final_list_2_5(:,8),final_list_2_5(:,6))
% 92475.4221575 s = 25.7 hours
time_3_climb = trapz(final_list_3(:,8),final_list_3(:,6))
% 164693.805554 s = 45.7 hours

% ------ Plotting graph (X vs Y)------ %

figure     % 1/ROC vs Altitude
plot(final_list_2_5(:,6),final_list_2_5(:,8))
hold on
plot(final_list_3(:,6),final_list_3(:,8))
xlabel('Altitude [m]')
ylabel('1/Rate of climb [s/m]')
title('1/Rate of climb vs Altitude')
legend('PS = 2.5','PS = 3')
grid on

figure     % Velocity vs Altitude
plot(final_list_2_5(:,6),final_list_2_5(:,1))
hold on
plot(final_list_3(:,6),final_list_3(:,1))
xlabel('Altitude [m]')
ylabel('Velocity [m/s]')
title('Velocity vs Altitude')
legend('PS = 2.5','PS = 3')
grid on


%% ------ Atmos function ------ %%
% Calculating T, p ,rho and speedsound for every altitude in the ISA atmosphere
function [rho, speedsound] = atmos(h)
h1 = 11000;    % Height of tropopause
h2 = 20000;    % End height of table
g = 9.81;
R = 287; 
c = 6.51e-3;       % temperature lapse dt/dh = - c = -6.51 degcelcius/km

T0 = 15+273.15;       % Temperature sea level
p0 = 101325;       % pressure sealevel
rho0 = 101325/R/T0;      % density sealevel = pressure / R*T, R=287, T = 15 degcelcius       

T1 = T0 - c*h1;    % Temperature at 11km
p1 = p0 * (T1/T0)^5.2506;   % Pressure at 11km
rho1 = rho0 * (T1/T0)^4.2506;   % Density at 11km
T2 = T1;    % Temperature at 20km
p2 = p1 * exp(-g/(R*T2)*(h2-h1));  % Pressure at 20km
rho2 = rho1 * exp(-g/(R*T2)*(h2-h1));  % Density at 20km


if h <= h1
%       disp('Troposphere');
    T = T0 - c*h;
    p = p0 * (T/T0)^5.2506;
    rho = rho0 * (T/T0)^4.2506;
    speedsound = (1.4*R*T)^0.5;
elseif h <= h2
%       disp('Tropopause');
    T = T1;
    p = p1 * exp(-g/(R*T)*(h-h1));
    rho = rho1 * exp(-g/(R*T)*(h-h1));
    speedsound = (1.4*R*T)^0.5;
end

return 

end

%% ------ AOA guess function ------ %%
% To guess the range of AOA at each altitude
function [aoa_start,aoa_end] = aoa_guess(PS)
h = 0;
W = 1248.5*9.81;        % Weight of UAV
S = 17.1;       % Wing area        

%     CL = (6.44*aoa + 0.355*eledefl);
%     CD = (0.03 + 0.05*(6.44*aoa + 0.355*eledefl)^2);
%     Cm = (0.05 - 0.683*aoa - 0.923*eledefl);

%     thrust = (PS * ( (7 + V/speedsound )*200/3  + h*(2*(V/speedsound) - 11) ));

%%
% V = x(1);
% aoa = x(2);
% eledefl = x(3);
[rho, speedsound] = atmos(h);

qbar = @(x) 0.5*rho*x(1)^2*S;
CL = @(x) (6.44*x(2) + 0.355*x(3));
CD = @(x) (0.03 + 0.05*(CL(x))^2);
Cm = @(x) 0.05 - 0.683*x(2) - 0.923*x(3);
thrust = @(x) ( PS * ( 200/3*( 7 + x(1)/speedsound ) + h/1000*( 2*x(1)/speedsound - 11) ) );

func_1 = @(x) qbar(x)*CL(x) + thrust(x)*sin(x(2)) - W;
func_2 = @(x) qbar(x)*CD(x) - thrust(x)*cos(x(2));
func_3 = @(x) Cm(x);

FUNC = @(x) [func_1(x); func_2(x); func_3(x)];

% ------ Guess --> Range of AOA ------ %
Xlow = fsolve(@(x) FUNC(x), [10 0 -30/180*pi]); % high AOA
aoa_end = Xlow(1,2);

Xhigh = fsolve(@(x) FUNC(x), [100 25/180*pi -30/180*pi]); % low AOA
aoa_start = Xhigh(1,2);
end

Answer 1

[Torsten]

% ------ Finding time to climb to each altitude ------ %
time_2_5_climb = trapz(final_list_2_5(:,6),final_list_2_5(:,8))
% 92475.4221575 s = 25.7 hours <-- this is the result I got 
time_3_climb = trapz(final_list_3(:,6),final_list_3(:,8))
% 164693.805554 s = 45.7 hours <-- this is the result I got

instead of

% ------ Finding time to climb to each altitude ------ %
time_2_5_climb = trapz(final_list_2_5(:,8),final_list_2_5(:,6))
% 92475.4221575 s = 25.7 hours <-- this is the result I got 
time_3_climb = trapz(final_list_3(:,8),final_list_3(:,6))
% 164693.805554 s = 45.7 hours <-- this is the result I got

And one further hint:

fsolve will perform much faster in your loops if you choose the last solution as initial guess for the next call:

x = fsolve(@(x) func(x), x);

instead of

x = fsolve(@(x) func(x), [40 0 0]);

Dont forget to initialize

x = [40 0 0] 

in front of the two loops.

[Rachel Ong]

Hello Torsten,

Thank you. I've tried that and the answer seems to be a bit more correct. Thanks for your help. Did my previous answer find the area between the curve and the y axis instead of the x axis?

Also, I edited what you suggested. Could I check if it is in the right place?

[Torsten]

Did my previous answer find the area between the curve and the y axis instead of the x axis?

As long as final_list_2_5(:,8) and final_list_3(:,8) are strictly increasing, this should be the case, yes.

Also, I edited what you suggested. Could I check if it is in the right place?

It should be placed where the input functions for fsolve change much and continuation with the last obtained solution makes no sense because the next solution will be far off from the last. You must decide, but my guess is that x = [40 0 0] should be placed at the beginning of the inner for-loop:

for j = 1:length(hlist_3)
    PS = 3; 
    h = hlist_3(j);
    [rho, speedsound] = atmos(h);
    [aoa_start,aoa_end] = aoa_guess(PS);
    aoalist = linspace(aoa_start,aoa_end,150);
    store_v_de_gam_aoa_roc = zeros(length(aoalist),5); 
    % 1: V   2:deltae   3:gamma   4:aoa   5:roc   6. thrust
    x = [40 0 0];
    for i = 1:length(aoalist)
        qbar = @(x) 0.5*rho*x(1)^2*S;
        CL = @(x) 6.44*aoalist(i) + 0.355*x(2);
        L = @(x) CL(x)*qbar(x);
        CD = @(x) 0.03 + 0.05*(CL(x))^2;
        D = @(x) CD(x)*qbar(x);
        Cm = @(x) 0.05 - 0.683*aoalist(i) - 0.923*x(2);
        T = @(x) PS *( (7+x(1)/speedsound)*200/3 + h/1000*(2*x(1)/speedsound -11));
        f1 = @(x) (1/m)*(T(x)-D(x)-m*g*sin(x(3)));
        f2 = @(x) (1/(m*x(1)))*(L(x)-m*g*cos(x(3)));
        f3 = @(x) Cm(x);
        func = @(x) [f1(x);f2(x);f3(x)];
        x = fsolve(@(x) func(x), x);
        
        store_v_de_gam_aoa_roc(i,1) = x(1); % V
        store_v_de_gam_aoa_roc(i,2) = x(2); % delta e
        store_v_de_gam_aoa_roc(i,3) = x(3); % gamma
        store_v_de_gam_aoa_roc(i,4) = aoalist(i); % aoa
        
        thrust = PS * ( 200/3*(7 + x(1)/speedsound) + h/1000*(2*x(1)/speedsound - 11)  );
        liftcoeff = 6.44*aoalist(i) + 0.355*x(2);
        drag = 0.5*rho*x(1)^2*S*(0.03 + 0.05*liftcoeff^2);
        ROC = (thrust-drag)*x(1)/W;
        store_v_de_gam_aoa_roc(i,5)= ROC; % ROC
        store_v_de_gam_aoa_roc(i,6)= thrust; % thrust
    end     
    
    [max_roc, index] = max(store_v_de_gam_aoa_roc(:,5));
    final_list_3(j,5) = max_roc;  % ROC
    final_list_3(j,1) = store_v_de_gam_aoa_roc(index,1); % V
    final_list_3(j,2) = store_v_de_gam_aoa_roc(index,2); % deltae
    final_list_3(j,3) = store_v_de_gam_aoa_roc(index,3); % gamma
    final_list_3(j,4) = store_v_de_gam_aoa_roc(index,4); % aoa
    final_list_3(j,6) = hlist_3(j); % h
    final_list_3(j,7) = store_v_de_gam_aoa_roc(index,6); % thrust
    final_list_3(j,8) = 1/max_roc; % 1/ROC
end

(same for the other loop).

Anyway: It should only speed up the execution, but shouldn't change the result.

I hope the x(3) value in your calculations is not that important because it can only be determined up to multiples of 2*pi, I guess.

[Rachel Ong]

I see, thank you for your explanations! Really appreciate it :)

[Rachel Ong]

Hi Torsten, is there anyway for me to find the point on the x-axis values where the area under the curve = 1000 for example?

[Torsten]

f = @(t)interp1(final_list_2_5(:,6),final_list_2_5(:,8),t);
fun = @(x)integral(f,final_list_2_5(1,6),x) - 1000;
x_1000 = fzero(fun,(final_list_2_5(1,6)+final_list_2_5(end,6))/2)

[Rachel Ong]

Using this value of x_1000 can I still use integral function to find the area under the graph of another set of values instead of using final_list_2_5(:,8)) I use (final_list_2_5(:,1)) for y axis up tilll this x point?

[Torsten]

You mean

f = @(t)interp1(final_list_2_5(:,6),final_list_2_5(:,1),t);
val = integral(f,final_list_2_5(1,6),x_1000) 

?

To see how the area under the curves develop with altitude, you may also use "cumtrapz":

S = cumtrapz(final_list_2_5(:,6),final_list_2_5(:,8));
plot(final_list_2_5(:,6),S)

[Rachel Ong]

I see, I'll try that and check out cumtrapz as well. Thank you!