close all;
clear all;
clc;
%---------------------------------------------------------
%VSTUPNÍ PARAMETRY

B = 2;                      % [-] poèet listù 
R = 1.6/2;                  % [m] polomìr listu 
ro = 1.225;                 % [kg/m3] hustota vzduchu
mi = 1.59e-5;               % [] viskozita vzduchu
V = 55;                     % [m/s] max. rychlost letu 
Vz = 331.4;                 % [m/s] rychlost zvuku
omega = 2300*pi/30;         % [rad/s] uhlova rychlost;
P = 15000;                  % [W] výkon 

clD = 0.5;                  % [-] návrhový souèinitel vztlaku
alfaD = 1;                  % [°] návrhový úhel nábìhu
k = 0.0125;                 % k=cd/cl [-] pomìr odporu ku vztlaku ;
%-------------------------------------------------------------------
%VÝPOÈET POMÌRNÝCH VELIÈIN


r = 0.008:0.008:R;          % [m] polomìr elementu;
xi = r/R;                   % [-] pomìrný polomìr;
x = (omega.*r)./V;          % [-] pøevrácená hodnota rychlostního pomìru
lambda = V/(omega*R);       % [-] rychlostní pomìr
%---------------------------------------------------------------
%VÝPOÈET CIRKULACE


f = (B/2).*(((lambda^2+1)^0.5)./lambda).*(1.-xi);
% [-] souèinitel Prandtlovy ztrátové fukce
F = (2/pi).*acos(exp(-f));              % [-] Prandtlova ztrátová funkce
G = (F.*x.^2)./((x.^2)+1);              % [-] bezrozmìrná cirkulace

figure                                  % graf cirkulace funkcí pomìrného polomìru
plot(xi,G,'k')
grid on
xlabel('\xi [-]')
ylabel('G [-]')
axis ([0 1 0 1])
%--------------------------------------------------------------
%VÝPOÈET INTEGRANTÙ ZATÍŽENÍ


i1 = 4.*G.*xi.*(1.-k./x);
% integrant integrálu zatížení 

i2 = 2.*G.*(1.-k./x).*(1./(x.^2.+1)).*xi;
% integrant integrálu zatížení 

j1 = 4.*G.*xi.*(1+k.*x);
% integrant integrálu zatížení 

j2 = 2.*G.*(1.+k.*x).*(x.^2./(x.^2.+1)).*xi;
% integrant integrálu zatížení 

figure                                  %graf fce integrantù v závislosti na pomìrném polomìru
plot(xi,i1/4,'r--', xi,i2/2,'r',xi,j1/4,'b--',xi,j2/2,'b')
grid on
xlabel('\xi [-]')
ylabel('i1/4, i2/2, j1/4, j2/2 [-]')
legend('i1/4', 'i2/2', 'j1/4', 'j2/2')
axis ([0 1 0 1])

%---------------------------------------------------------------
%INTEGRACE OBDELNÍKOVOU METODOU POMOCÍ FOR-CYKLU


N=length(xi);   % velikost vektoru xi

I1 = 0;         % pøíprava
for i=1:N-1
    I1 = I1 + (xi(i+1)-xi(i))*i1(i);  
    % pøidání 1 obdelníku
end
I1;

I2 = 0;         % pøíprava
for i=1:N-1
    I2 = I2 + (xi(i+1)-xi(i))*i2(i);
    % pøidání 1 obdelníku
end
I2;

J1 = 0;         % pøíprava
for i=1:N-1
    J1 = J1 + (xi(i+1)-xi(i))*j1(i);
    % pøidání 1 obdelníku
end
J1;

J2 = 0;         % pøíprava
for i=1:N-1
    J2 = J2 + (xi(i+1)-xi(i))*j2(i);
    % pøidání 1 obdelníku
end
J2;
%---------------------------------------------------------------
%VÝPOÈET PØES MAX. VÝKON MOTORU


Pc = 2*P/(ro*V^3*pi*R^2);           % [-] koeficient výkonu

zeta = (J1./(2.*J2)).*(((1.+(4.*Pc.*J2./J1.^2)).^0.5)-1);
% [-] pomìr hybné rychlosti ku rychlosti letu

Tc = I1.*zeta - I2.*zeta.^2;        % [-] koeficient tahu
Pc = J1.*zeta + J2.*zeta.^2;        % [-] koeficient výkonu

T = Tc*ro*pi*V^2*R^2/2;             % [N] tah

eta = Tc./Pc;                       % [-] úèinnost

fi = atan((lambda./xi).*(1.+zeta./2))*180/pi;
% [°] rychlostní nastavení
beta = fi + alfaD;                  % [°] úhel nastavení 
zkrouceni = beta-beta(75);          % [°] zkroucení listu

WkuV = (x.^2.+1.-(zeta.*cos(fi.*pi./180)./2).^2).^0.5;
% [-] pomìr výsledné rychlosti a rychlosti letu
W = WkuV.*V;                        % [m/s] výsledná rychlost vrtulového listu

ckuR = (4.*pi.*lambda./B).*(G./WkuV).*(zeta./clD);
% [-] pomìrná šíøka listu 
c = ckuR.*R;                        % [m] šíøka listu

Re = W.*c.*ro./mi;                  % [-] Reynoldsovo èíslo
Ma = (W./Vz);                       % [-] Machovo èíslo

figure                              % graf šíøky listu v závislosti na pomìrném polomìru
plot(xi,c.*1000,'k')
grid on
xlabel('\xi [-]')
ylabel('c [mm]')

figure                              % graf zkroucení listu v závislosti na pomìrném polomìru
plot(xi,zkrouceni,'k')
grid on
xlabel('\xi [-]')
ylabel('\beta - \beta75 [°]')

figure                              % graf úhlu nastavení v závislosti na pomìrném polomìru
plot(xi,beta,'k')
grid on
xlabel('\xi [-]')
ylabel('\beta [°]')
axis([0 1 0 90])

%-------------------------------------------------------------

%VYSLEDKY
D=2*R
B
eta
c
beta
W
Re
Ma
