### Install Celmech from Source

Source: https://github.com/shadden/celmech/blob/master/docs/quickstart.md

Install the celmech Python package after cloning the repository. Navigate into the top celmech directory before running this command.

```default
python setup.py install
```

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

### Setup REBOUND Simulation

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/CoordinatesConvertingToFromNbody.ipynb

Initializes a REBOUND simulation with a central mass and two smaller bodies, then moves to the center of mass frame and plots the orbits.

```python
import celmech
import rebound
%matplotlib inline
```

```python
sim = rebound.Simulation()
sim.add(m=1.)
sim.add(m=1e-3, a=1, e=0.4, f = 0.3)
sim.add(m=1e-3, a=1.5, e=0.2, f = 1)
sim.move_to_com()
rebound.OrbitPlot(sim)
```

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

### Install Celmech using Pip

Source: https://github.com/shadden/celmech/blob/master/docs/quickstart.md

Install the celmech Python package directly from PyPI using pip. This is the recommended method for most users.

```default
pip install celmech
```

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

### Import necessary libraries

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/OverstableLibrations.ipynb

Imports required libraries for numerical simulations, resonance modeling, and plotting. Ensure these are installed before running.

```python
import numpy as np
import rebound as rb
from celmech.numerical_resonance_models import PlanarResonanceEquations
import matplotlib.pyplot as plt
%matplotlib inline
```

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

### Laplace-Lagrange Secular System Setup

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/Laplace-LagrangeSecularTheory.ipynb

This snippet illustrates the mathematical formulation for setting up a Laplace-Lagrange secular system, relevant for N-body simulations. It involves complex terms related to masses, gravitational constants, and orbital parameters.

```latex
\begin{array}{c}
\frac{2}{ \Lambda_4 } \left( \frac{ m_4 \mu_4 }{2} - \frac{ 2 C_{(0,0,0,0,0,0)}^{(0,0,1,0),(0,0)}_{2,4} \alpha G M_4 m_2 m_4 \mu_4}{3} - \frac{ 2 C_{(0,0,0,0,0,0)}^{(0,0,1,0),(0,0)}_{3,4} \alpha G M_4 m_3 m_4 \mu_4}{3} \right)
\end{array}
```

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

### Import necessary libraries

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/Laplace-LagrangeSecularTheory.ipynb

Imports essential libraries for N-body simulations, plotting, and secular theory calculations. Ensure these libraries are installed before running.

```python
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from celmech.nbody_simulation_utilities import set_time_step,align_simulation
from celmech.nbody_simulation_utilities import get_simarchive_integration_results
import rebound as rb
from os.path import isfile
from sympy import init_printing
init_printing()
```

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

### Initialize SecularSystemSimulation from REBOUND Archive

Source: https://github.com/shadden/celmech/blob/master/docs/secular.md

Initializes a SecularSystemSimulation object from the first snapshot of a REBOUND simulation archive. This is the starting point for integrating nonlinear secular equations of motion.

```python
# grab first simulation archive snapshot
sim0=sa[0]
```

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

### Set planet masses and migration parameters

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/OverstableLibrations.ipynb

Configures planet masses, migration timescales (tau_alpha), and the parameter 'p'. The instability requires p != 0. K1 and K2 are set equal in this example.

```python
res_eqs.m1 = 1e-4
res_eqs.m2 = 1.5e-4
res_eqs.tau_alpha = 3e5
res_eqs.p = 1
def set_Kvals(res_eqs,K):
    res_eqs.K1=K
    res_eqs.K2=K
```

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

### Setup Poincare Hamiltonian with Disturbing Terms

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/AddingDisturbingFunctionTerms.ipynb

Sets up a Poincare Hamiltonian and adds Mean Motion Resonance (MMR) terms and secular terms. This is a prerequisite for transforming to mean variables.

```python
max_order = 2
sim=get_sim()
poincare_variables = Poincare.from_Simulation(sim)
Hp = PoincareHamiltonian(poincare_variables)
Hp.add_MMR_terms(p=3, q=1, max_order=max_order, indexIn=1, indexOut=2)
for idx1 in range(1,sim.N):
    for idx2 in range(idx1+1,sim.N):
        Hp.add_secular_terms(max_order=max_order, indexIn=idx1, indexOut=idx2)
```

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

### Initialize Hamiltonian and Add Terms (2nd Order)

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/AddingDisturbingFunctionTerms.ipynb

Initializes a PoincareHamiltonian object and adds Mixed Resonance (MMR) and secular terms up to the second order. This setup is used for integrating the system and observing orbital evolution.

