Hamiltonian Class

class galport.Hamiltonian(Htype, ham=None, dJdt=None, dthetadt=None, t=None, **kwargs)

Bases: object

A class for working with Hamiltonian systems.

This class provides a unified interface for various 1D Hamiltonian models, including pendulum, Taylor series expansions, and resonant axisymmetric Hamiltonians. It supports time-dependent parameters and offers methods for evaluation, integration, and analysis of fixed points.

Parameters:

• Htype (str) –

  Type of the Hamiltonian. Must be one of:

  • pendulum : Simple pendulum Hamiltonian

  • taylor : Generalized Hamiltonian using Taylor series in J

  • sqrt_taylor : Generalized Hamiltonian using Taylor series in √J

  • axisym_res : Resonant axisymmetric Hamiltonian (2D)

  • my_fun : User-defined Hamiltonian functions

• ham (callable, optional) – User-defined Hamiltonian function when Htype=’my_fun’. Signature: ham(J, theta, **kwargs)

• dJdt (callable, optional) – User-defined dJ/dt function when Htype=’my_fun’. Signature: dJdt(J, theta, **kwargs)

• dthetadt (callable, optional) – User-defined dθ/dt function when Htype=’my_fun’. Signature: dthetadt(J, theta, **kwargs)

• t ((N,) numpy array, optional) – Time array for time-dependent parameters. If provided, parameters in **kwargs will be interpolated using cubic splines.

• **kwargs (dict) – Additional parameters passed to the Hamiltonian functions. For time-dependent case (if t is defined), each parameter should be an array of same length as t.

Hamiltonian Types

• pendulum — Simple pendulum Hamiltonian

  \[H = a(J-J_0)^2 + b\cos\theta\]

  Required parameters: coef=[a, b], J0

• taylor — Generalized Hamiltonian using Taylor series in J

  \[H = h_0 + \sum_i \left[h_i^c\cos(i\theta) + h_i^s\sin(i\theta)\right]\]

  where \(h_i = a_0 + a_1 p + a_2 p^2 + \ldots\) and \(p = (J-J_0)\)

  Parameters:

  • n : (M,) int array — Fourier modes to include.

  Positive integers → cosine terms \(h_i^c\cos(i\theta)\)

  Negative integers → sine terms \(h_i^s\sin(i\theta)\)

  Zero → constant term \(h_0\)

  • coef : (M, degree+1) array — Polynomial coefficients for each mode in n.

  Each row contains \([a_0, a_1, a_2, \ldots]\) for the corresponding mode.

  • J0 : float — Reference action for the Taylor expansion.

• sqrt_taylor — Generalized Hamiltonian using Taylor series in \(\sqrt{J}\)

  \[H = h_0 + \sum_i \left[h_i^c\cos(i\theta) + h_i^s\sin(i\theta)\right]\]

  where \(h_i = a_1 p + a_2 p^2 + a_3 p^3 + \\ldots\) and \(p = \sqrt{J}\)

  Parameters:

  • n : (M,) int array — Fourier modes to include (same format as taylor).

  • coef : (M, degree+1) array — Polynomial coefficients for each mode.

• axisym_res — Resonant axisymmetric Hamiltonian (2D)

  \[H = a_0 + a_R J_R + a_{R2} J_R^2 + a_z J_z + a_{z2} J_z^2 + a_{zR} J_z J_R + a_{zR}^{res} J_z J_R \cos(\theta_z - \theta_R)\]

  Parameters:

  coef=[a0, aR, aR2, az, az2, azR, azR_res]

• my_fun — User-defined Hamiltonian functions

  Parameters:

  ham (callable), dJdt (callable), dthetadt (callable)

dJdt(J, theta, t=None)

Find dJ/dt

Parameters:

• J ((N,) numpy array) – The array of actions

• theta ((N,) numpy array) – The array of angles

• t (float, optional) – time, by default None

Returns:

dJ/dt – The array of actions’ time derivative

Return type:

(N,) numpy array

derivative(J, theta, t=None)

Find time derivatives of action and angle dJ/dt = -dH/dθ, dθ/dt = dH/dJ

Parameters:

• J ((N,) numpy array) – The array of actions

• theta ((N,) numpy array) – The array of angles

• t (float, optional) – time, by default None

Returns:

• dJdt ((N,) numpy array) – The array of actions’ time derivative dJ/dt

