import numpy as np
import matplotlib.pyplot as plt
from lagrange import maxminv


#Exercício 3:  Variação do período com comprimento do pêndulo
#a) Faça um gráfico do logaritmo do período 𝑇 em funçaõ do logaritmo do comprimento 𝐿.


def acelera_pendulo(theta, L, g=9.8):
    return -g / L * np.sin(theta)


def euler_cromer(theta0, v0, L, g, dt, t):
    n = len(t)
    theta = np.zeros(n)
    v = np.zeros(n)
    theta[0] = theta0
    v[0] = v0
    for i in range(1, n):
        a = acelera_pendulo(theta[i-1], L, g)
        v[i] = v[i-1] + a * dt
        theta[i] = theta[i-1] + v[i] * dt
    return theta

def medir_periodo(theta, t):
    maximos_tempos = []
    for i in range(1, len(theta)-1):
        if theta[i-1] < theta[i] and theta[i] > theta[i+1]:
            t_max, _ = maxminv(t[i-1], t[i], t[i+1], theta[i-1], theta[i], theta[i+1])
            maximos_tempos.append(t_max)
    if len(maximos_tempos) >= 2:
        return maximos_tempos[1] - maximos_tempos[0]
    else:
        return None


g = 9.8
theta0 = 0.1
v0 = 0
dt = 0.001
t = np.arange(0, 10, dt)

# Comprimentos a testar
L_values = np.array([0.5, 0.8, 1.0, 1.5, 2.0, 3.0])
T_values = []

# Calcula períodos para cada L
for L in L_values:
    theta = euler_cromer(theta0, v0, L, g, dt, t)
    T = medir_periodo(theta, t)
    if T:
        T_values.append(T)
    else:
        T_values.append(np.nan)

# Remove Nans
L_valid = L_values[~np.isnan(T_values)]
T_valid = np.array(T_values)[~np.isnan(T_values)]

# Gráfico log-log
logL = np.log10(L_valid)
logT = np.log10(T_valid)

# Ajuste linear manual
N = len(logL)
sum_x = np.sum(logL)
sum_y = np.sum(logT)
sum_xx = np.sum(logL**2)
sum_xy = np.sum(logL * logT)

# Coeficientes
a = (N * sum_xy - sum_x * sum_y) / (N * sum_xx - sum_x**2)
b = (sum_y - a * sum_x) / N

# Erro do declive
residuals = logT - (a * logL + b)
sigma_y = np.sqrt(np.sum(residuals**2) / (N - 2))
sigma_a = sigma_y / np.sqrt(np.sum((logL - np.mean(logL))**2))

# Gráfico
plt.figure(figsize=(8, 5))
plt.plot(logL, logT, 'o', label='Dados')
plt.plot(logL, a * logL + b, '-', label=f'Ajuste: slope={a:.2f}±{sigma_a:.2f}')
plt.xlabel('log10(L)')
plt.ylabel('log10(T)')
plt.title('log(T) vs log(L)')
plt.grid()
plt.legend()
plt.show()

# Resultado
print(f"Declive (slope): {a:.4f} ± {sigma_a:.4f}")
print("O valor esperado para o declive é 0.5, pois T ∝ √L.")