```python
sim=get_sim()
poincare_variables = Poincare.from_Simulation(sim)

max_order = 2
Hp = PoincareHamiltonian(poincare_variables)
Hp.add_MMR_terms(p=3, q=1, max_order=max_order, indexIn=1, indexOut=2)
for idx1 in range(1,sim.N):
    for idx2 in range(idx1+1,sim.N):
        Hp.add_secular_terms(max_order=max_order, indexIn=idx1, indexOut=idx2)
chi = FirstOrderGeneratingFunction(poincare_variables)
chi.add_zeroth_order_term()
chi.add_MMR_terms(p=2,q=1,l_max=1, indexIn=1, indexOut=2)
chi.add_MMR_terms(p=4,q=1,l_max=1, indexIn=1, indexOut=2)
chi.osculating_to_mean()
```

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

### Initialize Laplace-Lagrange System from Simulation

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/Laplace-LagrangeSecularTheory.ipynb

Load simulation results and initialize the Laplace-Lagrange system. This requires a simulation archive and the integration results.

```python
nb_results = get_simarchive_integration_results(sa)
sim = sa[0]
lsys = LaplaceLagrangeSystem.from_Simulation(sim)
sec_soln = lsys.secular_solution(nb_results['time'])
```

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

### Set up Simulation with GR Potential

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/Laplace-LagrangeTheoryWithGR.ipynb

Initializes a `rebound` simulation, adds particles, moves to the center of mass, and configures the integrator. It then uses `reboundx` to load and add the `gr_potential` force, setting the speed of light parameter.

```python
import rebound
import numpy as np
from celmech.nbody_simulation_utilities import get_simarchive_integration_results
from celmech.secular import LaplaceLagrangeSystem
import matplotlib.pyplot as plt
import reboundx
from reboundx import constants

# Set up a rebound simulation
sim = rebound.Simulation()
sim.add(m=1)
sim.add(m=1e-5, a=0.1, e=0.001)
sim.add(m=1e-5, a=0.2, e=0.001, pomega=0.5*np.pi)
sim.move_to_com()
sim.integrator = 'whfast'
sim.dt = sim.particles[1].P / 25.

# Use reboundx to add GR. 
# This example uses the simplest implementation, `gr_potential`.
rebx = reboundx.Extras(sim)
gr = rebx.load_force("gr_potential")
rebx.add_force(gr)
gr.params["c"] = constants.C
```

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

### Set up and run REBOUND simulation

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/Laplace-LagrangeSecularTheory.ipynb

Sets up a REBOUND simulation, initializes particles with specified masses and orbital elements, and integrates the system. The simulation is saved to a binary file. This block includes error handling for existing simulation files.

```python
try:
    sa = rb.Simulationarchive('sa_laplace_lagrange.bin')
except RuntimeError:
    np.random.seed(1)
    sim = rb.Simulation()
    sim.add(m=1)
    m = 5e-5
    sigma_e = 0.02
    sigma_inc = sigma_e / 2
    for a in np.linspace(1,8,4):
        sim.add(
            m=m / np.sqrt(a),
            a=a,
            e=np.random.rayleigh(sigma_e),
            inc=np.random.rayleigh(sigma_inc),
            pomega = "uniform",
            Omega = "uniform",
            l="uniform"
        )
    sim.integrator = 'whfast'
    set_time_step(sim,1/40)
    sim.ri_whfast.safe_mode = 0
    sim.move_to_com()
    align_simulation(sim)
    
    lsys = LaplaceLagrangeSystem.from_Simulation(sim)
    eigvals = lsys.eccentricity_eigenvalues()
    Tsec = np.max(np.abs(2 * np.pi / eigvals))
    sim.save_to_file('sa_laplace_lagrange.bin',interval=Tsec/2048,delete_file=True)
    sim.integrate(0.5*Tsec)
    sa = rb.Simulationarchive('sa_laplace_lagrange.bin')
```

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

### Disturbing Function Argument Example 3

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/DisturbingFunctionArguments.ipynb

This example illustrates a disturbing function argument with a focus on specific orbital elements and their interactions.

```mathematica
C_{(0,0,0,2,-2,0)}_{1,2})__{(0,0,0,0),(0,0)}(\alpha⋅cos(2⋅Ω₁ - 2⋅varpi₂)⋅│X₂│ 

     2
⋅│Y₁│ 
```

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

### Disturbing Function Argument Example 1

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/DisturbingFunctionArguments.ipynb

This example shows a specific form of a disturbing function argument involving trigonometric terms and symbolic variables.

```mathematica
C_{(0,0,2,0,-2,0)}_{1,2})__{(0,0,0,0),(0,0)}(\alpha⋅cos(2⋅Ω₁ - 2⋅varpi₁)⋅│X₁│ 

     2
⋅│Y₁│ 
```

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

### Initialize Rebound simulation

Source: https://github.com/shadden/celmech/blob/master/docs/quickstart.md

