### N-Body Simulation Setup and Initialization

Source: https://github.com/seap-udea/pymcel/blob/main/ejemplos/cuadernos-libro/mcel_zuluaga-10-ncuerpos_solucion_numerica.ipynb

Initializes an N-body simulation with specific masses and positions, then converts this system into the format required by the numerical integration function (`edm_ncuerpos`). It sets up the time points for integration and prepares the plot axes, but does not perform the integration or plotting itself.

```python
# Definición de las condiciones iniciales 
sistema_ejemplo=[
    dict(m=1.0,r=[0.0,0.0,0.0],v=[0.0,0.0,0.0]),
    dict(m=0.2,r=[1.0,0.0,0.0],v=[0.0,1.0,0.0]),
    dict(m=0.0,r=[-1.0,0.0,0.0],v=[0.0,-1.0,0.0]),
]
from pymcel import sistema_a_Y
N,mus,Yo=sistema_a_Y(sistema_ejemplo)

#Tiempo de integración
import numpy as np
Nt=50
ts=np.linspace(0.0,10.0,Nt,endpoint=True)

# Solución al sistema de ecuaciones diferenciales
from scipy.integrate import odeint
solucion_ejemplo=odeint(edm_ncuerpos,Yo,ts,args=(N,mus))
rs,vs=solucion_a_estado(solucion_ejemplo,N,Nt)

# Componente gráfica del algoritmo
import matplotlib.pyplot as plt
fig=plt.figure()
ax=fig.gca()

puntos=[]
lineas=[]
for i in range(N):
    punto,=ax.plot(rs[i,0,0],rs[i,0,1],marker='o',lw=0);
    linea,=ax.plot(rs[i,0,0],rs[i,0,1],marker=None,lw=1,color='k');
    puntos+=[punto]
    lineas+=[linea]

titulo=ax.set_title(f"t = {ts[0]}")
vmax=abs(rs).max()
ax.set_ylim((-vmax,vmax))
ax.set_xlim((-vmax,vmax))
ax.grid()

```

--------------------------------

### Initialize System for N-body Simulation (Python)

Source: https://github.com/seap-udea/pymcel/blob/main/docs/examples/mcel_zuluaga-19-trescuerpos_aplicaciones.ipynb

Defines the initial state of celestial bodies for an N-body simulation. It includes mass, position, and velocity for each body. This setup is used as input for the `ncuerpos_solucion` function.

```python
sistema=[
    #Partícula de prueba
    dict(
        m=1e-5,
        r=[6.0,0,-0.1],
        v=[0,2.0,0.1]
    ),
    # Particula 1
    dict(
        m=1000.0,
        r=[-0.005,0,0],
        v=[0,-7.075,0]
    ),
    # Particula 2
    dict(
        m=1.0,
        r=[4.995,0,0],
        v=[0,+7.075,0]
    ),
]

```

--------------------------------

### N-Body Simulation Initialization and Plotting Setup (Python)

Source: https://github.com/seap-udea/pymcel/blob/main/docs/examples/mcel_zuluaga-10-ncuerpos_solucion_numerica.ipynb

Initializes an N-body simulation with specific masses and positions, then sets up the integration time. It performs the numerical integration and prepares for plotting by creating plot elements for each body's trajectory, including points and lines, and setting up the initial plot title and axes limits.

```python
# Definición de las condiciones iniciales 
sistema_ejemplo=[
    dict(m=1.0,r=[0.0,0.0,0.0],v=[0.0,0.0,0.0]),
    dict(m=0.2,r=[1.0,0.0,0.0],v=[0.0,1.0,0.0]),
    dict(m=0.0,r=[-1.0,0.0,0.0],v=[0.0,-1.0,0.0]),
]
from pymcel import sistema_a_Y
N,mus,Yo=sistema_a_Y(sistema_ejemplo)

#Tiempo de integración
import numpy as np
Nt=50
ts=np.linspace(0.0,10.0,Nt,endpoint=True)

# Solución al sistema de ecuaciones diferenciales
from scipy.integrate import odeint
solucion_ejemplo=odeint(edm_ncuerpos,Yo,ts,args=(N,mus))
rs,vs=solucion_a_estado(solucion_ejemplo,N,Nt)

# Componente gráfica del algoritmo
import matplotlib.pyplot as plt
fig=plt.figure()
ax=fig.gca()

puntos=[]
lineas=[]
for i in range(N):
    punto,=ax.plot(rs[i,0,0],rs[i,0,1],marker='o',lw=0);
    linea,=ax.plot(rs[i,0,0],rs[i,0,1],marker=None,lw=1,color='k');
    puntos+=[punto]
    lineas+=[linea]

titulo=ax.set_title(f"t = {ts[0]}")
vmax=abs(rs).max()
ax.set_ylim((-vmax,vmax))
ax.set_xlim((-vmax,vmax))
ax.grid()
```

--------------------------------

### Initialize Particle System and Convert to Y0s

Source: https://github.com/seap-udea/pymcel/blob/main/docs/examples/mcel_zuluaga-10-ncuerpos_solucion_numerica.ipynb

Initializes a system of N particles with random masses, positions, and velocities. It then converts this system into the Y0s format required by the `pymcel` library using the `sistema_a_Y` function. This setup is crucial for simulating the system's dynamics.

