vys = input('Zadaj vysku v metroch: ');
vzd = input('Zadaj vzdialenost v metroch: ');

% prealokacia
a_max_vys=zeros(4,1);
a_max_vzd=zeros(4,1);
a_vys=cell(1,1);
a_vzd=cell(1,1);
F_vys_a=cell(1,1);
m_vys_a=cell(1,1);
m_vys_1_a=cell(1,1);
F_vys_b=cell(1,1);
m_vys_b=cell(1,1);
m_vys_1_b=cell(1,1);

F_vzd_a=cell(1,1);
F_vzd_h_a=cell(1,1);
m_vzd_a=cell(1,1);
m_vzd_1_a=cell(1,1);
F_vzd_b=cell(1,1);
F_vzd_h_b=cell(1,1);
m_vzd_b=cell(1,1);
m_vzd_1_b=cell(1,1);

F_kl_a=cell(1,1);
m_kl_a=cell(1,1);
m_kl_1_a=cell(1,1);
F_kl_b=cell(1,1);
m_kl_b=cell(1,1);
m_kl_1_b=cell(1,1);

t_vys_a=zeros(4,1);
t_vzd_a=zeros(4,1);
t_kl_a=zeros(4,1);

t1_s=zeros(4,1);
t2_s=zeros(4,1);
t1_l=zeros(4,1);
t2_l=zeros(4,1);
t1_kl=zeros(4,1);
t2_kl=zeros(4,1);


I_1_vys_a=cell(1,1);
I_vys_a=cell(1,1);
I_1_vys_b=cell(1,1);
I_vys_b=cell(1,1);
I_1_vzd_a=cell(1,1);
I_vzd_a=cell(1,1);
I_1_vzd_b=cell(1,1);
I_vzd_b=cell(1,1);
I_1_kl_a=cell(1,1);
I_kl_a=cell(1,1);
I_1_kl_b=cell(1,1);
I_kl_b=cell(1,1);


C_vystup=zeros(4,1);
C_let=zeros(4,1);
C_klesanie=zeros(4,1);
C_nova=zeros(4,1);
C_nova2=cell(1,1);
Cas2=cell(1,1);

t1_s=cell(1,1);
t2_s=cell(1,1);
t1_l=cell(1,1);
t2_l=cell(1,1);
t1_kl=cell(1,1);
t2_kl=cell(1,1);

j=0;