Sets up a Rebound simulation with a star and two Earth-mass planets in a 3:2 period ratio. Includes plotting the initial configuration.

```python
sim = rb.Simulation()
sim.add(m=1,hash='star')
sim.add(m=3e-6,P = 1, e = 0.03,l=0)
sim.add(m=3e-6,P = 3/2, e = 0.03,l=np.pi/5,pomega = np.pi)
sim.move_to_com()
rb.OrbitPlot(sim,periastron=True)
```

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

### Disturbing Function Argument Example 2

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/DisturbingFunctionArguments.ipynb

This example presents another variation of a disturbing function argument, featuring different angular relationships and variable combinations.

```mathematica
C_{(0,0,1,1,-2,0)}_{1,2})__{(0,0,0,0),(0,0)}(\alpha⋅cos(-2⋅Ω₁ + varpi₁ + varpi

                 2
₂)⋅│X₁│⋅│X₂│⋅│Y₁│ 
```

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

### Setting up a REBOUND simulation

Source: https://github.com/shadden/celmech/blob/master/docs/canonical_transformations.md

Initialize a REBOUND simulation with a star and two planets, setting their masses and orbital parameters.

```python
sim = rebound.Simulation()
sim.add(m=1)
sim.add(m=1e-5,P = 1, e=0.04)
sim.add(m=1e-5,P = (3/2) * ( 1 + .03), e=0.02,pomega = np.pi/2,l=0)
sim.move_to_com()
```

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

### Figure Size Example

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/OverstableLibrations.ipynb

Indicates the size and structure of a matplotlib figure, showing it contains 2 axes.

```text
<Figure size 504x504 with 2 Axes>
```

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

### Plot Output Example

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/OverstableLibrations.ipynb

Represents the output range for a plot's y-axis, indicating values between -180.0 and 180.0.

```text
(-180.0, 180.0)
```

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

### Initialize Simulation Parameters

Source: https://github.com/shadden/celmech/blob/master/docs/quickstart.md

Sets up time arrays, result keys, and particle lists for both N-body and CelMech simulations. Define the number of output steps and the time span for the simulation.

```python
# Here we define the times at which we'll get simulation outputs
Nout = 150
times = np.linspace(0 , 3e3, Nout) * sim.particles[1].P

# These are the quantites we'll track in our rebound and celmech integrations
keys = ['l1','l2','pomega1','pomega2','e1','e2','a1','a2']

# These dictionaries will hold our results
rebound_results= {key:np.zeros(Nout) for key in keys}
celmech_results= {key:np.zeros(Nout) for key in keys}

# These are the lists of particles in both simulations
# for which we'll save quantities.
rb_particles = sim.particles
cm_particles = pvars.particles
```

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

### Get Current State of Hamiltonian

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/CanonicalTransformations.ipynb

Retrieve the current state of the Hamiltonian object, which may be useful for subsequent calculations or analysis.

```python
cart_state = new_ham_cart.state
```

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

### Set up REBOUND Simulation with Planets

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/AddingDisturbingFunctionTerms.ipynb

Initializes a REBOUND simulation with a star and three Earth-mass planets. Sets up orbital parameters like mass, period, eccentricity, and inclination. Includes moving to the center of mass and setting the integrator. Use this to create a multi-planet system for simulation.

```python
import rebound
def get_sim():
    np.random.seed(1)
    sim = rebound.Simulation()
    sim.add(m=1)
 
    sim.G = 4*np.pi*np.pi
    mass=3e-6
    e=0.01
    inc=0.05
    
    sim.add(m=mass,P=1, e=0.02, pomega=0, theta=np.pi, inc=0.02)
    sim.add(m=mass,P=3/2, e=0.03, pomega=0, theta=np.pi, inc=0.03)
    sim.add(m=mass,P=3.2, e=0, pomega=0, theta=0)
    sim.move_to_com() 
    
    sim.dt = sim.particles[1].P / 30.
    sim.integrator='whfast'
    return sim
```

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

### Initialize simulation with two planets

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/FirstOrderGeneratingFunction.ipynb

Sets up a REBOUND simulation with a central mass and two smaller planets, positioned near a 3:2 mean motion resonance. Includes a plot of the initial configuration.

```python
def get_sim():
    sim = rb.Simulation()
    sim.add(m=1)
    sim.add(m=1e-5,P = 1, e=0.04)
    sim.add(m=1e-5,P = (3/2) * ( 1 + .03), e=0.02,pomega = np.pi/2,l=0)
    sim.move_to_com()
    return sim
sim = get_sim()
rb.OrbitPlot(sim);
```

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

### Count Resonance Terms

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/DisturbingFunctionArguments.ipynb

Get the total number of resonance terms generated by `list_resonance_terms` by checking the length of the returned list.

