Calculate the difference between minimum values of a parabola and straight line (from a plot)

Hello everyone, please could you help me with a code.

How I can calculate the values from a plot? I need the difference between straight line (P) and between the minimum value of a parabola (P) for each curve.

So, for example in this curve, the difference between the straight blue curve and the first parabola, then the next straight line and green parabola and so on.

My code is:

clear all; close all

W = 60000;
S = 28.2;
AR=7;
cd0 = 0.02;
k = 0.04;
RC=0.51;
clalpha = 2*pi;
Psl=741000; 
hv=0:1:10;
cdminp=4*cd0;
clminp=sqrt(3*cd0/k);
Vmin=sqrt(2*W/(1.225*28.2*clminp));
D=0.5*1.225*Vmin^2*S*cdminp;
Pminreq=D*Vmin;
deltaPgiven=RC*W;
figure(1);hold on; xlabel('V');ylabel('P')
hv=0:1:10;
for k1 = 1:numel(hv)
    h = hv(k1);
    i=0;
    for alpha = 1:0.25:15
        i=i+1;
        rho(i)=1.225*exp(-h/10.4);
        cl(i) = clalpha * alpha * pi/180;
        V(i) = sqrt(2*W/rho(i)/S/cl(i));
        L(i) = 0.5 * rho(i) * V(i) * V(i) * S * cl(i);
        cd(i) = cd0 + k * cl(i) * cl(i);
        D(i) = 0.5 * rho(i) * V(i) * V(i) * S * cd(i);
        clcd(i) = cl(i)/cd(i);
        p(i) = D(i)*V(i);
        Ph(i)=Psl*(rho(i)/1.225).^0.75;
    end
    figure(1); plot(V,p)
    hold on 
    plot(V,Ph);
end

Answer 1

[Les Beckham]

Maybe this?

W = 60000;
S = 28.2;
AR=7;
cd0 = 0.02;
k = 0.04;
RC=0.51;
clalpha = 2*pi;
Psl=741000; 
hv=0:1:10;
cdminp=4*cd0;
clminp=sqrt(3*cd0/k);
Vmin=sqrt(2*W/(1.225*28.2*clminp));
D=0.5*1.225*Vmin^2*S*cdminp;
Pminreq=D*Vmin;
deltaPgiven=RC*W;
figure(1);hold on; xlabel('V');ylabel('P')
hv=0:1:10;
for k1 = 1:numel(hv)
    h = hv(k1);
    i=0;
    for alpha = 1:0.25:15
        i=i+1;
        rho(i)=1.225*exp(-h/10.4);
        cl(i) = clalpha * alpha * pi/180;
        V(i) = sqrt(2*W/rho(i)/S/cl(i));
        L(i) = 0.5 * rho(i) * V(i) * V(i) * S * cl(i);
        cd(i) = cd0 + k * cl(i) * cl(i);
        D(i) = 0.5 * rho(i) * V(i) * V(i) * S * cd(i);
        clcd(i) = cl(i)/cd(i);
        p(i) = D(i)*V(i);
        Ph(i)=Psl*(rho(i)/1.225).^0.75;
    end
    figure(1); plot(V,p)
    hold on 
    plot(V,Ph);
    [pmin, imin] = min(p); % find the min p
    deltas(k1) = Ph(1) - pmin; % calculate the difference
    tolerance = 5000; % or whatever you want
    if (abs(deltas(k1) - 300000) < tolerance)
        fprintf('delta = %8.1f at h = %4.1f, rho = %.5f, V = %.2f, Ph = %.1f, p = %.1f\n', ...
            deltas(k1), h, rho(imin), V(imin), Ph(imin), p(imin))
    end
end

delta = 302330.3 at h =  4.0, rho = 0.83387, V = 64.31, Ph = 555317.0, p = 252986.7

legend(compose('h = %.1f', hv), 'location', 'northwest')
grid on

[Alina Abdikadyr]

Yes, thank you! Also, how can I show points of P that brings me an exact delta on graph? For example, delta =3.0233. Also, can I display the values of h, rho, Vmin, Pmin for this exact delta in table?

[Les Beckham]

See edited code above

[Alina Abdikadyr]

Thank you very much! But If I will have more h values, so may I display the rows where delta is almost equal to 300000

[Les Beckham]

You are quite welcome.

I don't understand the new question. You can display whatever you want.

Maybe you mean to only display the results if delta is near 30000? How close?

You could put the fprintf statement inside an if condition. Something like

tolerance = 1000; % or whatever you want
if (abs(deltas(k1) - 30000) < tolerance)
    fprintf('delta = %8.1f at h = %4.1f, rho = %.5f, V = %.2f, Ph = %.1f, p = %.1f\n', ...
        deltas(k1), h, rho(imin), V(imin), Ph(imin), p(imin))
end

[Alina Abdikadyr]