```python
#Número de partículas
N=5

#Generación de las condiciones para cada partícula
from numpy.random import uniform,seed
seed(7)

#Condiciones iniciales
sistema=[]
for i in range(N):
    particula=dict(
        m=uniform(0.0,1.0),
        r=uniform(-1.0,1.0,size=3),
        v=uniform(-1.0,1.0,size=3)
    )
    sistema+=[particula]

from pymcel import sistema_a_Y
N,mus,Y0s=sistema_a_Y(sistema)

#Tiempos
from numpy import linspace
Nt=100
ts=linspace(0.0,10.0,Nt)

print(f"N = {N}")
print(f"mus = {mus}")
print(f"Y0s = {Y0s}")
```

--------------------------------

### Animate Conical Pendulum Motion in Coordinate and Phase Spaces (Python)

Source: https://github.com/seap-udea/pymcel/blob/main/ejemplos/cuadernos-libro/mcel_zuluaga-22-hamiltoniano.ipynb

This script simulates the conical pendulum's motion and generates an animated visualization. It reuses the simulation setup from previous examples and creates an animation function that updates the positions of points and lines in both the coordinate and phase space plots. The animation is then embedded as HTML.

```python
#Propiedades
g=9.81
L=1.0
m=1.0

#Condiciones iniciales
teta_0=15*grados
fi_0=0.0*grados
pteta_0=0.0
alfa_fi=0.8

#Tiempos de integración
T=2*pi*sqrt(L/g)
Nt=200
ts=linspace(0,3*T,Nt)

#Solución numérica
from scipy.integrate import odeint
solucion=odeint(edm_penduloconico_hamiltoniano,
                [teta_0,fi_0,pteta_0],ts,
                args=(alfa_fi,L,m))
tetas=solucion[:,0]
fis=solucion[:,1]
ptetas=solucion[:,2]

#Coordenadas cartesianas
xs=L*sin(tetas)*cos(fis)
ys=L*sin(tetas)*sin(fis)
zs=-L*cos(tetas)

#Preparación del gráfico
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig=plt.figure(figsize=(8,4))
ax_coord=fig.add_subplot(121,projection='3d')
ax_fase=fig.add_subplot(122,projection='3d')

#Gráfico en el espacio de coordenadas
ax_coord.plot(xs,ys,zs,'k-',alpha=0.3)
punto_coord,=ax_coord.plot([],[],[],'ko',ms=5)
cuerda,=ax_coord.plot([],[],[],'k-')

#Gráfico en el espacio de fases
ax_fase.plot(tetas,fis,ptetas,'b-',alpha=0.3)
punto_fase,=ax_fase.plot([],[],[],'bo',ms=5)

#Decoración
ax_coord.set_xlabel("$x$")
ax_coord.set_ylabel("$y$")
ax_coord.set_zlabel("$z$")
ax_coord.set_title("Espacio coordenado")
from pymcel.plot import fija_ejes3d_proporcionales
fija_ejes3d_proporcionales(ax_coord)
ax_fase.set_xlabel("$\theta$")
ax_fase.set_ylabel("$phi$")
ax_fase.set_zlabel("$p_\theta$")
ax_fase.set_title("Espacio de fase")
fig.tight_layout()

def animacion(it):
    punto_coord.set_data_3d([xs[it]],[ys[it]],[zs[it]])
    cuerda.set_data_3d([0,xs[it]],[0,ys[it]],[0,zs[it]])
    punto_fase.set_data_3d([tetas[it]],[fis[it]],[ptetas[it]])
    return punto_coord,cuerda,punto_fase
                          
from matplotlib import animation
anim=animation.FuncAnimation(fig,animacion,frames=Nt,interval=50,blit=True,repeat=False);
plt.close('all')

from matplotlib import rcParams
rcParams['animation.embed_limit']=2**128
from IPython.display import HTML
HTML(anim.to_jshtml())
```

--------------------------------

### Install geopandas Package

Source: https://github.com/seap-udea/pymcel/blob/main/ejemplos/atico.ipynb

This command installs or upgrades the 'geopandas' Python package. Geopandas is used for working with geospatial data in Python and may be used in conjunction with orbital data for mapping or analysis.

```bash
!pip install -Uq geopandas
```

--------------------------------

### Instalar pymcel en Colab

Source: https://github.com/seap-udea/pymcel/blob/main/docs/examples/mcel_zuluaga-01-introduccion.ipynb

Comando para instalar la librería pymcel en un entorno de Google Colab. Asegura que la versión más reciente esté disponible para su uso.

```python
!pip install -Uq pymcel
```

--------------------------------

### Install sgp4 Package

Source: https://github.com/seap-udea/pymcel/blob/main/ejemplos/atico.ipynb

This command installs or upgrades the 'sgp4' Python package. The 'sgp4' library is essential for propagating satellite orbits using the SGP4/SDP4 models, which are commonly used with TLE data.

```bash
!pip install -Uq sgp4
```

--------------------------------

### Setting up a Simulated System

Source: https://github.com/seap-udea/pymcel/blob/main/docs/examples/mcel_zuluaga-15-trescuerpos_constante_jacobi.ipynb

This code prepares a simulated system by defining initial properties and conditions. It sets the alfa parameter, initial position (ro), initial velocity (vo), and defines the time points for integration.

