% balloon floating simulation
%
% input: altitude vector
% output: resulting force vector 

clear all

% launch altitude in meters
launch_altitude = 300;
% free lift in kg
initial_free_lift = 0.002;

g = 9.81;

% altitude dependent data from ISA
alt = 1:9000;
% temperature in kelvin
% pressure in pascal
% density in kg / m^3
[T, ~, p, rho] = atmosisa(alt);

% initial density of helium and air
rho_he_init = p(launch_altitude) / (2077.058 * T(launch_altitude));
rho_air_init = rho(launch_altitude);

% balloon and payload weigths
m_balloon = 0.038;
m_payload = 0.0209;

% initial balloon volume, derived from F_A - F_G = F_free_lift
V_init = (initial_free_lift / 100) / ((rho_air_init - rho_he_init)*g) ...
    + (m_balloon + m_payload) / (rho_air_init - rho_he_init);

% helium filling weight
m_he = V_init * rho_he_init;

F_g = g*(m_balloon + m_payload + m_he);

% balloon volume over altitude
% assumption: pressure internal = pressure external
% p * V / T = const
const = p(launch_altitude) * V_init / T(launch_altitude);
% V = const * T / p
V_balloon(alt) = const * T(alt) ./ p(alt);

V_max = 0.125;
V_balloon = min(V_balloon, V_max);
F_a(alt) = g .* rho(alt) .* V_balloon(alt);

% calculate the resulting force
F_r = F_a - F_g;

% find zero crossing of the resulting force
tmp = diff(sign(F_r));
float_altitude = find(tmp<0)

hold on
[ax, p1, p2] = plotyy(alt, F_r, alt, V_balloon);
axis(ax(1), [0, float_altitude*1.1, -max(F_r)*0.5, max(F_r)*1.5]);
axis(ax(2), [0, float_altitude*1.1, 0, V_max * 1.2]);
xlabel('Altitude in m ASL');
ylabel(ax(1),'Resulting Force in N');
ylabel(ax(2),'Balloon Volume in m^3');
grid(ax(2), 'on');








