1. ¿Por qué modelar esto? 🌠

Porque con \(N>2\) no hay solución analítica — aparece el caos determinista. La compu se vuelve telescopio: podemos contar historias orbitales en cámara lenta y aprender de estabilidad, energía y transferencia de momento.

Mini-glosario express

• Cuerpo/partícula: agente con masa \(m\).
• \(\vec r_{ij}\): vector desde \(i\) a \(j\); \(r_{ij}=||\vec r_{ij}||\).
• \(G\): constante gravitacional (usamos u.a.: \(G=1\)).
• \(\vec v, \vec a\): velocidad y aceleración.
• COM: centro de masa. Quitamos su movimiento para que la “escena” no se nos vaya del encuadre.
• \(\varepsilon\) (softening): acolchado numérico: \(\sqrt{r^2+\varepsilon^2}\) evita singularidades en choques muy cercanos.
• \(\Delta t\): paso de tiempo; muy grande ⇒ inestabilidad, muy chico ⇒ más lento pero más fiel.

Dato musical random: Walking on the Moon (The Police, 1979) es perfecto de fondo; en nuestro universo nadie camina… ¡todos orbitan! 🎶🌙

---

2. Física exprés (pero bien puesta) 🧠🔭

• Segunda ley de Newton para cada cuerpo \(i\):  \[  m_i\,\vec a_i \;=\; \sum_{j\neq i} \vec F_{ij}.  \]

• Integración temporal (Velocity-Verlet): con aceleración \(\vec a(t)\),  \[  \vec x_{t+\Delta t} = \vec x_t + \vec v_t\,\Delta t + \tfrac12\vec a_t\,\Delta t^2,\qquad  \vec v_{t+\Delta t} = \vec v_t + \tfrac12\left(\vec a_t + \vec a_{t+\Delta t}\right)\Delta t.  \]  Es casi conservativo de energía (mucho más que Euler).

• Energía total (con suavizado \(\varepsilon\) para evitar singularidades numéricas):  \[  E = \underbrace{\tfrac12\sum_i m_i |\vec v_i|^2}_{\text{Cinética}}      \;-\; \underbrace{\sum_{i<j}\frac{G\,m_i m_j}{\sqrt{r_{ij}^2+\varepsilon^2}}}_{\text{Potencial}}.  \]

• Bound vs. escape (intuición 2-cuerpos): la energía específica  \[  \varepsilon_{\text{orb}}=\frac{v^2}{2}-\frac{G M}{r}  \]  < 0 → elipse (ligado), = 0 → parábola, > 0 → hipérbola (escape). En N-cuerpos, los slingshots redistribuyen energía y alguno sale volando.

---

3. Algoritmo estilo “bullet-time” 👾

1. Inicializa masas, posiciones y velocidades.
2. Quita el momento lineal total para que el centro de masa quede quietecito.
3. Calcula fuerzas (con softening \(\varepsilon\)).
4. Integra con Velocity-Verlet.
5. Guarda historia de posiciones, energía y haz snapshots.
6. (Opcional) Anima y exporta GIF/MP4.

Dato cultura-latam: en 1846, Le Verrier predijo Neptuno sin verlo, solo con lápiz, papel y dinámica gravitacional. OG del data science orbital. 📜🔭

---

4. Código listo para despegar 🚀🐍

Corre tal cual. Hace: simulación, panel de snapshots (4 escenas) y curva de energía.

• Guarda medios si SAVE_MEDIA=True (rutas al final del script).
• Activa el GIF con MAKE_GIF=True (requiere Pillow).

import numpy as np
import matplotlib.pyplot as plt
 
# =================== Parámetros ===================
G         = 1.0
DT        = 0.002          # paso temporal
STEPS     = 1500           # total de pasos
N         = 5              # número de cuerpos
EPS       = 1e-2           # suavizado gravitacional
SEED      = 20910          # para reproducibilidad
PRESET    = "random"       # "random" o "solarish"
MAKE_GIF  = False          # True para exportar GIF (requiere Pillow)
SAVE_MEDIA= False          # True para guardar PNGs de snapshots/energía
TRAIL     = 250            # largo de estela en la animación
 
# rutas de salida (ajusta a tu servidor)
PATH_GIF       = "nbody.gif"
PATH_SNAPSHOTS = "snapshots.png"
PATH_ENERGY    = "energia.png"
 