```python
#Propiedades del sistema y condiciones iniciales
alfa=0.3
ro=[1.0,0.0,0.0]
vo=[0.0,0.45,0.0]
#Tiempos de integración
from numpy import linspace
Nt=1000
ts=linspace(0,10,Nt)
#Resuelve numéricamente la ecuación de movimiento
from pymcel import crtbp_solucion
solucion=crtbp_solucion(alfa,ro,vo,ts)
#Extrae las posiciones y velocidades en el sistema rotante
rs=solucion[0]
vs=solucion[1]
```

--------------------------------

### Initialize System and Solve N-body Equations (Python)

Source: https://github.com/seap-udea/pymcel/blob/main/ejemplos/cuadernos-libro/mcel_zuluaga-19-trescuerpos_aplicaciones.ipynb

Sets up a system of celestial bodies with their masses, initial positions, and velocities. It then solves the N-body equations of motion over a specified time interval and plots the trajectories relative to the center of mass. Dependencies include numpy and pymcel.

```python
sistema=[
    #Partícula de prueba
    dict(
        m=1e-5,
        r=[6.0,0,-0.1],
        v=[0,2.0,0.1]
    ),
    # Particula 1
    dict(
        m=1000.0,
        r=[-0.005,0,0],
        v=[0,-7.075,0]
    ),
    # Particula 2
    dict(
        m=1.0,
        r=[4.995,0,0],
        v=[0,+7.075,0]
    ),
]

T=20.0
Nt=1000
from numpy import linspace
ts=linspace(0,T,Nt)
from pymcel import ncuerpos_solucion
rs,vs,rps,vps,constantes=ncuerpos_solucion(sistema,ts)
from pymcel import plot_ncuerpos_3d
fig=plot_ncuerpos_3d(rps,vps);
ax=fig.gca()
```

--------------------------------

### Plot 3D Periodic Orbit Example Near L1 (Python)

Source: https://github.com/seap-udea/pymcel/blob/main/ejemplos/cuadernos-libro/mcel_zuluaga-19-trescuerpos_aplicaciones.ipynb

Example demonstrating a periodic trajectory related to the L1 point in 3D within the CRTBP. It uses `scipy.optimize.bisect` to find the L1 point and `pymcel.orbitas_crtbp3d` for visualization. The output is a 3D matplotlib figure.

```python
#Propiedades del sistema
alfa=0.0121505856
xL1=bisect(funcion_puntos_colineales,0,1-2*alfa,args=(alfa,)) # Assuming funcion_puntos_colineales is defined elsewhere
#Condiciones iniciales
ro=[0.8329,0,0.1304]
vo=[0,0.2437,0]
fig=orbitas_crtbp3d(alfa,ro,vo,
                T=3,Nt=100,
                xlim=(-0.1,1.1),ylim=(-0.5,0.5),zlim=(-0.5,0.5),
                xL=xL1,yL=0,zL=0
               )
```

--------------------------------

### Autoreload Extension Setup

Source: https://github.com/seap-udea/pymcel/blob/main/ejemplos/cuadernos-libro/pymcel/data/mcel_zuluaga-plantilla.ipynb

Configures the IPython environment to automatically reload modules before executing code. This is useful during development to see changes without restarting the kernel.

```python
%load_ext autoreload 
%autoreload 2
```

--------------------------------

### Instalar PymCel desde PyPI

Source: https://github.com/seap-udea/pymcel/blob/main/README.md

Instala la última versión del paquete PymCel utilizando pip. Asegúrate de tener pip actualizado para una instalación sin problemas.

```bash
pip install -U pymcel
```

--------------------------------

### Solving ODE with SciPy (solve_ivp)

Source: https://github.com/seap-udea/pymcel/blob/main/docs/examples/mcel_zuluaga-02-vectores_calculo.ipynb

Solves the ordinary differential equation for a simple harmonic oscillator using `scipy.integrate.solve_ivp`. This is an alternative to `odeint`. The output is then plotted similarly to the previous example.

```python
from scipy.integrate import solve_ivp
solucion=solve_ivp(fun=lambda t,Y:ode_simple(Y,t,k),t_span=[ts[0],ts[-1]],y0=Yos,t_eval=ts)

#Extraemos los valores de la función F
Fs=solucion.y[0]

plt.figure();
plt.plot(ts,Fs,marker='o',linewidth=0);

#--hide--
plt.xlabel("$t$");
plt.ylabel("$F(t)$");
```

--------------------------------

### Obtener estado inicial de la Tierra y Marte con SPICE

Source: https://github.com/seap-udea/pymcel/blob/main/docs/examples/mcel_zuluaga-13-doscuerpos_sintesis_y_especiales.ipynb

Este fragmento de código utiliza la biblioteca `spiceypy` para cargar efemérides y obtener los vectores de estado (posición y velocidad) de la Tierra y Marte con respecto al Sol en la fecha J2000.0. Requiere los archivos de datos `de430.tpc` y `de430.bsp`.

```python
from spiceypy import furnsh
furnsh("pymcel/data/de430.tpc")
furnsh("pymcel/data/de430.bsp")

#Fecha de referencia J2000.0
t0=0

#Parametro gravitacional del sistema
from spiceypy import bodvrd
mu=bodvrd("SUN","GM",1)[1][0]

#Vectores de estado de la Tierra y Marte
from spiceypy import spkezr
tierra,tluz=spkezr("EARTH_BARYCENTER",t0,"ECLIPJ2000","None","SUN")
rtierra0=tierra[:3]
vtierra0=tierra[3:]
marte,tluz=spkezr("MARS_BARYCENTER",t0,"ECLIPJ2000","None","SUN")
rmarte0=marte[:3]
vmarte0=marte[3:]

print("Estado inicial de la Tierra (SPICE):")
print(f"	Posición: {rtierra0}")
print(f"	Velocidad: {vtierra0}")
print("Estado inicial de Marte (SPICE):")
print(f"	Posición: {rmarte0}")
print(f"	Velocidad: {vmarte0}")
```

