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

# ---------------------------
# Parámetros del sistema
# ---------------------------
m = 1.0
l = 2.0
R = 1.0
g = 9.81
w = 1.74  # frecuencia angular (omega)
T = 2 * np.pi / w
tf = 1000 * T

# ---------------------------
# Condiciones iniciales (rejilla tipo Mathematica)
# ---------------------------
phi_degrees = np.arange(0, 181, 10)
v_values    = np.arange(-5, 6, 1)
CI = [(np.deg2rad(phi0), v0) for phi0 in phi_degrees for v0 in v_values]

# ---------------------------
# Wrap a [-pi, pi)
# ---------------------------
def wrap(theta):
    return (theta + np.pi) % (2.0 * np.pi) - np.pi

# ---------------------------
# EDO: x' = v
#      v' = (R*w^2*cos(x - w t) - g*sin(x))/l
# ---------------------------
def f(t, y):
    x, v = y
    dxdt = v
    dvdt = (R * w**2 * np.cos(x - w * t) - g * np.sin(x)) / l
    return [dxdt, dvdt]

# ---------------------------
# Tiempos estroboscópicos: t_n = n T, n = 0..1000
# ---------------------------
nmax = int(np.floor(tf / T))
t_samples = (np.arange(0, nmax + 1) * T)

# ---------------------------
# Integración para cada CI y construcción de puntos (wrap[x(t_n)], xdot(t_n))
# ---------------------------
PE = []  # lista de arreglos (N,2), uno por cada CI

# Control numérico: max_step ayuda a no "saltar" dinámica
max_step = T / 50.0

for (x0, v0) in CI:
    sol = solve_ivp(
        f,
        t_span=(0.0, tf),
        y0=[x0, v0],
        t_eval=t_samples,
        method="DOP853",      # buena opción para integración larga
        rtol=1e-9,
        atol=1e-12,
        max_step=max_step
    )
    x = sol.y[0]
    v = sol.y[1]
    pts = np.column_stack((wrap(x), v))
    PE.append(pts)

# ---------------------------
# Gráfico estroboscópico
# ---------------------------
fig, ax = plt.subplots(figsize=(8, 4.8))
for pts in PE:
    ax.scatter(pts[:, 0], pts[:, 1], s=1)  # puntos pequeños
ax.text(
    0.98, 0.98, rf"$\omega = {w:.2f}$",
    transform=ax.transAxes,
    ha="right", va="top",
    fontsize=16,
    bbox=dict(facecolor="white", alpha=0.7, edgecolor="none", pad=3)
)
ax.set_xlabel(r"$\phi$", fontsize=16)
ax.set_ylabel(r"$\dot{\phi}$", fontsize=16)
ax.grid(True, alpha=0.5)
ax.tick_params(axis="both", which="major", labelsize=16)
fig.tight_layout()
plt.show()