rng = np.random.default_rng(SEED)
 
# =================== Estado inicial ===================
def init_state(preset="random"):
    if preset == "solarish":
        # Estrella gorda al centro + 4 planetitas casi circulares con ligera perturbación
        m = np.array([12.0, 2.5, 2.0, 1.5, 1.2])
        pos = np.zeros((N,2)); vel = np.zeros((N,2))
        radii = np.array([1.2, 1.8, 2.7, 3.5])   # radios
        for k, r in enumerate(radii, start=1):
            ang = rng.uniform(0, 2*np.pi)
            pos[k] = r * np.array([np.cos(ang), np.sin(ang)])
            v_circ = np.sqrt(G*m[0]/r)
            tang = np.array([-np.sin(ang), np.cos(ang)])
            vel[k] = 1.0*rng.normal(1.0, 0.05) * v_circ * tang
        # estrella al centro, con pequeña velocidad para cerrar momento total
        vel[0] = rng.normal(0, 0.02, 2)
    else:
        # Random pinball para ver slingshots y eyecciones
        m   = rng.uniform(1.0, 5.0, N)
        m[0] += 5.0   # que haya una "estrella" razonable
        pos = rng.uniform(-4, 4, (N,2))
        vel = rng.uniform(-0.4, 0.4, (N,2))
    # quitar momentum total (usar el vector, no el escalar)
    P = (m[:,None]*vel).sum(axis=0)              # (2,)
    vel -= P / m.sum()
    return m, pos, vel
 
m, pos, vel = init_state(PRESET)
 
# =================== Fuerzas y energía ===================
def forces(pos, m, G=G, eps=EPS):
    N = len(m)
    F = np.zeros_like(pos)
    for i in range(N):
        r = pos[i] - pos                         # (N,2)
        dist2 = (r*r).sum(axis=1) + eps**2       # (N,)
        inv32 = 1.0 / (dist2**1.5)               # (N,)
        mask = np.ones(N, bool); mask[i] = False
        contrib = -G * m[i] * m[mask][:,None] * r[mask] * inv32[mask,None]
        F[i] = contrib.sum(axis=0)               # (2,)
    return F
 
def total_energy(pos, vel, m, G=G, eps=EPS):
    KE = 0.5 * np.sum(m * np.sum(vel**2, axis=1))
    PE = 0.0
    for i in range(len(m)-1):
        r = pos[i+1:] - pos[i]
        dist = np.sqrt((r*r).sum(axis=1) + eps**2)
        PE += -G * m[i] * np.sum(m[i+1:] / dist)
    return KE + PE
 
# =================== Simulación (Velocity-Verlet) ===================
pos_hist = np.zeros((STEPS, N, 2))
E_hist   = np.zeros(STEPS)
 
F   = forces(pos, m)
acc = F / m[:,None]
 
for t in range(STEPS):
    pos += vel*DT + 0.5*acc*(DT**2)
    F_new   = forces(pos, m)
    acc_new = F_new / m[:,None]
    vel += 0.5*(acc + acc_new)*DT
    acc  = acc_new
    pos_hist[t] = pos
    E_hist[t]   = total_energy(pos, vel, m)
 