--------------------------------

### Set Parameters and Initial Conditions (Python)

Source: https://github.com/seap-udea/pymcel/blob/main/docs/examples/mcel_zuluaga-21-lagrangiano_fuerzas_centrales.ipynb

Initializes system parameters such as gravitational constant 'mu', potential exponent 'n', initial radial position 'r', angle 'teta', and their velocities. It also calculates the specific angular momentum 'h' and specific energy 'epsilon', and defines the time points for integration.

```python
#Parámetros del potencial
mu=1
n=1.1

#Condiciones iniciales
r=1.0
teta=0.0
r_dot=0.0
teta_dot=0.5

#Momento angular específico
h=r**2*teta_dot

#Energía específica
from pymcel import Vfuerza,Vcen
epsilon=0.5*r_dot**2+Vfuerza(r,mu=mu,n=n)+Vcen(r,h=h)

#Tiempos de integración
Nt=1000
from numpy import linspace
ts=linspace(0.0,10.0,Nt)
```

--------------------------------

### Example 3D CRTBP Orbit near L1 (Python)

Source: https://github.com/seap-udea/pymcel/blob/main/docs/examples/mcel_zuluaga-19-trescuerpos_aplicaciones.ipynb

Demonstrates a periodic trajectory related to the L1 Lagrange point in 3D CRTBP. It calculates the L1 point using `scipy.optimize.bisect` and `pymcel.funcion_puntos_colineales`, then plots the orbit using `orbitas_crtbp3d`.

```python
#Propiedades del sistema
alfa=0.0121505856
xL1=bisect(funcion_puntos_colineales,0,1-2*alfa,args=(alfa,))
#Condiciones iniciales
ro=[0.8329,0,0.1304]
vo=[0,0.2437,0]
fig=orbitas_crtbp3d(alfa,ro,vo,
                T=3,Nt=100,
                xlim=(-0.1,1.1),ylim=(-0.5,0.5),zlim=(-0.5,0.5),
                xL=xL1,yL=0,zL=0
               )
```

--------------------------------

### CRTBP Orbits near L3 (Python)

Source: https://github.com/seap-udea/pymcel/blob/main/docs/examples/mcel_zuluaga-19-trescuerpos_aplicaciones.ipynb

Illustrates a family of trajectories near the L3 Lagrange point in the CRTBP. This example uses `scipy.optimize.bisect` to find the L3 point and `pymcel.funcion_puntos_colineales` to assist in the calculation. It then plots the orbit using `orbitas_crtbp`.

```python
#Propiedades del sistema
alfa=1e-4
from scipy.optimize import bisect
from pymcel import funcion_puntos_colineales
xL3=bisect(funcion_puntos_colineales,-2,-0.5,args=(alfa,)) 
#Condiciones iniciales
ro=[-1.112349859300,0,0]
vo=[0,+0.202041957868,0]
fig=orbitas_crtbp(alfa,ro,vo,
              T=180,Nt=500,
              xlim=(-1.2,1.2),ylim=(-1.2,1.2),
              xL=xL4,yL=yL4)
```

--------------------------------

### Instalar pymcel

Source: https://github.com/seap-udea/pymcel/blob/main/docs/examples/mcel_zuluaga-04-calculo_variacional.ipynb

Instala la librería pymcel usando pip. Se recomienda ejecutar este comando en entornos como Google Colab antes de usar las funcionalidades del cuaderno.

```python
!pip install -Uq pymcel

```

--------------------------------

### Preparar gráfico de diferencias de posición con SPICE

Source: https://github.com/seap-udea/pymcel/blob/main/docs/examples/mcel_zuluaga-13-doscuerpos_sintesis_y_especiales.ipynb

Este fragmento de código inicializa un bucle para calcular las diferencias de posición entre las predicciones de PyMCEL y SPICE a lo largo de un período sinódico de Marte. Utiliza `numpy.linspace` para generar los puntos de tiempo y `spiceypy.prop2b` para la propagación con SPICE.

```python
from numpy import linspace

#Tiempos por un período sinódico completo de Marte
ts=linspace(0,780*86400,100)
errores=[]

for t in ts:
    #Propaga la posición de la Tierra
    from spiceypy import prop2b
    prediccion_tierra=prop2b(mu,list(rtierra0)+list(vtierra0),t-t0)
    rtierra=prediccion_tierra[:3]
    vtierra=prediccion_tierra[3:]

    #Propaga la posición de Marte
    from spiceypy import prop2b
    prediccion_marte=prop2b(mu,list(rmarte0)+list(vmarte0),t-t0)
    rmarte=prediccion_marte[:3]
    vmarte=prediccion_marte[3:]
    
    #Calcula las coordenadas predichas
    dist_aprox,long_aprox,lat_aprox=reclat(rmarte-rtierra)
    
    #Obtiene las coordenadas reales de Marte y la Tierra

```

--------------------------------

### Example CRTBP Orbit near L4 (Python)

Source: https://github.com/seap-udea/pymcel/blob/main/docs/examples/mcel_zuluaga-19-trescuerpos_aplicaciones.ipynb

