import numpy as np
import matplotlib.pyplot as plt

# Datos experimentales J.A. Joaqui, B.D. Tovar and H.L. Valencia
#tiempos = np.array([0.448, 0.587, 0.553, 0.485, 0.680, 0.629, 0.731, 0.705])
#distancias = np.array([0.97, 1.68, 1.50, 1.16, 2.25, 1.93, 2.62, 2.41])

#Esperimento PhyPhox
tiempos = np.array([0.382, 0.434, 0.488, 0.539, 0.586, 0.629, 0.66, 0.702])
distancias = np.array([0.70, 0.95, 1.2, 1.45, 1.70, 1.95, 2.2, 2.45])

# Calcular tiempos al cuadrado
tiempo_cuadrado = tiempos**2

# Ajuste lineal para altura vs. tiempo^2
coef_lineal, cov_lineal = np.polyfit(tiempo_cuadrado, distancias, 1, cov=True)
m, b = coef_lineal
error_m, error_b = np.sqrt(np.diag(cov_lineal))

# Calcular alturas predichas por el modelo lineal
distancias_predichas = m * tiempo_cuadrado + b

# Calcular R²
ss_total = np.sum((distancias - np.mean(distancias))**2)
ss_residual = np.sum((distancias - distancias_predichas)**2)
r2 = 1 - (ss_residual / ss_total)

# Crear la curva de ajuste
tiempo_ajuste = np.linspace(min(tiempo_cuadrado), max(tiempo_cuadrado), 100)
distancia_ajuste = m * tiempo_ajuste + b

# Gráfico altura vs. tiempo^2
plt.figure(figsize=(10, 5))
plt.scatter(tiempo_cuadrado, distancias, color='blue', s=50)
plt.plot(tiempo_ajuste, distancia_ajuste, color='red', linestyle='--', 
         label=f'Ajuste Lineal: $y = {m:.3f}x + {b:.3f}$\n $R^2 = {r2:.3f}$')
plt.xlabel('Tiempo al cuadrado (s$^2$)', fontsize=16)
plt.ylabel('Altura (m)', fontsize=16)
plt.xticks(fontsize=14)
plt.yticks(fontsize=14)
plt.legend(fontsize=12)
plt.grid(True)
plt.tight_layout()
plt.show()

# Generar tabla con repetición, altura, tiempo y tiempo²
print("Toma, Altura, Tiempo, Tiempo al cuadrado:")
print("="*42)
print(f"{'Toma':<8}{'h(m)':<12}{'t(s)':<12}{'t²(s²)':<12}")
print("="*42)
for i, (h, t, t2) in enumerate(zip(distancias, tiempos, tiempo_cuadrado), start=1):
    print(f"{i:<8}{h:<12.3f}{t:<12.3f}{t2:<12.6f}")
print("="*42)

# Mostrar resultados del ajuste
print("\n")
print("="*60)
print("Resultados del ajuste lineal - Altura vs. Tiempo²")
print("="*60)
print(f" Pendiente (g/2): {m:.3f} ± {error_m:.3f} m/s²")
print(f" Ordenada al origen: {b:.3f} ± {error_b:.3f} m")
print(f" Coeficiente de determinación (R²): {r2:.3f}")
# Calcular y mostrar la aceleración de la gravedad
aceleracion = 2 * m
error_aceleracion = 2 * error_m
aex = aceleracion + error_aceleracion
ade = aceleracion - error_aceleracion
print(f" Aceleración de la gravedad: g = {aceleracion:.3f} ± {error_aceleracion:.3f} m/s²")
print(f" Aceleración por exceso: g = {aex:.3f} m/s²")
print(f" Aceleración por defecto: g = {ade:.3f} m/s²")
print("="*60)