for v_d = 1:1:4
	j = j+1;
	p_m = 4;  %pocet motorov
	min = 60;   %pocet sekund za minutu
	S = 0.052572;   %plocha UAV
	g = 9.81;   %gravitacne zrychlenie
	cd_v = 0.09;    %koeficient odporu vo vertikalnom smere
    cd_h = 0.05;    %koeficient odporu v horizontalnom smere
	Ro = 1.225;     %hustota vzduchu na jeden meter stvorcovy
	p = 100;    %100%
	C = 240;    %kapacita baterie v amperminutach
	m = 1583;   %hmotnost UAV v gramoch
    max_m = 0.72;   %maximalny tah jedneho motora v kilogramoch
    a =(p_m*max_m*g-(m/1000)*g-(Ro/2)*cd_v*S*(v_d)^2)/(m/1000); %maximalne zrychlenie motorov pri vahe 720g na jeden motor
    
    t1_s{j} = 0:0.01:(v_d/a);   %cas, za ktory sa UAV dostane na pozadovanu ryhlost pocas stupania
    t2_s{j} = 0:0.01:(vys/v_d - t1_s{j}(end));  %cas, ktory je potrebny na vystupanie do zadanej vysky s pozadovanou rychlostou pocas stupania
    
    t1_l{j} = 0:0.01:(v_d/a);   %cas, za ktory sa UAV dostane na pozadovanu ryhlost pocas letu
    t2_l{j} = 0:0.01:(vzd/v_d - t1_l{j}(end));  %cas, ktory je potrebny na prekonanie zadanej vdzialenosti s pozadovanou rychlostou pocas letu    
    
    t1_kl{j} = 0:0.01:(v_d/a);  %cas, za ktory sa UAV dostane na pozadovanu ryhlost pocas klesania
    t2_kl{j} = 0:0.01:(vys/v_d - t1_kl{j}(end));    %cas, ktory je potrebny na pristatie zo zadanej vysky s pozadovanou rychlostou pocas klesania  
    
    F_vys_a{j} = ((m/1000)*a+(m/1000)*g+(Ro/2)*cd_v*S*(a*t1_s{j}).^2);  %tahova sila UAV pri zrychlovani na pozadovanu rychlost pocas stupania
    F_vys_b{j} = ((m/1000)*g+(Ro/2)*cd_v*S*(v_d)^2)*ones(1,length(t2_s{j})); %tahova sila UAV pri pohybe s konstantnou rychlostou pocas stupania a vyrobenie vektora so samymi jednotkami o dlzke t2_s
    
    m_vys_a{j} = F_vys_a{j}/g;  %hmotnost UAV pri zrychlovani na pozadovanu rychlost pocas stupania
    m_vys_b{j} = F_vys_b{j}/g;  %hmotnost UAV pri pohybe s konstantnou rychlostou pocas stupania
    
    m_vys_1_a{j}= m_vys_a{j}/p_m;   %hmotnost na jeden motor pri zrychlovani na pozadovanu rychlost pocas stupania
    m_vys_1_b{j}= m_vys_b{j}/p_m;   %hmotnost na jeden motor pri pohybe s konstantnou rychlostou pocas stupania

    F_vzd_h_a{j} = ((m/1000)*a+(Ro/2)*S*cd_h*(a*t1_l{j}).^2);   %tahova sila UAV pri zrychlovani na pozadovanu rychlost pocas pohybu v horizontalnom smere
    F_vzd_h_b{j} = ((Ro/2)*S*cd_h*v_d^2)*ones(1,length(t2_l{j}));   %tahova sila UAV pri pohybe s konstantnou rychlostou pocas pohybu v horizontalnom smere a vyrobenie vektora so samymi jednotkami o dlzke t2_l
    F_vzd_v = (m/1000)*g;   %tahova sila UAV pri pohybe vo vertikalnom smere 
    
    F_vzd_a{j} = sqrt((F_vzd_v).^2 + (F_vzd_h_a{j}).^2);	%vysledna tahova sila UAV pri zrychlovani na pozadovanu rychlost pocas letu
    F_vzd_b{j} = sqrt((F_vzd_v).^2 + (F_vzd_h_b{j}).^2);    %vysledna tahova sila UAV pri pohybe s konstantnou rychlostou pocas letu
    
    m_vzd_a{j} = F_vzd_a{j}/g;  %hmotnost UAV pri zrychlovani na pozadovanu rychlost pocas letu
    m_vzd_b{j} = F_vzd_b{j}/g;  %hmotnost UAV pri pohybe s konstantnou rychlostou pocas letu
    
    m_vzd_1_a{j}= m_vzd_a{j}/p_m;   %hmotnost na jeden motor pri zrychlovani na pozadovanu rychlost pocas letu
    m_vzd_1_b{j}= m_vzd_b{j}/p_m;   %hmotnost na jeden motor pri pohybe s konstantnou rychlostou pocas letu
	
    %Sila je na malu chvilu zaporna z dovodu ze motory musia zrychlit na
    %pozadovanu rychlost a potom ju len udrzuju, pretoze volnym padom by
    %dosiahol pozadovanu rychlost v case t=v_d/g
    
    F_kl_a{j} = ((m/1000)*g-(m/1000)*a-(Ro/2)*cd_v*S*(a*t1_kl{j}).^2);  %tahova sila UAV pri zrychlovani na pozadovanu rychlost pocas klesania
    F_kl_b{j} = ((m/1000)*g-(Ro/2)*cd_v*S*v_d^2)*ones(1,length(t2_kl{j}));  %tahova sila UAV pri pohybe s konstantnou rychlostou pocas klesania a vyrobenie vektora so samymi jednotkami o dlzke t2_kl
    
    m_kl_a{j} = abs(F_kl_a{j})/g;   %hmotnost UAV pri zrychlovani na pozadovanu rychlost pocas klesania
    m_kl_b{j} = abs(F_kl_b{j})/g;   %hmotnost UAV pri pohybe s konstantnou rychlostou pocas klesania
    
    m_kl_1_a{j} = m_kl_a{j}/p_m;    %hmotnost na jeden motor pri zrychlovani na pozadovanu rychlost pocas klesania
    m_kl_1_b{j} = m_kl_b{j}/p_m;    %hmotnost na jeden motor pri pohybe s konstantnou rychlostou pocas klesania

	m_1 = m/p_m;    %celkova hmotnost podelena poctom motorov
 	I_1 = 0.0012 * (m_1^1.3685);    %urcenie prudoveho odberu jedneho motora v zavislosti na hmotnosti
	I=I_1*p_m;  %urcenie prudoveho odberu celeho UAV
	T = linspace(-10,25,100);   %naskalovanie vonkajsej teploty prostredia
	C_p = 0.8814 + 0.0091*T - 0.0001 * T.^2;    %urcenie straty kapacity na vonkajsej teplote prostredia v percentach

	I_1_vys_a{j} = 0.0012 * ((m_vys_1_a{j}*1000).^1.3685);  %prudovy odber jedneho motora pri zrychlovani pocas stupania
	I_1_vys_b{j} = 0.0012 * ((m_vys_1_b{j}*1000).^1.3685);  %prudovy odber jedneho motora pri pohybe konstantnou rychlostou pocas stupania
    
	I_vys_a{j} = I_1_vys_a{j}*p_m;  %prudovy odber celeho UAV pri zrychlovani pocas stupania 
	I_vys_b{j} = I_1_vys_b{j}*p_m;  %prudovy odber celeho UAV pri pohybe konstantnou rychlostou pocas stupania

	I_1_vzd_a{j} = 0.0012 * ((m_vzd_1_a{j}*1000).^1.3685);  %prudovy odber jedneho motora pri zrychlovani pocas letu
	I_1_vzd_b{j} = 0.0012 * ((m_vzd_1_b{j}*1000).^1.3685);  %prudovy odber jedneho motora pri pohybe konstantnou rychlostou pocas letu
    
	I_vzd_a{j} = I_1_vzd_a{j}*p_m;  %prudovy odber celeho UAV pri zrychlovani pocas letu
	I_vzd_b{j} = I_1_vzd_b{j}*p_m;  %prudovy odber celeho UAV pri pohybe konstantnou rychlostou pocas letu

	I_1_kl_a{j} = 0.0012 * ((m_kl_1_a{j}*1000).^1.3685);    %prudovy odber jedneho motora pri zrychlovani pocas klesania
	I_1_kl_b{j} = 0.0012 * ((m_kl_1_b{j}*1000).^1.3685);    %prudovy odber jedneho motora pri pohybe konstantnou rychlostou pocas klesania
    
	I_kl_a{j} = I_1_kl_a{j}*p_m;    %prudovy odber celeho UAV pri zrychlovani pocas klesania
	I_kl_b{j} = I_1_kl_b{j}*p_m;    %prudovy odber celeho UAV pri pohybe konstantnou rychlostou pocas klesania
    
    % lichobeznikova metoda integrovania, integrujeme prud podla casu   
	C_vystup(j)=trapz(t1_s{j},I_vys_a{j})/min+trapz(t2_s{j},I_vys_b{j})/min;
	C_let(j)=trapz(t1_l{j},I_vzd_a{j})/min+trapz(t2_l{j},I_vzd_b{j})/min;
	C_klesanie(j)=trapz(t1_kl{j},I_kl_a{j})/min+trapz(t2_kl{j},I_kl_b{j})/min;
    
	C_nova(j) = C - (C_vystup(j) + C_let(j) + C_let(j) + C_klesanie(j)); %odcitanie straty kapacit pocas manevrov od maximalnej kapacity akumulatora
	C_nova2{j} = C_nova(j)-C_nova(j)*(1-C_p); %urcenie vyslednej kapacity v zavislosti na teplote prostredia, ktora sa vyuzije pri viseni UAV

	Cas2{j} = (C_nova2{j}./I);
end
figure()
plot(T,Cas2{1},'r')
hold on
plot(T,Cas2{2},'g')
plot(T,Cas2{3},'b')
plot(T,Cas2{4},'k')
xlabel('Teplota vonkajšieho prostredia [°C]')
ylabel('Vydrz baterie [min]')
title('Vydrz baterie v zavislosti na teplote')
legend('1m/s','2m/s','3m/s','4m/s')
grid on;