Demonstrates plotting a periodic trajectory near the L4 Lagrange point in the CRTBP, as found by Grebow (2006). It sets up system properties and initial conditions, then calls the `orbitas_crtbp` function to generate the plot.

```python
#Propiedades del sistema
alfa=0.0121505856
xL4=0.5-alfa
yL4=3**0.5/2
#Condiciones iniciales
ro=[0.6867,yL4,0]
vo=[0.1126,-0.2040,0]
from numpy import pi
fig=orbitas_crtbp(alfa,ro,vo,
              T=2*pi,Nt=1000,
              xlim=(-0.1,1.1),ylim=(-0.1,1.1),
              xL=xL4,yL=yL4)
```

--------------------------------

### Preparar el entorno SPICE en pymcel

Source: https://github.com/seap-udea/pymcel/blob/main/ejemplos/pymcel-primeros-pasos.ipynb

Carga todos los kernels de SPICE necesarios para utilizar datos de efemérides y posicionamiento. Este comando prepara el entorno para que `pymcel` pueda acceder y procesar datos de SPICE.

```python
pc.prepara_spice()
```

--------------------------------

### Integrar sistema completo N-cuerpos en Python

Source: https://github.com/seap-udea/pymcel/blob/main/docs/examples/mcel_zuluaga-13-doscuerpos_sintesis_y_especiales.ipynb

Integra el sistema completo de N cuerpos utilizando la función `ncuerpos_solucion` de pymcel. Se define un array de tiempos para la integración. La salida incluye las posiciones y velocidades de todas las partículas en cada instante de tiempo.

```python
from pymcel import ncuerpos_solucion
ts=[t0,t]
rs,vs,rps,vps,constantes=ncuerpos_solucion(sistema,ts)
```

--------------------------------

### Example 3D Halo Orbit around L2 (Python)

Source: https://github.com/seap-udea/pymcel/blob/main/docs/examples/mcel_zuluaga-19-trescuerpos_aplicaciones.ipynb

Shows a halo orbit around the L2 Lagrange point in the CRTBP. The L2 point is found using `scipy.optimize.bisect` and `pymcel.funcion_puntos_colineales`. The `orbitas_crtbp3d` function is then used to plot this 3D orbit with specified viewing angles.

```python
#Propiedades del sistema
alfa=0.0121505856
xL2=bisect(funcion_puntos_colineales,1-0.5*alfa,2,args=(alfa,))
#Condiciones iniciales
ro=[1.1003,0,0]
vo=[0,-0.3217,0.5973]
fig=orbitas_crtbp3d(alfa,ro,vo,
                T=10,Nt=100,
                xlim=(-0.1,1.1),ylim=(-0.5,0.5),zlim=(-0.5,0.5),
                xL=xL2,yL=0,zL=0,
                elevation=10,azimuth=-60
               )
```

--------------------------------

### Visualizar aproximación de dos cuerpos vs N-cuerpos en Python

Source: https://github.com/seap-udea/pymcel/blob/main/docs/examples/mcel_zuluaga-13-doscuerpos_sintesis_y_especiales.ipynb

Genera una visualización 3D comparando la solución N-cuerpos con la aproximación de dos cuerpos. Utiliza las funciones `ncuerpos_solucion`, `propaga_estado` y `plot_ncuerpos_3d` de pymcel, junto con numpy para la generación de puntos de tiempo y arrays. Dependencias: pymcel, numpy, matplotlib.

```python
from pymcel import propaga_estado
from numpy import linspace

#Ventana de integración
t0=0.0
T=5.0

#Solución al problema de los N cuerpos
from pymcel import ncuerpos_solucion
ts=linspace(t0,T,200)
rs_num,vs_num,rps_num,vps_num,constantes=ncuerpos_solucion(sistema,ts)

from numpy import zeros_like
rs_aprox=zeros_like(rs_num)
Nt=20
ts=linspace(t0,T,Nt)
for i,t in enumerate(ts):
    #Sistema A
    ra1,va1,ra2,va2,ravecA,vavecA=propaga_estado(sistemaA,t0,t)

    #Sistema B
    raA,vaA,ra3,va3,ravecB,vavecB=propaga_estado(sistemaB,t0,t) 
    
    rs_aprox[0,i]=ra1
    rs_aprox[1,i]=ra2
    rs_aprox[2,i]=ra3
    
#Grafíco
from pymcel import plot_ncuerpos_3d
fig=plot_ncuerpos_3d(rs_num,vs_num)
ax=fig.gca()

ax.plot(rs_aprox[0,:Nt,0],rs_aprox[0,:Nt,1],rs_aprox[0,:Nt,2],'b+')
ax.plot(rs_aprox[1,:Nt,0],rs_aprox[1,:Nt,1],rs_aprox[1,:Nt,2],'r+')
ax.plot(rs_aprox[2,:Nt,0],rs_aprox[2,:Nt,1],rs_aprox[2,:Nt,2],'g+')

plt.show()
```

--------------------------------

### Acceder al directorio de instalación de pymcel

Source: https://github.com/seap-udea/pymcel/blob/main/ejemplos/pymcel-primeros-pasos.ipynb

Este ejemplo demuestra cómo acceder a la variable global `ROOTDIR` de `pymcel`, que indica la ruta de instalación del paquete en el sistema de archivos. Es útil para localizar archivos o recursos del paquete.