```python
len(res_terms)
```

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

### Initialize and Run REBOUND Simulation

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/OverstableLibrations.ipynb

Sets up and runs a REBOUND simulation to observe the system's evolution, including dissipation effects. The simulation archive is configured for a specific duration and output frequency.

```python
from celmech.nbody_simulation_utilities import set_time_step, get_simarchive_integration_results
Tfin = 1.85 * res_eqs.tau_alpha
Nout = 500
sim.integrator='whfast'
set_time_step(sim,1/50.)
sim.automateSimulationArchive("overstable_libration.sa",Tfin/Nout,deletefile=True)
sim.integrate(Tfin)
```

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

### Import Poincare Functions from Celmech

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/AddingDisturbingFunctionTerms.ipynb

Imports necessary classes for working with Poincare variables and Hamiltonians from the celmech package. Ensure celmech is installed.

```python
from celmech import Poincare, PoincareHamiltonian
```

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

### Get Number of Degrees of Freedom

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/ResonantChain.ipynb

Retrieve the initial and reduced number of degrees of freedom after canonical transformation. This helps understand the extent of reduction achieved.

```python
(pham.N_dof,kam.N_dof)
```

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

### Initialize Hamiltonian System

Source: https://github.com/shadden/celmech/blob/master/docs/hamiltonian.md

Initialize a Hamiltonian system using symbolic expressions, numerical parameters, and a phase space state. Displays the symbolic Hamiltonian.

```python
import numpy as np
from sympy import symbols,cos,sqrt
from celmech.hamiltonian import Hamiltonian,PhaseSpaceState
p,theta,m,g,l=symbols("p,theta,m,g,l")
H = p*p/2/(m*l*l) + (m*g*l) * (1 - cos(theta))
theta0,p0=np.pi/2,0
state = PhaseSpaceState([theta,p],[theta0,p0])
pars=dict()
pars[g]=9.8
pars[l] = 1
pars[m] = 0.15
ham = Hamiltonian(H,pars,state)
ham.H
```

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

### Polar to Cartesian Transformation

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/CanonicalTransformations.ipynb

This code initializes a polar to Cartesian canonical transformation, specifying which QP pairs to transform by their indices. Indices start at 0.

```python
to_cart_tr = CanonicalTransformation.polar_to_cartesian(new_ham.qp_vars,indices=[2,3,4,5])
```

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

### Initialize Hamiltonian and Phase Space

Source: https://github.com/shadden/celmech/blob/master/docs/canonical_transformations.md

Sets up a Hamiltonian object and its associated phase space state using sympy symbols and numpy arrays. This is the initial step before applying transformations.

```python
from celmech import CanonicalTransformation, Hamiltonian, PhaseSpaceState
import sympy as sp
import numpy as np
p1,p2,q1,q2 = sp.symbols("p(1:3),q(1:3)")
H = p1**2/2 + p2 +  p2**4 + sp.cos(3*q1 + 2 * q2)
state=PhaseSpaceState([q1,q2,p1,p2],np.random.uniform(-0.5,0.5,4))
ham = Hamiltonian(H,{},state)
```

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

### Integrating Initial Conditions with Celmech

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/PlanarResonanceEquations.ipynb

Sets up and integrates initial conditions using the celmech library. This involves defining initial orbital elements, creating a simulation object, and then integrating the equations of motion.

```python
Delta0 = 0.007
e1_0 = 0.01
e2_0 = 0.01
theta1_0 = np.random.uniform(0,2*np.pi)
sim = get_sim(res_eqs,Delta0,e1_0,e2_0,theta1_0)
dyvars = res_eqs.dyvars_from_rebound_simulation(sim)
times = np.linspace(0,1e5,512)
soln = res_eqs.integrate_initial_conditions(dyvars,times,dissipation=True)
els_res_eqs = soln['orbital_elements']

# Store Delta
Pratio = (els_res_eqs['a2']/els_res_eqs['a1'])**(1.5)
Pratio0 = res_eqs.j / (res_eqs.j-res_eqs.k) 
els_res_eqs['Delta'] = Pratio/Pratio0  - 1
```

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

### Set up Outer Solar System Simulation

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/FrequencyAnalysis.ipynb

Initializes a REBOUND simulation for the outer solar system, adding major planets and setting simulation parameters.

```python
import rebound as rb
import numpy as np
from celmech.nbody_simulation_utilities import set_time_step
sim = rb.Simulation()
sim.units = ("Msun","AU","yr")
bodies = ["Sun","Jupiter","Saturn","Uranus","Neptune"]
for body in bodies:
    sim.add(body)
sim.move_to_com()
sim.integrator = 'whfast'
sim.ri_whfast.safe_mode=0
set_time_step(sim,1/20)
```

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

