import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# Parámetros del sistema
alpha = -0.2
beta = 1.0
gamma = 0.9
omega = 0.6  # Cambia a 0.3, 0.7, 0.8, etc.

# Sistema de ecuaciones
def forced_pendulum(t, z):
    x, y = z
    dxdt = y
    dydt = alpha * y + beta * np.sin(x) + gamma * np.cos(omega * t)
    return [dxdt, dydt]

# Tiempo de integración
n_cycles = 50
points_per_cycle = 100
T = 2 * np.pi / omega
t_max = n_cycles * T
t_eval = np.linspace(0, t_max, n_cycles * points_per_cycle)

# Condiciones iniciales
z0 = [0.1, 0.0]

# Integrar
sol = solve_ivp(forced_pendulum, [0, t_max], z0, t_eval=t_eval, rtol=1e-10, atol=1e-10)

# Extraer soluciones
x = sol.y[0]
y = sol.y[1]
t = sol.t

# Diagrama de fase
plt.figure(figsize=(10, 4))
plt.plot(x, y, lw=0.5)
plt.title(f"Diagrama de fase para ω = {omega}")
plt.xlabel("x")
plt.ylabel("y")
plt.grid(True)
plt.tight_layout()
plt.show()

# Poincaré: cada T segundos
poincare_times = np.arange(0, t_max, T)
poincare_indices = [np.argmin(np.abs(t - pt)) for pt in poincare_times]
x_poincare = x[poincare_indices]
y_poincare = y[poincare_indices]

# Diagrama de Poincaré
plt.figure(figsize=(6, 6))
plt.plot(x_poincare, y_poincare, 'o', markersize=1)
plt.title(f"Diagrama de Poincaré para ω = {omega}")
plt.xlabel("x")
plt.ylabel("y")
plt.grid(True)
plt.tight_layout()
plt.show()