# =================== Snapshots (4 escenas) ===================
def plot_snapshots(
    pos_hist,
    m,
    steps_list,
    title_prefix="Evolución N-cuerpos",
    save_path=None,
    shared=False,          # <= False: límites por-subplot (zoom por escena). True: límites globales
    pad_frac=0.08          # margen alrededor de los datos en cada subplot
): 
    # Figura: compartir ejes solo si shared=True
    if shared:
        fig, axes = plt.subplots(2, 2, figsize=(10,7), sharex=True, sharey=True)
        # límites globales (en todo el run) si se comparten
        Rg = pos_hist.reshape(-1, 2)
        padg = pad_frac * (Rg.max(axis=0) - Rg.min(axis=0) + 1e-9)
        xlim_g = (Rg[:,0].min()-padg[0], Rg[:,0].max()+padg[0])
        ylim_g = (Rg[:,1].min()-padg[1], Rg[:,1].max()+padg[1])
    else:
        fig, axes = plt.subplots(2, 2, figsize=(10,7))
 
    axes = axes.ravel()
    for ax, s in zip(axes, steps_list):
        s = min(s, len(pos_hist)-1)
 
        # trayectorias hasta el paso s
        for i in range(len(m)):
            ax.plot(pos_hist[:s, i, 0], pos_hist[:s, i, 1], lw=1.2)
 
        # posiciones en el paso s
        ax.scatter(
            pos_hist[s, :, 0], pos_hist[s, :, 1],
            s=30*(m/np.min(m)), edgecolor='k', linewidth=0.6, alpha=0.9
        )
        ax.set_title(f"t = {s} pasos")
        ax.set_aspect('equal', 'box')
 
        if shared:
            ax.set_xlim(*xlim_g); ax.set_ylim(*ylim_g)
        else:
            # límites SOLO con datos hasta s (más zoom y detalle)
            Rs = pos_hist[:s+1].reshape(-1, 2)
            pads = pad_frac * (Rs.max(axis=0) - Rs.min(axis=0) + 1e-9)
            ax.set_xlim(Rs[:,0].min()-pads[0], Rs[:,0].max()+pads[0])
            ax.set_ylim(Rs[:,1].min()-pads[1], Rs[:,1].max()+pads[1])
 
    fig.suptitle(f"{title_prefix}: snapshots", fontsize=14)
    for ax in axes[-2:]: ax.set_xlabel("x")
    for ax in axes[::2]: ax.set_ylabel("y")
    plt.tight_layout()
    if save_path:
        plt.savefig(save_path, dpi=200, bbox_inches='tight')
    plt.show()
 
steps_list = [int(0.1*STEPS), int(0.3*STEPS), int(0.6*STEPS), int(0.9*STEPS)]
plot_snapshots(pos_hist, m, steps_list,
               title_prefix="Órbitas y trayectorias — Verlet",
               save_path=(PATH_SNAPSHOTS if SAVE_MEDIA else None), shared=False)
 
# =================== Chequeo de energía ===================
plt.figure(figsize=(7,3.2))
E0 = E_hist[0]
plt.plot(np.arange(len(E_hist))*DT, (E_hist - E0)/abs(E0))
plt.xlabel("Tiempo"); plt.ylabel("Deriva de energía relativa")
plt.title("Estabilidad numérica (Velocity-Verlet)")
plt.tight_layout()
if SAVE_MEDIA:
    plt.savefig(PATH_ENERGY, dpi=200, bbox_inches='tight')
plt.show()
 
# =================== Animación (opcional) ===================
if MAKE_GIF:
    from matplotlib.animation import FuncAnimation, PillowWriter
    fig, ax = plt.subplots(figsize=(6,4))
    ax.set_aspect('equal', 'box')
    ax.set_xlim(np.min(pos_hist[:,:,0])-1, np.max(pos_hist[:,:,0])+1)
    ax.set_ylim(np.min(pos_hist[:,:,1])-1, np.max(pos_hist[:,:,1])+1)
    lines = [ax.plot([], [], lw=1.2)[0] for _ in range(N)]
    pts = ax.scatter([], [], s=30*(m/np.min(m)), edgecolor='k', linewidth=0.6)
    def init():
        for ln in lines: ln.set_data([], [])
        pts.set_offsets(np.zeros((N,2))); return lines + [pts]
    def update(frame):
        lo = max(0, frame-TRAIL)
        for i, ln in enumerate(lines):
            ln.set_data(pos_hist[lo:frame, i, 0], pos_hist[lo:frame, i, 1])
        pts.set_offsets(pos_hist[frame-1]); return lines + [pts]
    ani = FuncAnimation(fig, update, frames=len(pos_hist), init_func=init, interval=20, blit=True)
    ani.save(PATH_GIF, writer=PillowWriter(fps=40))   # guarda GIF
    plt.close(fig)
    print(f"GIF guardado como {PATH_GIF} ✅")

---

5. Resultados y discusión — Caso semilla 20910 📊🛰️

El siguiente panel resume la “película” de este caso (snapshots con zoom por escena). Reproduce este panel: usa SEED=20910, STEPS=1500, DT=0.002, EPS=1e-2, PRESET="random" y shared=False en los snapshots.

Izq→Der, Arriba→Abajo: acercamiento → close pass → tríada → dos eyecciones; colores coherentes por cuerpo.

Lectura guiada (qué está pasando en cada cuadro):

• t = 150 pasos — “La reunión incómoda”
  Cinco cuerpos, una masa dominante (círculo más grande) cerca del cuadrante inferior derecho. Geometría abierta: varios pares encaminados a un primer flyby.