### Default Output of df_coefficient_Ctilde

Source: https://github.com/shadden/celmech/blob/master/docs/disturbing_function.md

The output of df_coefficient_Ctilde for a given set of arguments is a dictionary representing a sum of Laplace coefficients and their derivatives. This example shows the typical format.

```default
{(1, (1.5, 3, 0)): 2.5, (2, (1.5, 3, 1)): 0.25, ('indirect', 1): 0}
```

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

### Get Integration Results with Dense Output

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/AndoyerResonanceCapture.ipynb

Performs numerical integration and returns the solution object with dense output enabled, allowing for interpolation at arbitrary times.

```python
def get_integration_results(pars):
    x0,y0,delta0,tau,tfin = pars
    fun = lambda t,y: Nflow(*y,tau).reshape(-1)
    Dfun  = lambda t,y: Njac(*y,tau)
    soln = solve_ivp(fun,(0,tfin),y0=[y0,x0,delta0],jac=Dfun,method="Radau",dense_output=True)
    return soln
```

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

### Initialize Phase Space State and Hamiltonian Parameters

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/SimplePendulum.ipynb

Sets up the initial conditions for the pendulum (angle and momentum) and defines the physical parameters (gravity, length, mass) for the Hamiltonian.

```python
theta0,p0=np.pi/2,0
state = PhaseSpaceState([theta,p],[theta0,p0])
pars=dict()
pars[g]=9.8
pars[l] = 1
pars[m] = 0.15
ham = Hamiltonian(H,pars,state)
```

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

### Set up Simulation for Resonance Analysis

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/Laplace-LagrangeSecularTheory.ipynb

Defines a helper function `get_sim` to create a simulation object with a specified planet-star mass ratio and orbital parameters. This function is used to set up systems for analyzing resonance effects.

```python
def get_sim(pratio):
    sim = rb.Simulation()
    sim.add(m=1)
    sim.add(m = 3e-5, P = 1 ,e = 0.01, inc = 0.01, pomega = 0)
    sim.add(m = 3e-5, P = pratio ,e = 0.01, pomega = np.pi/2)
    sim.integrator = 'whfast'
    sim.ri_whfast.safe_mode = 0 
    set_time_step(sim,1/27.)    
    sim.move_to_com()
    align_simulation(sim)
    return sim
```

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

### Generate Trajectories for the Standard Map

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/ChirikovStandardMap.ipynb

Generates multiple trajectories for the Standard Map. Each trajectory consists of a specified number of points, starting from different initial conditions.

```python
Ntraj = 20 # number of trajectories
Npts = 8000 # number of points per trajectory
all_points = np.zeros((Ntraj,Npts,2))

for i,p0 in enumerate(np.linspace(0,2*np.pi,Ntraj)):
    theta0 = np.random.uniform(0,2*np.pi)
    x = (theta0,p0)
    for j in range(Npts):
        all_points[i,j] = x
        x = smap(x)
```

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

### Initialize LaplaceLagrangeSystem and Calculate Solution

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/Laplace-LagrangeTheoryWithGR.ipynb

Retrieves the initial simulation state from the `Simulationarchive`, sets up a `LaplaceLagrangeSystem` object, and calculates the secular eccentricity solution without GR corrections.

```python
sim0 = sa[0] # get initial simulation
ll_sys = LaplaceLagrangeSystem.from_Simulation(sim0) # set up analytic Laplace-Lagrange solution
soln_no_GR  = ll_sys.secular_eccentricity_solution(nbody_soln['time']) # calculate analytic solution 

ll_sys.eccentricity_matrix # show the eccentricity matrix
```

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

### Symbolic Definition of Flow and Jacobian

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/AndoyerResonanceCapture.ipynb

Defines the symbolic flow vector and Jacobian matrix for the dynamical system using sympy. This is a setup step for numerical integration.

```python
x,y,delta,tau = sympy.symbols("x,y,delta,tau")
dyvars=[y,x,delta]
flow = sympy.Matrix([ham.flow[0],ham.flow[1],1/tau])
jac = sympy.Matrix(3,3,lambda i,j: sympy.diff(flow[i],dyvars[j]))
Nflow = sympy.lambdify(dyvars + [tau],flow)
Njac = sympy.lambdify(dyvars + [tau],jac)
Nham=sympy.lambdify(dyvars,ham.H)
Hsxfn = lambda delta:  Nham(0,xunst(delta),delta)
```

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

### Initialize Simulation and Laplace-Lagrange System

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/Laplace-LagrangeSecularTheory.ipynb

Initializes an N-body simulation (`sim`) with a specific orbital period ratio close to a 3:2 resonance, and then creates a `LaplaceLagrangeSystem` object from this simulation. It also calculates the maximum simulation time based on the system's eigenvalues.