Yes, you inderstand it correctly, where I should put this piece of code? I implmented it here:

clear all; close all
W = 60000;
S = 28.2;
AR=7;
cd0 = 0.02;
k = 0.04;
RC=0.51;
clalpha = 2*pi;
Psl=741000; 
hv=0:1:10;
cdminp=4*cd0;
clminp=sqrt(3*cd0/k);
Vmin=sqrt(2*W/(1.225*28.2*clminp));
D=0.5*1.225*Vmin^2*S*cdminp;
Pminreq=D*Vmin;
deltaPgiven=RC*W;
figure(1);hold on; xlabel('V');ylabel('P')
hv=0:1:10;
for k1 = 1:numel(hv)
    h = hv(k1);
    i=0;
    for alpha = 1:0.25:15
        i=i+1;
        rho(i)=1.225*exp(-h/10.4);
        cl(i) = clalpha * alpha * pi/180;
        V(i) = sqrt(2*W/rho(i)/S/cl(i));
        L(i) = 0.5 * rho(i) * V(i) * V(i) * S * cl(i);
        cd(i) = cd0 + k * cl(i) * cl(i);
        D(i) = 0.5 * rho(i) * V(i) * V(i) * S * cd(i);
        clcd(i) = cl(i)/cd(i);
        p(i) = D(i)*V(i);
        Ph(i)=Psl*(rho(i)/1.225).^0.75;
    end
    figure(1); plot(V,p)
    hold on 
    plot(V,Ph);
    [pmin, imin] = min(p); % find the min p
    deltas(k1) = Ph(1) - pmin; % calculate the difference
tolerance = 1000; % or whatever you want
if (abs(deltas(k1) - 30000) < tolerance)
    fprintf('delta = %8.1f at h = %4.1f, rho = %.5f, V = %.2f, Ph = %.1f, p = %.1f\n', ...
        deltas(k1), h, rho(imin), V(imin), Ph(imin), p(imin))
end
legend(compose('h = %.1f', hv), 'location', 'northwest')
grid on
end

[Les Beckham]

That's close to what I would do, but I would keep the legend and grid on outside of the for loop. I've edited my answer to reflect this change.

[Alina Abdikadyr]

Thank you! What does the tolerance value mean?

[Les Beckham]

It just represents how close to 300000 the delta has to be to make the code print the results. Adjust as desired.

Answer 2

[Torsten]

syms h V W S rho cd0 k cl Psl
eqn = V == sqrt(2*W/rho/S/cl);
cl = solve(eqn,cl);

Warning: Solutions are only valid under certain conditions. To include parameters and conditions in the solution, specify the 'ReturnConditions' value as 'true'.

cd = cd0 + k * cl^2;
D = 0.5 * rho * V * V * S * cd;
p = D*V;
Vmin = solve(diff(p,V)==0,V);
pmin = subs(p,V,Vmin);
Ph = Psl*(rho/1.225)^0.75;
hnum = 0:0.01:10;
Wnum = 60000;
Snum = 28.2;
rhonum = 1.225*exp(-hnum/10.4);
cd0num = 0.02;
knum = 0.04;
Pslnum = 741000;
for i = 1:numel(hnum)
  Phnum(i) = double(subs(Ph,[rho Psl],[rhonum(i),Pslnum]));
  pm = double(subs(pmin,[W S rho cd0 k],[Wnum Snum rhonum(i) cd0num,knum]));
  pminnum(i) = pm(end);
  deltaP(i) = Phnum(i)-pminnum(i);
end
format long
deltaP.'

ans = 1001×1

1.0e+05 * 5.322767552311991 5.316422010284389 5.310079836471793 5.303741027866013 5.297405581460750 5.291073494251587 5.284744763235998 5.278419385413340 5.272097357784856 5.265778677353669

idx = deltaP > 2.8e5 & deltaP < 3.2e5;  % select those deltaP with 2.8e5 <= deltaP < = 3.2e5
[hnum(idx).' deltaP(idx).']                 % Show these values together with the corresponding h values

ans = 77×2

1.0e+05 * 0.000036700000000 3.196913122178178 0.000036800000000 3.191616062324753 0.000036900000000 3.186321382077705 0.000037000000000 3.181029079029593 0.000037100000000 3.175739150774377 0.000037200000000 3.170451594907418 0.000037300000000 3.165166409025479 0.000037400000000 3.159883590726722 0.000037500000000 3.154603137610705 0.000037600000000 3.149325047278384

[Alina Abdikadyr]

Thank you very much! But If I will have more h values, so may I display the deltaP values that are almost equal to 3?

[Torsten]

Add the lines

idx = deltaP > 2.8e5 & deltaP < 3.2e5;  % select those deltaP with 2.8e5 <= deltaP < = 3.2e5
[hnum(idx) deltaP(idx)]                 % Show these values together with the corresponding h values

at the end of the code.