pos_fig = figure;
radi_fig = figure;
vel_fig = figure;
accel_fig = figure;
for ro = 0.1:0.5:5
snowball_sim(ro, 100, 0.05, pos_fig, radi_fig, vel_fig, accel_fig)
end
function snowball_sim(ro, h, C, pos_fig, radi_fig, vel_fig, accel_fig)
% Set the constants
g = 9.81;
density = 600;
k1 = 30/7;
k2 = (5 * g * C) / (14 * pi); % Assuming g = 9.81 m/s^2
% Set the initial conditions
initial_conditions = [ro; 0];
% Set the time span for the solution
tspan = [0 30]; % Adjust the end time as needed
% Set the maximum step size
maxStep = 0.1; % Adjust this value as needed
% Solve the system of ODEs using ode45
options = odeset('MaxStep', maxStep);
[t, y] = ode45(@(t, y) myODE(t, y, k1, k2), tspan, initial_conditions, options);
% Extract the results
r_t = y(:, 1);
% Calculate y(t)
y_t = ((2 * pi) / C) * (r_t.^2 / 2 - ro^2 / 2);
% Calculate the velocity u(t)
u_t = ((2 * pi) / C) * r_t .* y(:, 2);
% Differentiate the velocity to get acceleration a(t)
a_t = diff(u_t) ./ diff(t);
% Plot individual figures
figure(pos_fig);
% Figure 1: Snowball's vertical position y
plot(t, y_t, 'LineWidth', 2);
xlabel('t(s)');
ylabel("y(t)");
title("Snowball's Vertical Position as a Function of Time");
grid on;
hold on;
figure(radi_fig);
% Figure 2: Snowball's Radius r
plot(t, r_t, 'LineWidth', 2);
xlabel('t(s)');
ylabel("r(t)");
title("Snowball's Radius as a Function of Time");
grid on;
hold on;
figure(vel_fig);
% Figure 3: Snowball's Velocity u
plot(t, u_t, 'LineWidth', 2);
xlabel('t(s)');
ylabel("u(t)");
title("Snowball's velocity as a Function of time");
grid on;
hold on;
figure(accel_fig);
% Snowball's acceleration a
plot(t(1:end-1), a_t, 'LineWidth', 2);
xlabel('t(s)');
ylabel('a(t)');
title("Snowball's Acceleration");
grid on;
hold on;
%Calculate time for snowball to hit road
indx = 1;
while r_t(indx) < sqrt((h*C)/pi + ro^2)
to = t(indx);
indx = indx+1;
end
disp("____________________________________________________")
disp("ro = "+ro +", C = " +C + ", h = " + h)
disp("____________________________________________________")
to
disp("----------------------------------------------------")
%Snowball's mass ratio
mass_ratio = (r_t(indx)/ro).^3
disp("----------------------------------------------------")
mass_final = (4/3)*pi*density*(r_t(indx))^3
disp("----------------------------------------------------")
%Snowball's total kinetic energy
K = (28*pi*density/15)*((r_t(indx))^3)*(u_t(indx))^2
% Define the function representing the system of ODEs
function dydt = myODE(t, y, k1, k2)
dydt = zeros(2, 1);
dydt(1) = y(2);
dydt(2) = (k2 - k1*y(2)^2) / y(1);
end
end
