% skript pro vypocet tahu idealniho propulzoru
% neuvazuje se vliv rotace proudu vzduchu
% zadani v zakladnich jednotkach SI

clear all
close all

%==============================================================
% zadani
%==============================================================

% vykon motoru [W]
Pmot = 783e3;

% hustota vzduchu [kg/m3]
rho = 0.653;

% prumer vrtule [m]
D = 2.667;

% rozsah rychlosti pro vypocet [m/s]
V = 111.5:0.1:112.5;

%==============================================================
% pomocne operace
S = pi*D^2/4;
T = zeros(size(V));

for ii = 1:length(V)
    % numericky vypocet urychleni vzduchu dv v rovine vrtule
    % pouzita metoda puleni intervalu
    % odhad horni meze
    dvh = 10;
    Ph = 2*rho*S*(V(ii)^2*dvh + 2*V(ii)*dvh^2 + dvh^3);
    while Ph<Pmot
        dvh = 2*dvh;
        Ph = 2*rho*S*(V(ii)^2*dvh + 2*V(ii)*dvh^2 + dvh^3);
    end
    % vlastni numericky vypocet
    PocetIteraci = 20;
    dvl = 0;
    Pl = 2*rho*S*(V(ii)^2*dvl + 2*V(ii)*dvl^2 + dvl^3);
    Ph = 2*rho*S*(V(ii)^2*dvh + 2*V(ii)*dvh^2 + dvh^3);
    for jj = 1:PocetIteraci
        dvm = (dvh+dvl)/2;
        Pm = 2*rho*S*(V(ii)^2*dvm + 2*V(ii)*dvm^2 + dvm^3);
        
        if (Pm - Pmot) < 0
            Pl = Pm;
            dvl = dvm;
        else
            Ph = Pm;
            dvh = dvm;
        end
        dv = dvm;
    end
    
    % vypocet tahu
    T(ii) = rho*S*(V(ii)+dv)*2*dv;
    
end

% graf
figure
plot(V,T)
grid on
xlabel('Rychlost letu [m/s]')
ylabel('Tah [N]')

T = 112.5