```python
Delta0 = 0.02
sim = get_sim(kres * (1 + 0.02) / (kres - 1))
llsys = LaplaceLagrangeSystem.from_Simulation(sim)
Tmax = 4 * np.max( 2 * np.pi / np.abs(llsys.eccentricity_eigenvalues()) )
```

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

### Import necessary libraries for celmech

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/CanonicalTransformations.ipynb

Imports essential libraries including numpy, matplotlib, rebound, celmech components, and sympy for symbolic manipulation. Ensure these are installed before use.

```python
import numpy as np
import matplotlib.pyplot as plt
import rebound as rb
from celmech import PoincareHamiltonian,Poincare
from celmech.canonical_transformations import CanonicalTransformation

import sympy
from sympy import init_printing
init_printing()
```

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

### Plot Stability of Equilibria

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/OverstableLibrations.ipynb

Plots the stable and unstable equilibria based on the calculated maximum real part of eigenvalues. It also initializes a REBOUND simulation for an example unstable equilibrium.

```python
stable_mask = eig_max<0
plt.plot(e1[stable_mask],e2[stable_mask],color='k',label="Stable")
unstable_mask = eig_max>0
plt.plot(e1[unstable_mask],e2[unstable_mask],color='r',label="Unstable")

# Select and example equilibrium and initialize a simulation there
I = np.argmax(eig_max) + 100
inst_time = 1 / eig_max[I]
print(inst_time/res_eqs.tau_alpha)
Kexample = Kvals[I]
set_Kvals(res_eqs,Kexample)
dyvars  = eq_vars[I]
sim,extras = res_eqs.dyvars_to_rebound_simulation(dyvars,include_dissipation=True,osculating_correction=False)
els_eq = res_eqs.dyvars_to_orbital_elements(dyvars)
plt.scatter(els_eq['e1'],els_eq['e2'],color='gray',zorder=99,s=80,label="Simulation")

plt.legend(fontsize=18)
plt.tick_params(direction='in',labelsize=18,size=8)
plt.xlabel(r"$e_1$",fontsize=24)
plt.ylabel(r"$e_2$",fontsize=24)
```

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

### Integrating Initial Conditions with REBOUND and REBOUNDx

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/PlanarResonanceEquations.ipynb

Integrates initial conditions using REBOUND and REBOUNDx, including dissipation. This involves converting celmech variables to a REBOUND simulation and saving the state for later retrieval.

```python
sim,rebx = res_eqs.dyvars_to_rebound_simulation(dyvars,include_dissipation=True,osculating_correction=False)
rebx.save("rebxarchive.bin")
nb_results_by_planet = get_results_from_sim(
    res_eqs,
    sim,
    times,
    sa_filename='sa_dis.bin',
    extras_filename='rebxarchive.bin'
)
```

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

### Generate initial conditions from simulation

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/SpatialResoanceEquations.ipynb

Converts a rebound simulation into dynamical variables using `dyvars_from_rebound_simulation`. This sets up the initial state for the resonance equations.

```python
sim = get_sim(res_eqs,0.001,0.05,-0.05,np.pi/3,0.15,3*np.pi/2)
dyvars = res_eqs.dyvars_from_rebound_simulation(sim,iIn=1,iOut=2)
```

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

### Get Orbital Angles and Eccentricities

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/ResonantChain.ipynb

Integrates a REBOUND simulation over specified times and extracts orbital angles and eccentricities for each particle. Requires a REBOUND simulation object and a list of times.

```python
# Run a rebound simulation and get orbital angles and eccentricities
def get_angles_and_eccentricities(sim,times):
    N = len(times) 
    ps = sim.particles
    Npl = len(ps) - 1
    angs = np.zeros((N,3*Npl))
    eccs = np.zeros((N,Npl))
    for i,t in enumerate(times):
        sim.integrate(t)
        for j,p in enumerate(ps[1:]):
            eccs[i,j] = p.e
            angs[i,j] = p.l
            angs[i,j+Npl] = -1 * p.pomega
            angs[i,j+2*Npl] = -1 * p.Omega
    return angs,eccs
```

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

### Initialize Simulation and Predict Instability Times

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/Lammers24DynamicalModels.ipynb

Initializes a simulation with a given period ratio and predicts instability times by constructing a Hamiltonian model including first and second-order mean-motion resonances.