• dthetadt ((N,) numpy array) – The array of angles’ time derivative dθ/dt

dthetadt(J, theta, t=None)

Find the value of the Hamiltonian

Parameters:

• J ((N,) numpy array) – The array of actions

• theta ((N,) numpy array) – The array of angles

• t (float, optional) – time, by default None

Returns:

dJdt – The array of actions’ time derivative dJ/dt

Return type:

(N,) numpy array

find_J(H, theta, J0=0.1, t=None)

Find an action over a known angle and Hamiltonian

Parameters:

• H (float) – the value of the Hamiltonian

• theta ((N,) numpy array) – The array of actions

• J0 (float, optional) – the initial action, by default 0.1

• t (float or None, optional) – time, by default None

find_theta(H, J, theta0=0.0, t=None)

Find an action over a known action and Hamiltonian

Parameters:

• H (float) – the value of the Hamiltonian

• J ((N,) numpy array) – The array of actions

• theta0 (float, optional) – the initial angle, by default 0.

• t (float or None, optional) – time, by default None

fix_points(J_range=[0.001, 0.5], Nstart=100, t=None, tol=0.01)

Parameters:

• J_range (list, optional) – range of actions, by default [0, 0.5]

• Nstart (int, optional) – number of initial random points, by default 100

• t (float, optional) – time, by default None

• tol (float, optional) – tolerance between roots, by default 10**-2

Returns:

fix_point – J, theta, type(stable, unstable), H, librating time

Return type:

dict

hamiltonian(J, theta, t=None)

Find the value of the Hamiltonian

Parameters:

• J ((N,) numpy array) – The array of actions

• theta ((N,) numpy array) – The array of angles

• t (float, optional) – time, by default None

Returns:

H – The array of the Hamiltonian’s values

Return type:

(N,) numpy array

integrate(J0, theta0, t0=0.0, Tint=100, Nint=10000, rtol=1e-10, atol=1e-12)

Integrate Hamiltonian trajectories

Parameters:

• J0 ((N,) or (N,2) numpy array) –

  Initial actions

  For 2D Hamiltonians (axisym_res) array of shape (N,2) containing [JR, Jz]

• theta0 ((N,) numpy array) – Initial angles

• t0 (float, optional) – Initial time. Default: 0.0

• Tint (int, optional) – Total integration time. Default: 100

• Nint (int, optional) – Number of time steps for output. Default: 10000

• rtol (float, optional) – Relative/absolute tolerance for solve_ivp. Default: 1e-10, 1e-12

• atol (float, optional) – Relative/absolute tolerance for solve_ivp. Default: 1e-10, 1e-12

Returns:

• t_eval ((M,) array) – Time points from t0 to t0+Tint.

• (J, theta) (tuple of arrays) – Evolution of actions and angles.

  • 1D: J.shape = (N, M), theta.shape = (N, M)

  • 2D: J.shape = (2, N, M) ([JR, Jz]), theta.shape = (N, M)

Notes

Uses scipy.integrate.solve_ivp with DOP853 method.

jacobian(J, theta, t=None, eps=1e-06)

Compute second derivatives (Jacobian elements) of the Hamiltonian using finite differences.

The returned values correspond to:

\[\frac{\partial^2 H}{\partial J^2} = \frac{d}{dJ}\left(\frac{d\theta}{dt}\right)\]

\[\frac{\partial^2 H}{\partial \theta^2} = -\frac{d}{d\theta}\left(\frac{dJ}{dt}\right)\]

\[\frac{\partial^2 H}{\partial J \partial \theta} = -\frac{d}{dJ}\left(\frac{dJ}{dt}\right)\]

These derivatives are useful for stability analysis and determining the nature of fixed points (elliptic/hyperbolic).

Parameters:

• J ((N,) numpy array) – Action values.

• theta ((N,) numpy array) – Angle values.

• t (float, optional) – Time for time-dependent Hamiltonians.

• eps (float, optional) – Step size for finite difference approximation. Default: 1e-6

Returns:

• d2HdJ2 ((N,) numpy array) – Second derivative with respect to action.

• d2Hdtheta2 ((N,) numpy array) – Second derivative with respect to angle.

• d2HdJdtheta ((N,) numpy array) – Mixed second derivative.

Notes

Uses central finite differences: \(\dfrac{f(x+\epsilon) - f(x-\epsilon)}{2\epsilon}\)