intersection from a plot and loop

Hello everyone, could you please help me to find the intersection of two graphs, figure 1

My code is:

clear all; close all

W_takeoff = 10000;
W_landing=6000;
S = 20;
AR = 5;
cd0 = 0.02;
k = 1/pi/AR;
RC=0.51;
clalpha = 2*pi;
amin=2;
astall=12;
rho=1;
ct=0.001;
figure(1);hold on; xlabel('V');ylabel('D')
    i=0;
    for alpha = amin:0.001:astall
        i=i+1;
        cl(i) = clalpha * alpha * pi/180;
        V_takeoff(i) = sqrt(2*W_takeoff/rho/S/cl(i));
        V_landing(i) = sqrt(2*W_landing/rho/S/cl(i));
        cd(i) = cd0 + k * cl(i) * cl(i);
        D_takeoff(i) = 0.5 * rho * V_takeoff(i) * V_takeoff(i) * S * cd(i);
        D_landing(i) = 0.5 * rho * V_landing(i) * V_landing(i) * S * cd(i);
        p_takeoff(i) = D_takeoff(i)*V_takeoff(i);
        p_landing(i) = D_landing(i)*V_landing(i);
        P_takeoff(i)=20000+500*V_takeoff(i);
        T_takeoff(i)=20000/V_takeoff(i)+500;
        P_landing(i)=20000+500*V_landing(i);
        T_landing(i)=20000/V_landing(i)+500;
        cl_cd2(i)=cl(i)/(cd(i)*cd(i));
        ang(i)=alpha;
    end  
    figure(1); plot(V_takeoff,D_takeoff)
    hold on 
    plot(V_takeoff,T_takeoff)
    title(['Takeoff'])

Answer 1

[Antoni Garcia-Herreros]

Hello Alina,

You can check the fzero function

clear all; close all
W_takeoff = 10000;
W_landing=6000;
S = 20;
AR = 5;
cd0 = 0.02;
k = 1/pi/AR;
RC=0.51;
clalpha = 2*pi;
amin=2;
astall=12;
rho=1;
ct=0.001;
figure(1);hold on; xlabel('V');ylabel('D')
    i=0;
    for alpha = amin:0.001:astall
        i=i+1;
        cl(i) = clalpha * alpha * pi/180;
        V_takeoff(i) = sqrt(2*W_takeoff/rho/S/cl(i));
        V_landing(i) = sqrt(2*W_landing/rho/S/cl(i));
        cd(i) = cd0 + k * cl(i) * cl(i);
        D_takeoff(i) = 0.5 * rho * V_takeoff(i) * V_takeoff(i) * S * cd(i);
        D_landing(i) = 0.5 * rho * V_landing(i) * V_landing(i) * S * cd(i);
        p_takeoff(i) = D_takeoff(i)*V_takeoff(i);
        p_landing(i) = D_landing(i)*V_landing(i);
        P_takeoff(i)=20000+500*V_takeoff(i);
        T_takeoff(i)=20000/V_takeoff(i)+500;
        P_landing(i)=20000+500*V_landing(i);
        T_landing(i)=20000/V_landing(i)+500;
        cl_cd2(i)=cl(i)/(cd(i)*cd(i));
        ang(i)=alpha;
    end  
    figure(1); plot(V_takeoff,D_takeoff)
    hold on 
    plot(V_takeoff,T_takeoff)
    title(['Takeoff'])
    
    pp = spline(V_takeoff,D_takeoff); %Fit curve 1 to spline
    pp2 = spline(V_takeoff,T_takeoff); %Fit curve 2 to spline
    x1=xlim;
    x0 = fzero(@(x) (ppval(pp,x)-ppval(pp2,x)), 50); % Find the root of the difference between the two curves
    line([x1(1) x0],[ppval(pp,x0) ppval(pp,x0)],'LineStyle','--','Color','k')
    y1=ylim;
    line([x0 x0],[y1(1) ppval(pp,x0)],'LineStyle','--','Color','k')

Answer 2

[Star Strider]

clear all; close all

