# -*- coding: utf-8 -*-
"""ejercicio_2_primera_entrega.ipynb

Automatically generated by Colab.

Original file is located at
    https://colab.research.google.com/drive/1BcZDWaydxID0J39XNnQbUhJOVutdL7Di
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import newton


def f(Q):
    return 12 * np.sin(Q / 5)

def f_prime(Q):
    return 12 / 5 * np.cos(Q / 5)

def ff(Q):
    return f(f(Q))

def ff_prime(Q):
    return f_prime(f(Q)) * f_prime(Q)

# Encontrar puntos de equilibrio
def find_equilibrium_points():
    Q_values = np.linspace(0, 5 * np.pi, 1000)
    tolerance = 1e-5
    equilibrium_points = []
    for Q0 in Q_values:
        try:
            root = newton(func=lambda Q: f(Q) - Q, fprime=lambda Q: f_prime(Q) - 1, x0=Q0, tol=tolerance)
            if all(abs(root - eq) > tolerance for eq in equilibrium_points):
                equilibrium_points.append(root)
        except (RuntimeError, OverflowError):
            continue
    return equilibrium_points

# Clasificar puntos de equilibrio
def classify_equilibrium_points(points):
    classification = []
    for eq in points:
        if 0 <= eq <= 5 * np.pi:
            derivative_at_eq = f_prime(eq)
            stability = "estable" if abs(derivative_at_eq) < 1 else "inestable"
            classification.append((eq, derivative_at_eq, stability))
    return classification

# Encontrar y clasificar órbitas de periodo 2
def find_period_2_orbits():
    Q_values = np.linspace(0, 5 * np.pi, 1000)
    tolerance = 1e-5
    period_2_points = []
    for Q0 in Q_values:
        try:
            root = newton(func=lambda Q: ff(Q) - Q, fprime=ff_prime, x0=Q0, tol=tolerance)
            if all(abs(root - eq) > tolerance for eq in period_2_points):
                period_2_points.append(root)
        except (RuntimeError, OverflowError):
            continue
    period_2_classification = []
    for p2 in period_2_points:
        if 0 <= p2 <= 5 * np.pi:
            derivative_at_p2 = ff_prime(p2)
            stability = "atractora" if abs(derivative_at_p2) < 1 else "repulsora"
            period_2_classification.append((p2, derivative_at_p2, stability))
    return period_2_classification

# Diagramas de telaraña
def cobweb_plot(Q0, f, ax, iterations=50):
    Q_values = [Q0]
    for _ in range(iterations):
        Q_next = f(Q_values[-1])
        Q_values.append(Q_next)
    line = np.linspace(0, 5 * np.pi, 400)
    ax.plot(line, f(line), 'b', lw=2, label='f(Q)')
    ax.plot(line, line, 'k', lw=1)
    for i in range(iterations):
        Q_n, Q_np1 = Q_values[i], Q_values[i+1]
        ax.plot([Q_n, Q_n], [Q_n, Q_np1], 'r', lw=1)
        ax.plot([Q_n, Q_np1], [Q_np1, Q_np1], 'r', lw=1)
    ax.set_xlim([0, 5 * np.pi])
    ax.set_ylim([0, 15])
    ax.set_title(f"Cobweb Plot for Q0 = {Q0}")
    ax.set_xlabel('Q_n')
    ax.set_ylabel('Q_{n+1}')


equilibriums = find_equilibrium_points()
equilibrium_classes = classify_equilibrium_points(equilibriums)
print("Puntos de equilibrio y su estabilidad:")
for eq, _, stability in equilibrium_classes:
    print(f"Q = {eq:.2f} es {stability}")

period_2_orbits = find_period_2_orbits()
print("\nÓrbitas de periodo 2 y su estabilidad:")
for p2, _, stability in period_2_orbits:
    print(f"Órbita en Q = {p2:.2f} es {stability}")

fig, axes = plt.subplots(1, 3, figsize=(18, 6))
initial_conditions = [5, 10, 15]
for ax, Q0 in zip(axes, initial_conditions):
    cobweb_plot(Q0, f, ax)
plt.tight_layout()
plt.show()