```python
pc.ROOTDIR
```

--------------------------------

### Propagate Two-Body State Over Time with PyMCEL

Source: https://context7.com/seap-udea/pymcel/llms.txt

Propagates the state of a two-body system from an initial time to a final time using an analytical solution. Requires 'pymcel' and 'numpy'. Defines the system with masses, initial positions, and velocities. The example sets up a system with two bodies and prints the initial and final positions and shapes of the output arrays.

```python
import pymcel as pc
import numpy as np

# Definir sistema de dos cuerpos
sistema = [
    dict(m=1.0, r=np.array([-0.5, 0.0, 0.0]), v=np.array([0.0, -0.5, 0.0])),
    dict(m=1.0, r=np.array([0.5, 0.0, 0.0]), v=np.array([0.0, 0.5, 0.0]))
]

# Note: The actual propagation call is missing in the provided snippet.
```

--------------------------------

### Descargar datos con pymcel

Source: https://github.com/seap-udea/pymcel/blob/main/docs/examples/mcel_zuluaga-19-trescuerpos_aplicaciones.ipynb

Utiliza la función 'obtiene_datos()' de la librería 'pymcel' para descargar los kernels y archivos de datos necesarios para las simulaciones de mecánica celeste. Estos datos son esenciales para el correcto funcionamiento de los algoritmos.

```python
import pymcel as pc
pc.obtiene_datos()
```

--------------------------------

### Prepare and Query Local SPICE Data with PyMCEL

Source: https://context7.com/seap-udea/pymcel/llms.txt

Prepares the local SPICE environment by downloading necessary kernels and then queries SPICE for celestial body position and velocity data. Requires the 'pymcel' library. Handles single time points and time ranges.

```python
import pymcel as pc

# Preparar el entorno SPICE (descarga kernels si es necesario)
pc.prepara_spice(verbose=True)
# Salida:
# Cargando todos los kernels de SPICE...
# El entorno está listo para usar los datos de SPICE.

# Consultar posición de la Tierra usando SPICE
data, tiempos, df = pc.consulta_spice(
    id='399',
    location='@0',
    epochs='2024-01-01 12:00:00'
)
print(f"Posición (m): {data[:3]}")
print(f"Velocidad (m/s): {data[3:]}")

# Consulta con rango de tiempos
data, tiempos, df = pc.consulta_spice(
    id='399',
    location='@0',
    epochs=dict(start='2024-01-01', stop='2024-01-31', step='1d')
)
print(f"Forma de datos: {data.shape}")  # (31, 6)
print(df.head())
```

--------------------------------

### Calcular Energía en el Problema de N Cuerpos (Tierra-Luna) con Python y SpiceyPy

Source: https://github.com/seap-udea/pymcel/blob/main/docs/examples/mcel_zuluaga-09-ncuerpos_general.ipynb

Este snippet calcula la energía potencial, cinética y total para el sistema Tierra-Luna en un instante dado. Utiliza la biblioteca SpiceyPy para obtener las posiciones y velocidades de los cuerpos celestes y las constantes físicas necesarias. Las dependencias incluyen 'spiceypy' y 'numpy'.

```python
import spiceypy as spy

#Constante de gravitación universal 
G=6.67e-20 # km^3 / kg s^2

#Asumimos un tiempo cualquiera, en este caso t=J2000.0
tef=0

#Carga kernels con posiciones (bsp) y masas (tpc)
spy.furnsh('pymcel/data/de430.bsp')
spy.furnsh('pymcel/data/de430.tpc')

#Parámetro gravitacional, posiciones y velocidades
mutierra=spy.bodvrd("EARTH","GM",1)[1][0]
tierra,tluz=spy.spkezr("EARTH",tef,
                       "ECLIPJ2000","None","EARTH_BARYCENTER")
rtierra=tierra[:3]
vtierra=tierra[3:]

muluna=spy.bodvrd("MOON","GM",1)[1][0]
luna,tluz=spy.spkezr("MOON",tef,
                     "ECLIPJ2000","None","EARTH_BARYCENTER")
rluna=luna[:3]
vluna=luna[3:]

#Masas
mtierra=mutierra/G
mluna=muluna/G

#Energía potencial
from numpy.linalg import norm
U=-G*mluna*mtierra/norm(rtierra-rluna)

#Energía cinética
K=0.5*mutierra*norm(vtierra)**2+0.5*mluna*norm(vluna)**2

#Energía mecánica total
E=K+U

print(f"Energía potencial U: {U} kg km^2/s^2")
print(f"Energía cinética K: {K} kg km^2/s^2")
print(f"Energía mecánica total E: {E} kg km^2/s^2")
```

--------------------------------

### Instalar y Verificar PyMCel

Source: https://context7.com/seap-udea/pymcel/llms.txt

Instala o actualiza el paquete PyMCel usando pip y verifica la versión instalada y el directorio raíz del paquete. Requiere acceso a internet para la descarga desde PyPI.

```python
import pymcel as pc
import numpy as np

# Instalar o actualizar pymcel
# pip install -U pymcel

# Verificar versión instalada
print(f"Versión: {pc.version}")
# Salida: Bienvenido a PyMCel v0.9.4 ¡al infinito y más allá!
# Salida: Versión: 0.9.4

# Directorio de instalación
print(f"Directorio: {pc.ROOTDIR}")
```

--------------------------------

### Calcular Energía del Sistema Solar (Sol-Júpiter-Saturno) con Python y SpiceyPy

