clc; clear;

%funkční_kosnt.
p0 = 101325;            %Pa
rho_vzd0 = 1.225;       %kg/m3
g0 = 9.807;             %m/s2
Re = 6378000;           %m poloměr Země
M_He = 4.003*10^(-3);   %kg/mol molární hmotnost helia
M_vzd = 28.96*10^(-3);  %kg/mol    .   .   .    vzduchu
Rp = 8.314;             %J/kg/K univerzální plynová konst.

p11 = p0*(1-0.0065/288*11000)^(5.256);
rho_vzd11 = rho_vzd0*(1-0.0065/288*11000)^(4.256);

%letové hladiny
H = 23000;      %m výška horní 
D = 19000;      %m výška dolní 

%funkce (11-25) km
p_ok = @(h) p11*exp(-0.0001578*(h-11000));          %Pa
rho_vzd = @(h) rho_vzd11*exp(-0.0001578*(h-11000)); %kg/m3
g_h = @(h) g0*(Re/(Re+h))^2;                        %m\s2
rho_id = @(p, T, M) (p*M)/(Rp*T);                   %kg/m3

%tlaková bilance
dT = 15;           %K   superteplota vlivem infra záření ve dne, řekněme počasí
T_ok = 216.65;     %K   okolí
Td = T_ok + dT;    %K   celkem
f = 0.2;           % -  free lift (voleno)

ps_H = p_ok(H); %tlak okolí v horní 
g_H = g_h(H); %tíh. zrychení horní
rho_He_H = rho_id(ps_H, T_ok, M_He); 
rho_vzd_H = rho_vzd(H); 

pd_H = rho_He_H*Rp*Td/M_He;      %Pa tlak v balonu ve dne
dpf_H = f*ps_H;                  %Pa super-pressure od přebytku vztlaku
dpT_H = (dT/Td)*(ps_H + dpf_H);  %Pa superpressure od teploty
dp_H = dpf_H + dpT_H;            %celkový super-pressure 25km

%super pressure dole
ps_D = p_ok(D);                     %tlak okolí 19 km
g_D = g_h(D);                       %tíh. trychení 19 km
rho_He_D = rho_id(ps_D, T_ok, M_He); 
rho_vzd_D = rho_vzd(D);

pd_D = rho_He_D*Rp*Td/M_He;         %Pa tlak v balonu ve dne
dpf_D = f*ps_D;                     %Pa super-pressure od přebytku vztlaku
dpT_D = (dT/Td)*(ps_D + dpf_D);     %Pa superpressure od teploty
dp_D = dpf_D + dpT_D;               %celkový super-pressure dole

%rozměry, fixní hmotnosti
m_pay = 250;      %kg   payload
m_aku = 100;      %kg   akumulátor
t_f = 0.000038;   %m    tloušťka filmu 38 mikrometrů
t_p = 0.008;      %m    tloušťka portu
h_p = 0.2;        %m    šířka portu

%hustoty konstrukce balonu
rho_f = 924;      %kg/m3     film LLDPE
rho_t = 0.024;    %kg/m      tendon 6 mm 
rho_p = 2.2;      %kg/m      prut, uhlíkový kompozit
rho_l = 0.018;    %kg/m      lanoví konstrukce
p_E = 52;         %W/kg      elektrárna #soláry

%dovolená napětí
S_the = 30*10^6; %Pa dovolené namáhání ve folii
FT_max = 33300;  %N pevnost tendonu diag 6 mm
k_t = 2.25;      %bezpečnost dle předpisu
FD_t = FT_max/k_t; %dovolená síla jednoho tendonu
l_min = 0.2;     %m minimální délka tětivy goru na malém průměru toroidu

%počáteční hodnoty iterace
n_gor = 40;      %-    počet gorů
S = 10;          %m2     plocha filmu
R = 1;           %m      velký poloměr
r = 1;           %m      malý poloměr
V = 10;          %m3     objem He
k = 1.3;         %volený poměr mezi R/r = i
p_He = pd_D;     %tlak helia...
V_blst = 1000;   %objem balastu
m_c = 1;
m_blst = 1;

%iterační cyklus 
for i = 1:100
    
    %hmotnosti
    m_f = (1.2*S + (V_blst/r + 8*pi*r^2))*rho_f*t_f;  %hmotnost filmu včetně přepážek
    m_t = 2*pi*r*rho_t*n_gor;
    m_p = (16*R*sin(pi/8) + 5.2*r)*rho_p; %karbonový prut
    m_l = R/cos(pi/6)*rho_l;
    m_k = 8*m_l + m_p; %konstrukce
    m_b = m_f + m_t;
    m_He = rho_id(p_He, T_ok, M_He)*V;
    m_E = (m_c + m_blst)/p_E + m_aku; %elektrárna
    m_c = m_pay + m_He + m_b + m_E + m_k;

    %rádiusy, objem, plášť
    Fvz = m_c*g_H;                     %Fvz=F_G 
    V = Fvz/(rho_vzd_H*g_H);
    
    R = ((k^2)*V/(2*pi^2))^(1/3); 
    r = R/k;
    S = 4*pi^2*R*r;

    %balast
    m_blst = V*rho_vzd(D) - m_c;                    %kg hmotnost vzduchového balastu
    p_blst = (m_blst/M_vzd*Rp*T_ok)/(V);            %absolutní tlak balastu v balonu
    p_He = p_ok(D)+ dp_D - p_blst;                  %absolutní tlak He v balonu
    V_blst = (m_blst*Rp*T_ok)/((pd_D+dp_D)*M_vzd);  %maximální objem vzduchu uvnitř balonu

    %návrh goru CR, tendonů
    R_the = S_the*t_f/dp_D;                   %m max radius vypoukliny
    F_p = dp_D*pi*4*R*r;                      %celková síla přenášená tendony
    n_tmin = F_p/(2*FD_t);                    %minimální počet tendonů
    n_GminV = (2*pi*(R+r))/(2*R_the);         %ks minimální počet gorů z velkého průměru toroidu
    n_GmaxM = (2*pi)/(atan(l_min/(2*(R-r)))); %ks maximální počet goru z malého průměru toru
    n_G = round(max(n_tmin, n_GminV)+2);      %volba počtu gorů
    n_gor = n_G;
end

fprintf('R = %f m\n', R)
fprintf('r = %f m\n', r)
fprintf('V = %f m3\n', V)
fprintf('R_the = %f m\n', R_the)
fprintf('n_tmin = %f ks\n', n_tmin)
fprintf('n_GmaxM = %f ks\n', n_GmaxM)
fprintf('n_GminV = %f ks\n', n_GminV)
fprintf('n_gor = %f ks\n', n_gor)
fprintf('delka = %f m\n', 2*(R+r))

if p_He - p_ok(H)>dp_H && n_GmaxM - n_gor > 0
    disp ('Dobrý')
else 
    disp ('Není to dobrý')
end