```python
P_ratios = []
log_t_insts = []

# predict instability times for 25 systems
for _ in tqdm(range(25)):
    P_ratio = np.random.uniform(1.03, 1.11)
    sim = initialize_sim(P_ratio)
    
    # get indices of bordering first-order MMRs
    j1 = int(np.floor(1/(1 - 1/P_ratio)))
    j2 = int(np.ceil(1/(1 - 1/P_ratio)))
    j3 = j1 + j2 #second-order MMR that lies between the first-order MMRs
    
    res_inds1 = [[j1, 1, i, i+1] for i in range(1, 5)] # all j1 resonances
    res_inds2 = [[j2, 1, i, i+1] for i in range(1, 5)] # all j2 resonances
    res_inds3 = [[j3, 2, i, i+1] for i in range(1, 5)] # all j3 resonances
    
    pham = construct_Hamiltonian(sim, res_inds1 + res_inds2 + res_inds3)
    t_inst = predict_t_inst(pham)
    
    P_ratios.append(P_ratio)
    log_t_insts.append(np.log10(t_inst))
```

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

### Plot N-body Solutions from Equilibria

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/ResonantChain.ipynb

Generates plots of eccentricities and orbital angles for N-body simulations starting from stable equilibria. It integrates simulations and visualizes the results using matplotlib.

```python
# Simulation times to get N-body results for
times = np.linspace(0,1,256) * 2e3 * sim.particles[1].P

# Set up plots
fig,ax = plt.subplots(2,6,figsize=(15,5),sharex=True,sharey='row')

# Loop over 6 stable equilibria, plot N-body solutions
for i,soln in enumerate(initial_solns):
    sim = to_sim(soln,amd0,sim)
    angs,eccs = get_angles_and_eccentricities(sim,times)
    ax[0,i].plot(times,eccs)
    res_angs = np.mod(np.array([A @ ang for ang in angs]),2*np.pi)
    ax[1,i].plot(times,180*res_angs[:,0]/np.pi)
    ax[1,i].plot(times,180*res_angs[:,1]/np.pi)
    
# make things  prettier
ax[0,0].set_ylabel(r"$e_i$ [deg]",fontsize = 16)
ax[1,0].set_ylabel(r"$\\phi_i$ [deg]",fontsize = 16)
ax[1,0].set_yticks([0,90,180,270,360])
```

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

### Initialize CometMap for 20:1 MMR

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/CometMap.ipynb

Initializes the CometMap class with parameters relevant to the 20:1 MMR with Neptune. Requires Neptune's semi-major axis, particle pericenter distance, Neptune's mass, and the resonance order.

```python
from celmech.maps import CometMap
import numpy as np
from matplotlib import pyplot as plt
```

```python
aNep = 30
q = 44.1
mNep = 5.15e-5
Nres = 20
cmap = CometMap(mNep,Nres,q/aNep,max_kmax=64)
```

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

### Initialize SecularSystemSimulation

Source: https://github.com/shadden/celmech/blob/master/docs/secular.md

Initializes a `SecularSystemSimulation` object from an existing REBOUND simulation (`sim0`). Specifies the integration method, order of expansion, time-step fraction, and Runge-Kutta specific keyword arguments.

```python
sec_sim = SecularSystemSimulation.from_Simulation(
    sim0,
    method='RK',
    max_order = 4,
    dtFraction=1/10.,
    rk_kwargs={'rtol':1e-10,'rk_method':'GL6'}
)
```

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

### Plotting Action-Angle Variables vs. Frequency

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/AndoyerResonanceCapture.ipynb

Plots the upper action J against the frequency omega for different delta values. Requires setup of action and frequency calculation functions.

```python
eigvals_ell = lambda d: np.linalg.eigvals(Njac_cons(yell(d),0,0,d)[:2,:2])
omega_J0 = lambda d: np.imag(eigvals_ell(d))[0]

fvals = np.concatenate((np.linspace(.3,0.8,5),1-np.logspace(-3.5,-1,5)))
N = len(fvals)

for delta in [1,1.5,2]:
    Jsx = sx_upper_action(delta) + sx_lower_action(delta)
    omega,J = np.zeros((2,N))
    for i,f in enumerate(fvals):
        J[i],omega[i] = get_action_and_frequency(delta,get_ic(delta,f))
    
    Js = np.concatenate(([Jsx],J,[0]))
    omegas = np.concatenate(([0],omega,[omega_J0(delta)]))
    order = np.argsort(Js)
    plt.plot(
        Js[order],
        omegas[order],
        's-',label=r"$\delta={:.1f}$".format(delta)
    )
    
    plt.scatter(Jsx,0)
    
plt.xlabel("$J$",fontsize=20)
plt.ylabel("$\\omega$",fontsize=20)
plt.tick_params(labelsize=16,direction='in',size=8)
plt.ylim(ymin=0)

plt.legend()
```

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

### Integrate N-body Simulation and Get Results

Source: https://github.com/shadden/celmech/blob/master/docs/secular.md

Integrates the REBOUND simulation over multiple secular timescales and retrieves the simulation results. This function is used to obtain data for comparing with secular solutions.