W_takeoff = 10000;
W_landing=6000;
S = 20;
AR = 5;
cd0 = 0.02;
k = 1/pi/AR;
RC=0.51;
clalpha = 2*pi;
amin=2;
astall=12;
rho=1;
ct=0.001;
% figure(1);hold on; xlabel('V');ylabel('D')
    i=0;
    for alpha = amin:0.001:astall
        i=i+1;
        cl(i) = clalpha * alpha * pi/180;
        V_takeoff(i) = sqrt(2*W_takeoff/rho/S/cl(i));
        V_landing(i) = sqrt(2*W_landing/rho/S/cl(i));
        cd(i) = cd0 + k * cl(i) * cl(i);
        D_takeoff(i) = 0.5 * rho * V_takeoff(i) * V_takeoff(i) * S * cd(i);
        D_landing(i) = 0.5 * rho * V_landing(i) * V_landing(i) * S * cd(i);
        p_takeoff(i) = D_takeoff(i)*V_takeoff(i);
        p_landing(i) = D_landing(i)*V_landing(i);
        P_takeoff(i)=20000+500*V_takeoff(i);
        T_takeoff(i)=20000/V_takeoff(i)+500;
        P_landing(i)=20000+500*V_landing(i);
        T_landing(i)=20000/V_landing(i)+500;
        cl_cd2(i)=cl(i)/(cd(i)*cd(i));
        ang(i)=alpha;
    end  
    
    Takeoffc = {D_takeoff; p_takeoff; P_takeoff; T_takeoff};
    Landingc = {D_landing; p_landing; P_landing; T_landing};
    
    L = numel(V_takeoff);
    xidx = find(diff(sign(D_takeoff-T_takeoff)));                                           % Approximate Index Of Intersection
    for k1 = 1:numel(xidx)
        idxrng = max(xidx(k1)-1,1) : min(xidx(k1)+1,L);
        V_t(k1,:) = interp1(D_takeoff(idxrng)-T_takeoff(idxrng), V_takeoff(idxrng), 0);
        X_t(k1,:) = interp1(V_takeoff(idxrng), D_takeoff(idxrng), V_t(k1));
    end
    
    figure; 
    plot(V_takeoff,D_takeoff)
    hold on 
    plot(V_takeoff,T_takeoff)
    xl = xlim;
    yl = ylim;
    for k = 1:numel(xidx)
        hp = plot(V_t(k), X_t(k), 'sm', 'MarkerSize',10);
        plot([xl(1) V_t(k)], [1 1]*X_t(k), ':', 'Color',hp.Color)
        plot([1 1]*V_t(k), [yl(1) X_t(k)], ':', 'Color',hp.Color)
    end
    title(['Takeoff'])
    xlabel('V')
    ylabel('D')
    text(V_t,X_t,compose('    (%.2f, %.2f)',[V_t X_t]), 'Horiz','left', 'Vert','middle')

I suspect that you want to do more than one of these comparisons. I created two cell arrays for the takeoffs and landings to make this easier. All that is necessary is to index into them, and then make the appropriate substitutions in the loops. So comparing ‘D_takeoff’ and ‘T_takeoff’ would only involve referencing ‘Takeoffc{1}’ and ‘Takeoffc{4}’ . To make this more efficient, build a matrix of the comparisons you want to make and then just index into it, so here the first row might be [1 4] and whatever other comparisons you want to make. The code shoud adapt eassily to it, and it is also set up to text for more than one intersection in each comparison, if necessary.

EDIT — (31 Mar 2023 at 15:56)

That might go something like this —

clear all; close all

W_takeoff = 10000;
W_landing=6000;
S = 20;
AR = 5;
cd0 = 0.02;
k = 1/pi/AR;
RC=0.51;
clalpha = 2*pi;
amin=2;
astall=12;
rho=1;
ct=0.001;
% figure(1);hold on; xlabel('V');ylabel('D')
    i=0;
    for alpha = amin:0.001:astall
        i=i+1;
        cl(i) = clalpha * alpha * pi/180;
        V_takeoff(i) = sqrt(2*W_takeoff/rho/S/cl(i));
        V_landing(i) = sqrt(2*W_landing/rho/S/cl(i));
        cd(i) = cd0 + k * cl(i) * cl(i);
        D_takeoff(i) = 0.5 * rho * V_takeoff(i) * V_takeoff(i) * S * cd(i);
        D_landing(i) = 0.5 * rho * V_landing(i) * V_landing(i) * S * cd(i);
        p_takeoff(i) = D_takeoff(i)*V_takeoff(i);
        p_landing(i) = D_landing(i)*V_landing(i);
        P_takeoff(i)=20000+500*V_takeoff(i);
        T_takeoff(i)=20000/V_takeoff(i)+500;
        P_landing(i)=20000+500*V_landing(i);
        T_landing(i)=20000/V_landing(i)+500;
        cl_cd2(i)=cl(i)/(cd(i)*cd(i));
        ang(i)=alpha;
    end  
    
    Takeoffc = {D_takeoff; p_takeoff; P_takeoff; T_takeoff};
    Landingc = {D_landing; p_landing; P_landing; T_landing};
    Lbl = {'D','p','P','T'};
    Pairs = [1 4; 2 3];
    
    L = numel(V_takeoff);
    for k2 = 1:size(Pairs,1)
        V1 = Takeoffc{Pairs(k2,1)};
        V2 = Takeoffc{Pairs(k2,2)};
        xidx = find(diff(sign(V1-V2)));                                                                 % Approximate Index Of Intersection
        for k1 = 1:numel(xidx)
            idxrng = max(xidx(k1)-1,1) : min(xidx(k1)+1,L);
            V_t(k1,:) = interp1(V1(idxrng)-V2(idxrng), V_takeoff(idxrng), 0);
            X_t(k1,:) = interp1(V_takeoff(idxrng), V1(idxrng), V_t(k1));
        end
        
        figure; 
        plot(V_takeoff,V1)
        hold on 
        plot(V_takeoff,V2)
        xl = xlim;
        yl = ylim;
        for k = 1:numel(xidx)
            hp = plot(V_t(k), X_t(k), 'sm', 'MarkerSize',10);
            plot([xl(1) V_t(k)], [1 1]*X_t(k), ':', 'Color',hp.Color)
            plot([1 1]*V_t(k), [yl(1) X_t(k)], ':', 'Color',hp.Color)
        end
        title(['Takeoff'])
        xlabel('V')
        ylabel(Lbl{Pairs(k2,1)})
        text(V_t,X_t,compose('    (%.2f, %.2f)',[V_t X_t]), 'Horiz','left', 'Vert','middle')
    end

.