• t = 450 — Primer close pass
  Dos ligeros pasan cerca de la masa grande → curvaturas marcadas y primer intercambio de energía. Empieza el pasito slingshot 💃.

• t = 900 — Trío temporal
  Se forma un sistema jerárquico transitorio: la masa grande + un compañero quedan ligados en torno a \((x,y)\approx(0,-3)\) con órbitas excéntricas; los otros dos hacen pases a alta velocidad. Aquí ya se “decide” quién se va.

• t = 1350 — Resolución
  Dos cuerpos son eyectados en trayectorias abiertas (hipérbolas aparentes dentro del encuadre), uno hacia la derecha y otro hacia el cuadrante superior izquierdo.
  La masa grande mantiene un binario excéntrico con un compañero cercano. Misión: accomplished.

💡 Criterio de escape (por si quieres etiquetar quién huye): para un cuerpo de masa \(m\) respecto a otra de masa \(M\) a distancia \(r\) y velocidad relativa \(v\), \[ \varepsilon_{\text{orb}} \;=\; \frac{v^2}{2}\;-\;\frac{G M}{r}. \] Si \(\varepsilon_{\text{orb}}>0\) ⇒ trayectoria abierta (escape). Si \(\varepsilon_{\text{orb}}<0\) ⇒ ligado (elipse).

Chequeo de energía (este run):

Observa el salto en \(t\approx 1.6\) (un encuentro muy cercano). Dos lecturas útiles:

1. La dinámica cualitativa sigue razonable (Verlet ayuda), pero
2. La métrica \((E-E_0)/|E_0|\) se dispara cuando \(E_0\) es pequeño en magnitud o hay pases muy cercanos. Si quieres aplanar esa curva en este mismo caso, prueba cualquiera de estos tunings soft:
      • bajar DT a 0.001
      • subir EPS a 0.02
      • activar fusión si \(r<r_\text{merge}\) para evitar “casi-singularidades”.

---

6. Laboratorio express con esta misma semilla 🔬🎮

Objetivo: explorar con este mismo estado inicial cómo cambian los desenlaces.

1. ¿Eyectado o capturado?
   Calcula \(\varepsilon_{\text{orb}}\) de cada cuerpo respecto a la masa grande cuadro a cuadro. Marca en la gráfica cuando cruce de \(\varepsilon_{\text{orb}}<0\) a \(>0\) → momento de escape.

2. Afinando estabilidad (mismo drama, menos drama-queen):
   • DT=0.001, EPS=0.02: la historia macroscópica no cambia (dos escapes, un binario), pero la deriva relativa cae drásticamente.
   • DT=0.003 (no recomendado): verás más “serrucho” y posibles artefactos.

3. Comparativa “Solar-ish” vs. “Random pinball”:
   Cambia PRESET="solarish" (en el mismo script) y repite los snapshots: pasan de eyecciones frecuentes a órbitas casi circulares y resonancias \(2{:}1\), \(3{:}2\) si añades un 5–10 % de perturbación en \(v\).

4. Mini-reto visual:
   Vuelve a generar los snapshots con shared=True (misma escala) y compáralos con shared=False (zoom por escena). Discute pros/cons: detalle local vs. narrativa global.

Troubleshooting (si algo se “rompe”)

• Saltos de energía: baja DT → 0.001 o sube EPS → 0.02.
• Todo sale disparado recto: o hubo slingshot bueno (válido) o DT es muy grande.
• Nada orbita: sube la masa central o inicializa vel. tangenciales (\(v\approx\sqrt{GM/r}\)).
• Escala rara en snapshots: usa shared=False (zoom por escena).

Micro-pista musical: prueba poner de fondo Kraftwerk; el pulso electrónico combina increíble con el paso fijo de \( \Delta t \). 🕺

---

7. Cierre (versión caso 20910) 📎

Con este seed vimos el arco canónico del pinball gravitacional: acercamiento → slingshots → dos escapes + un binario remanente. Lo lindo es que no “forzamos” nada: la estadística de encuentros hace la coreografía. La próxima vez que escuches sobre una asistencia gravitatoria, recuerda: es literalmente intercambiar \(\frac{v^2}{2}\) por \(\frac{GM}{r}\) con estilo. 🎷🪐