```python
Tsec = llsystem.Tsec # Shortest secular mode period of the system
sim.automateSimulationArchive("secular_example.sa",interval=Tsec/100,deletefile=True)
sim.integrate(10 * Tsec)
sa = rb.Simulationarchive("secular_example.sa")
nbody_results = get_simarchive_integration_results(sa,coordinates='heliocentric')
```

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

### Simulate and predict instability times for multiple systems (Model 1)

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/Lammers24DynamicalModels.ipynb

Iterates through 25 systems, initializing each with a random period ratio, constructing a Hamiltonian with first-order MMRs, and predicting its instability time. This loop generates data for plotting the relationship between initial period ratio and instability.

```python
P_ratios = []
log_t_insts = []

# predict instability times for 25 systems
for _ in tqdm(range(25)):
    P_ratio = np.random.uniform(1.03, 1.11)
    sim = initialize_sim(P_ratio)
    
    j = int(np.round(1/(1 - 1/P_ratio))) # get index of nearest first-order MMR
    res_inds = [[j, 1, i, i+1] for i in range(1, 5)] # all resonances to be included
    
    pham = construct_Hamiltonian(sim, res_inds)
    t_inst = predict_t_inst(pham)
    
    P_ratios.append(P_ratio)
    log_t_insts.append(np.log10(t_inst))
```

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

### Set Integrator and Prepare for Simulation

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/SimplePendulum.ipynb

Configures the numerical integrator ('vode' with 'adams' method) and calculates the initial period of oscillation to set the simulation time span. It also initializes an array to store the solution.

```python
ham.integrator.set_integrator('vode',method='adams')

omega0 = float(sqrt(g/l).subs(ham.H_params))
T0 = 2 * np.pi / omega0
times  = np.linspace(0,4*T0,100)
soln = np.zeros((len(times),len(state.values)))
```

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

### Integrate Hamiltonian System

Source: https://github.com/shadden/celmech/blob/master/docs/hamiltonian.md

Integrate the Hamiltonian system over a specified time range, updating the internal phase space state. The integration loop includes an example of modifying a Hamiltonian parameter mid-simulation.

```python
omega0 = float(sqrt(g/l).subs(ham.H_params))
T0 = 2 * np.pi / omega0
times  = np.linspace(0,4*T0,100)
soln = np.zeros((len(times),len(ham.N_dim)))

for i,t in enumerate(times):
    # Double the length of of the pendulum at the  50th step
    if i==50:
        ham.H_params[l] *= 2
    ham.integrate(t)
    soln[i]=ham.state.values
```

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

### Set up REBOUND Simulation for Kepler-223

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/ResonantChain.ipynb

Initializes a REBOUND simulation with observed masses and orbital periods for the Kepler-223 system. Units are set to Msun, days, and AU.

```python
# Taken from Table 1 of Mills+ (2016)
Mstar = 1.125
Mearth = 3e-6

periods = [7.38449,9.84564,14.78869,19.72567] # Periods in days
masses = [7.4,5.1,8.0,4.8] # Masses in Earth masses


sim = rb.Simulation()
sim.units  = ("Msun","days","AU")
sim.add(m=Mstar)
for mass,per in zip(masses,periods):
    sim.add(m = mass*Mearth,P =per)
sim.move_to_com()
```

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

### Calculating Reference Values for Momenta

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/ResonantChain.ipynb

Defines a list of periods and calculates ratios of Lambda0, which are then used to determine the reference values (Phi0) for the canonical momenta. This setup is crucial for defining the deviation variables.

```python
from sympy import S
periods = [0,3,4,6,8] # absolute value doesn't matter, just need to be in ratio of 3:4:6:8 for this chain
Lambda0Ratios = [(p/S(3))**(1/S(3)) * S("mu{}".format(i))/S("mu1") for i,p in enumerate(periods)]
Lambda0Ratios
```

```python
Phi0 = [Lambda0Ratios[1], 
        2*Lambda0Ratios[1] + Lambda0Ratios[2], 
        8*Lambda0Ratios[1] + 6*Lambda0Ratios[2] + 4*Lambda0Ratios[3] + 3*Lambda0Ratios[4],
        Lambda0Ratios[1] + Lambda0Ratios[2] + Lambda0Ratios[3] + Lambda0Ratios[4]]
Phi0
```

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

### Instantiate CometMap

Source: https://github.com/shadden/celmech/blob/master/jupyter_examples/ScatteredDiskResonanceOverlapBoundary.ipynb

Creates an instance of the CometMap class with specified parameters for perturber mass, initial comet period, and initial pericenter distance.

```python
mNeptune = 5.15e-5
aN = 30
cmap = CometMap(mNeptune,10,37/aN,max_kmax=64)
```