Source: https://github.com/seap-udea/pymcel/blob/main/docs/examples/mcel_zuluaga-09-ncuerpos_general.ipynb

Este snippet calcula las energías cinética, potencial y total promedio para el sistema Sol-Júpiter-Saturno a lo largo de un período de tiempo. Utiliza SpiceyPy para obtener posiciones y velocidades, y NumPy para cálculos y promedios. Las dependencias incluyen 'spiceypy' y 'numpy'.

```python
#Constante de gravitación universal 
G=6.67e-20 # km^3 / kg s^2

#Carga kernels con posiciones (bsp) y masas (tpc)
import spiceypy as spy
spy.furnsh('pymcel/data/de430.bsp')
spy.furnsh('pymcel/data/de430.tpc')

#Número de valores de tiempo
Nt=100

#Lista de tiempos en los que calcularemos el virial:
#Tomamos 60 años que es aprox. 2 veces el período de Saturno
from numpy import linspace
tefs=linspace(0.0,60*356.25*86400,Nt)

#Masas de los cuerpos
msol=spy.bodvrd("SUN","GM",1)[1][0]/G
mjupiter=spy.bodvrd("JUPITER_BARYCENTER","GM",1)[1][0]/G
msaturno=spy.bodvrd("SATURN_BARYCENTER","GM",1)[1][0]/G
```

```python
from numpy.linalg import norm
Ks=[]
Us=[]
Es=[]
for tef in tefs:
    #Posiciones, velocidades, energías cinéticas
    sol,tluz=spy.spkezr("SUN",tef,
                        "ECLIPJ2000","None","SSB")
    rsol=sol[:3]
    vsol=sol[3:]
    K_sol=0.5*msol*norm(vsol)**2
    
    jupiter,tluz=spy.spkezr("JUPITER_BARYCENTER",tef,
                            "ECLIPJ2000","None","SSB")
    rjupiter=jupiter[:3]
    vjupiter=jupiter[3:]
    K_jup=0.5*mjupiter*norm(vjupiter)**2
    
    saturno,tluz=spy.spkezr("SATURN_BARYCENTER",tef,
                            "ECLIPJ2000","None","SSB")

    rsaturno=saturno[:3]
    vsaturno=saturno[3:]
    K_sat=0.5*msaturno*norm(vsaturno)**2

    #Distancias entre los cuerpos
    r_jup_sol=norm(rjupiter-rsol)
    r_sat_sol=norm(rsaturno-rsol)
    r_jup_sat=norm(rjupiter-rsaturno)
    
    #Energia potencial
    U_jup_sol=-G*mjupiter*msol/r_jup_sol
    U_sat_sol=-G*msaturno*msol/r_sat_sol
    U_jup_sat=-G*mjupiter*msaturno/r_jup_sat

    #Energía cinética, potencial y mecánica
    Ktot=K_sol+K_jup+K_sat
    U=U_jup_sol+U_sat_sol+U_jup_sat
    E=Ktot+U
    
    #Guarda valores en la lista
    Ks+=[Ktot]
    Us+=[abs(U)]
    Es+=[abs(E)]

from numpy import mean
print(f"<K> = {mean(Ks)} kg km^2/s^2")
print(f"<U>/2 = -{mean(Us)/2} kg km^2/s^2")
print(f"E = -{mean(Es)} kg km^2/s^2")
```

--------------------------------

### Verificar estado propagado con SPICE

Source: https://github.com/seap-udea/pymcel/blob/main/docs/examples/mcel_zuluaga-13-doscuerpos_sintesis_y_especiales.ipynb

Compara las predicciones obtenidas con `pymcel` contra los valores de estado de la Tierra y Marte calculados directamente por SPICE para el mismo instante futuro. Esto permite evaluar la precisión del método de propagación aproximado.

```python
#Vectores de estado de la Tierra y Marte
tierra,tluz=spkezr(t,"EARTH_BARYCENTER","ECLIPJ2000","None","SUN")
marte,tluz=spkezr(t,"MARS_BARYCENTER","ECLIPJ2000","None","SUN")

print("Estado propagado de la Tierra (SPICE):")
print(f"	Posición: {tierra[:3]}")
print(f"	Velocidad: {tierra[3:]}")
print("Estado inicial de Marte (SPICE):")
print(f"	Posición: {marte[:3]}")
print(f"	Velocidad: {marte[3:]}")
```

--------------------------------

### Instalar pymcel en Google Colab

Source: https://github.com/seap-udea/pymcel/blob/main/ejemplos/pymcel-descarga-kernels.ipynb

Instala la última versión de la biblioteca pymcel usando pip. Este comando es específico para entornos como Google Colab donde se requiere acceso de superusuario para la instalación global.

```python
import sys
if 'google.colab' in sys.modules:
    !sudo pip install -Uq pymcel
```

--------------------------------

### Set Up and Solve Conical Pendulum Motion (Python)

Source: https://github.com/seap-udea/pymcel/blob/main/docs/examples/mcel_zuluaga-22-hamiltoniano.ipynb

This script sets up the parameters, initial conditions, and time span for simulating a conical pendulum. It then uses `scipy.integrate.odeint` to numerically solve the equations of motion defined by `edm_penduloconico_hamiltoniano`. Finally, it extracts the calculated positions and momenta.

