% 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 = 1.225;

% prumer vrtule [m]
D = 2.667;

% rozsah rychlosti pro vypocet [m/s]
V = 0:10:200;

%==============================================================
% 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;
    % dle Alexandrova na str. 175 úplnì dole (zde ve vzorci je vytknutá
    % "2")
    Ph = 2*rho*S*(V(ii)^2*dvh*(1/2) + 2*V(ii)*dvh^2*(1/4) + dvh^3*(1/8));
    while Ph<Pmot
        dvh = 2*dvh;
        Ph = 2*rho*S*(V(ii)^2*dvh*(1/2) + 2*V(ii)*dvh^2*(1/4) + dvh^3*(1/8));
    end
    % vlastni numericky vypocet
    PocetIteraci = 20;
    dvl = 0;
    Pl = 2*rho*S*(V(ii)^2*dvl*(1/2) + 2*V(ii)*dvl^2*(1/4) + dvl^3*(1/8));
    Ph = 2*rho*S*(V(ii)^2*dvh*(1/2) + 2*V(ii)*dvh^2*(1/4) + dvh^3*(1/8));
    for jj = 1:PocetIteraci
        dvm = (dvh+dvl)/2;
        Pm = 2*rho*S*(V(ii)^2*dvm*(1/2) + 2*V(ii)*dvm^2*(1/4) + dvm^3*(1/8));
        
        if (Pm - Pmot) < 0
            Pl = Pm;
            dvl = dvm;
        else
            Ph = Pm;
            dvh = dvm;
        end
        dv = dvm;
    end
    
    % vypocet tahu dle Alexandrova. rce (91) str. 176 (opìt vytknutá "2")
    T(ii) = rho*S*(V(ii)*(1/2)+dv*(1/4))*2*dv;
    
end

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