% Constantes
m = 2.6;  % Massa em g
L = 1.0;  % Comprimento em m
r = 0.003; % Raio em m
g = 9.81; % Aceleração da gravidade em m/s^2
rho = 1.2; % Massa específica do ar em kg/m^3
cd = 0.1;  % Coeficiente de arrasto
A = pi * r^2; % Área transversal da esfera

% Calcula o coeficiente de arrasto aerodinâmico
b = 1/2 * rho * cd * A;

% Valores iniciais
theta0 = 0.05;  % Ângulo inicial em rad
w0 = 0.0;       % Velocidade angular inicial em rad/s
t0 = 0.0;       % Tempo inicial em s
h = 0.1;        % Passo de tempo em s

% Listas para armazenar os resultados
t = t0;
theta = theta0;
w = w0;
k1x = 0;
k1v = 0;
k2x = 0;
k2v = 0;

% Loop de tempo
while t(end) < 100.0
    % Cálculo dos valores intermediários
    k1x = w(end);
    k1v = -sign(w(end)) * (((b*L)/m) * w(end)^2) - (g / L) * theta(end);
    k2x = w(end) + k1v * h;
    k2v = -sign(w(end) + k1v h) * (((b*L)/m ) (w(end) + k1v *h)^2) - (g / L) * (theta(end)+k1x(end)*h);

    % Atualização dos valores de θ e w
    theta = [theta, theta(end) + ((k1x + k2x) / 2.0)*h];
    w = [w, w(end) + ((k1v + k2v) / 2.0)*h];

    % Atualização do tempo
    t = [t, t(end) + h];
end



% Cria o gráfico para θ com a linha vermelha
plot(t, theta, 'r');
% Adiciona uma linha para w com a linha azul
hold on;
plot(t, w, 'b');
% Adiciona uma legenda
legend('Theta', 'w');
% Adiciona títulos aos eixos
xlabel('tempo (s)');
ylabel('Theta (rad), w (rad/s)');
% Exibe o gráfico
hold off;