```python
#Propiedades
g=9.81
L=1.0
m=1.0

#Factores de conversión
from numpy import pi,linspace
grados=pi/180

#Condiciones iniciales
teta_0=15*grados
fi_0=0.0*grados
pteta_0=0.0
alfa_fi=0.8

#Tiempos de integración
from numpy import pi,sqrt,linspace
T=2*pi*sqrt(L/g)
Nt=300
ts=linspace(0,3*T,Nt)

#Solución numérica
from scipy.integrate import odeint
solucion=odeint(edm_penduloconico_hamiltoniano,
                [teta_0,fi_0,pteta_0],ts,
                args=(alfa_fi,L,m))

#Extrae momentos
tetas=solucion[:,0]
fis=solucion[:,1]
ptetas=solucion[:,2]

#Coordenadas cartesianas
from numpy import sin,cos
xs=L*sin(tetas)*cos(fis)
ys=L*sin(tetas)*sin(fis)
z s=-L*cos(tetas)
```

--------------------------------

### Simular y Graficar Sistema de Tres Cuerpos en Python

Source: https://github.com/seap-udea/pymcel/blob/main/docs/examples/mcel_zuluaga-11-doscuerpos_motivacion_a_orbita_osculatriz.ipynb

Simula la trayectoria 3D de un sistema de tres cuerpos y grafica el resultado. Utiliza la función 'ncuerpos_solucion' para obtener las posiciones y velocidades, y 'plot_ncuerpos_3d' para la visualización. Requiere la librería 'pymcel'.

```python
sistema=[
    dict(
        m=10.0,
        r=[1,0,0],
        v=[0,1,0.5]),
    dict(
        m=1.0,
        r=[1.5,0,0],
        v=[0,-3,1],
    ),
    dict(
        m=0.1,
        r=[-1,0,0],
        v=[0,3,0],
    )
]

from numpy import linspace

#Solución
from pymcel import ncuerpos_solucion 
rs,vs,rps,vps,constantes=ncuerpos_solucion(sistema,
                                           linspace(0.0,10.0,200))

#Gráfica en el sistema de referencia inercial original
from pymcel import plot_ncuerpos_3d
fig=plot_ncuerpos_3d(rs,vs)
```

```python
#Gráfica en el sistema de referencia del centro de masa
from pymcel import plot_ncuerpos_3d
fig=plot_ncuerpos_3d(rps,vps);
```

--------------------------------

### Configurar Condiciones Iniciales y Tiempos para Péndulo Elástico

Source: https://github.com/seap-udea/pymcel/blob/main/docs/examples/mcel_zuluaga-20-lagrangiano_generalidades.ipynb

Establece las propiedades del sistema (longitud, constante del resorte, masa) y las condiciones iniciales (posiciones y velocidades generalizadas) para el péndulo elástico. También define los puntos de tiempo en los que se calculará la solución.

```python
#Propiedades del sistema
L=1.0
k=20.0
m=1.0

#Condiciones iniciales
from numpy import pi
y=[pi/3,0.1,0.0,0.0]

#Tiempos de integración
from numpy import linspace
Nt=200
ts=linspace(0,10,Nt)
```

--------------------------------

### Obtener la versión de pymcel

Source: https://github.com/seap-udea/pymcel/blob/main/ejemplos/pymcel-primeros-pasos.ipynb

Muestra cómo acceder a la variable incorporada `version` del paquete `pymcel` para verificar la versión instalada. Esta es una forma rápida de confirmar la versión del software en uso.

```python
pc.version
```

--------------------------------

### Leer y graficar objetos cercanos a Júpiter con Tisserand Parameter - Python

Source: https://github.com/seap-udea/pymcel/blob/main/docs/examples/mcel_zuluaga-19-trescuerpos_aplicaciones.ipynb

Este script de Python lee datos orbitales de cuerpos cercanos a Júpiter desde un archivo CSV, calcula el semieje mayor 'a', y genera un gráfico de dispersión de 'a' vs. 'e'. Incluye la superposición de curvas de igual valor del parámetro de Tisserand (T_J=2 y T_J=3) para ayudar en la clasificación. Las dependencias incluyen pandas y matplotlib. Los inputs son un archivo CSV y los parámetros Tisserand. Las salidas son un gráfico interactivo.

```python
#Lee objetos
from pandas import read_csv
datos=read_csv("pymcel/data/JupiterComets.csv",
               index_col=["full_name"],skiprows=1)
datos["a"]=datos["q"]/(1-datos["e"])

#Gráfico de dispersión
fig=plt.figure(figsize=(5,5))
ax=fig.gca()
s=ax.scatter(datos["a"],datos["e"],
             s=5,c='b')

#Contorno de paraámetro de Tisserand
aJ=5.2044
from numpy import linspace,sqrt
aes=linspace(3.2,5.5,100)
T=2
eTs=sqrt(1-1/(4*(aes/aJ)**3)*(T*(aes/aJ)-1)**2)
ax.plot(aes,eTs,'k-',lw=2,label=f"$T_J={T}$")
T=3
eTs=sqrt(1-1/(4*(aes/aJ)**3)*(T*(aes/aJ)-1)**2)
ax.plot(aes,eTs,'r-',lw=3,label=f"$T_J={T}$")

#Decoración
ax.set_xlim((3.2,5.5))
ax.set_ylim((0.0,1.0))
ax.set_xlabel("$a$ [ua]")
ax.set_ylabel("$e$")
ax.legend(loc='lower left')
fig.tight_layout()
```