{ "cells": [ { "cell_type": "markdown", "id": "a1b2c3d4e5f6", "metadata": { "papermill": { "duration": 0.004996, "end_time": "2026-05-16T13:40:41.032029+00:00", "exception": false, "start_time": "2026-05-16T13:40:41.027033+00:00", "status": "completed" }, "tags": [] }, "source": [ "# Action-Angle Coordinates: Staeckel Approximation\n", "\n", "The Staeckel approximation (`actionAngleStaeckel`) is the most accurate general\n", "method for computing actions in axisymmetric potentials. It works by locally\n", "approximating the potential as a Staeckel potential (separable in confocal\n", "ellipsoidal coordinates), following the algorithm of\n", "[Binney (2012)](http://adsabs.harvard.edu/abs/2012MNRAS.426.1324B).\n", "\n", "For all intents and purposes, the Staeckel approximation makes the\n", "[adiabatic approximation](adiabatic.ipynb) obsolete: it is as fast and more precise. The only\n", "additional requirement is that the user must specify a *focal length* $\\Delta$, but this\n", "can be easily estimated from the potential.\n", "\n", "For an introduction to action-angle coordinates and simpler methods, see\n", "the [Introduction](introduction.ipynb). For the orbit-integration-based method\n", "that works for highly eccentric orbits, see [IsochroneApprox](isochroneapprox.ipynb)." ] }, { "cell_type": "code", "execution_count": 1, "id": "b2c3d4e5f6a7", "metadata": { "execution": { "iopub.execute_input": "2026-05-16T13:40:41.041676Z", "iopub.status.busy": "2026-05-16T13:40:41.041430Z", "iopub.status.idle": "2026-05-16T13:40:43.425310Z", "shell.execute_reply": "2026-05-16T13:40:43.424141Z" }, "papermill": { "duration": 2.389793, "end_time": "2026-05-16T13:40:43.426164+00:00", "exception": false, "start_time": "2026-05-16T13:40:41.036371+00:00", "status": "completed" }, "tags": [] }, "outputs": [], "source": [ "%matplotlib inline\n", "import numpy\n", "import matplotlib.pyplot as plt\n", "from galpy.potential import MWPotential2014\n", "from galpy.orbit import Orbit\n", "import warnings\n", "\n", "warnings.filterwarnings(\"ignore\", category=RuntimeWarning)\n", "warnings.filterwarnings(\"ignore\", category=UserWarning)" ] }, { "cell_type": "markdown", "id": "c3d4e5f6a7b8", "metadata": { "papermill": { "duration": 0.004187, "end_time": "2026-05-16T13:40:43.435054+00:00", "exception": false, "start_time": "2026-05-16T13:40:43.430867+00:00", "status": "completed" }, "tags": [] }, "source": [ "## Setup\n", "\n", "The key parameter is `delta`, the focus of the confocal coordinate system.\n", "You can estimate an appropriate value using `estimateDeltaStaeckel`, which\n", "estimates $\\Delta$ from the second derivatives of the potential\n", "(see [Sanders 2012](http://adsabs.harvard.edu/abs/2012MNRAS.426..128S)):" ] }, { "cell_type": "code", "execution_count": 2, "id": "d4e5f6a7b8c9", "metadata": { "execution": { "iopub.execute_input": "2026-05-16T13:40:43.445144Z", "iopub.status.busy": "2026-05-16T13:40:43.444719Z", "iopub.status.idle": "2026-05-16T13:40:43.455933Z", "shell.execute_reply": "2026-05-16T13:40:43.455140Z" }, "papermill": { "duration": 0.017386, "end_time": "2026-05-16T13:40:43.456935+00:00", "exception": false, "start_time": "2026-05-16T13:40:43.439549+00:00", "status": "completed" }, "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Estimated delta: 0.7013326670567281\n" ] } ], "source": [ "from galpy.actionAngle import actionAngleStaeckel, estimateDeltaStaeckel\n", "\n", "# Estimate delta for an orbit near R=1, z=0\n", "delta = estimateDeltaStaeckel(MWPotential2014, 1.0, 0.0)\n", "print(\"Estimated delta:\", delta)" ] }, { "cell_type": "markdown", "id": "e5f6a7b8c9d0", "metadata": { "papermill": { "duration": 0.004303, "end_time": "2026-05-16T13:40:43.465790+00:00", "exception": false, "start_time": "2026-05-16T13:40:43.461487+00:00", "status": "completed" }, "tags": [] }, "source": [ "For more accuracy, you can estimate `delta` along an integrated orbit by\n", "averaging (through the median) estimates at positions around the orbit:" ] }, { "cell_type": "code", "execution_count": 3, "id": "f6a7b8c9d0e1", "metadata": { "execution": { "iopub.execute_input": "2026-05-16T13:40:43.475731Z", "iopub.status.busy": "2026-05-16T13:40:43.475494Z", "iopub.status.idle": "2026-05-16T13:40:43.757963Z", "shell.execute_reply": "2026-05-16T13:40:43.757047Z" }, "papermill": { "duration": 0.288548, "end_time": "2026-05-16T13:40:43.758774+00:00", "exception": false, "start_time": "2026-05-16T13:40:43.470226+00:00", "status": "completed" }, "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Delta from orbit: 0.40272708628053305\n" ] } ], "source": [ "o = Orbit([1.0, 0.1, 1.1, 0.0, 0.25, 0.0])\n", "ts = numpy.linspace(0, 100, 1001)\n", "o.integrate(ts, MWPotential2014)\n", "\n", "# Estimate delta from the orbit's R and z range\n", "delta_orbit = estimateDeltaStaeckel(MWPotential2014, o.R(ts), o.z(ts))\n", "print(\"Delta from orbit:\", delta_orbit)" ] }, { "cell_type": "markdown", "id": "a7b8c9d0e1f2", "metadata": { "papermill": { "duration": 0.004232, "end_time": "2026-05-16T13:40:43.767642+00:00", "exception": false, "start_time": "2026-05-16T13:40:43.763410+00:00", "status": "completed" }, "tags": [] }, "source": [ "We will use $\\Delta = 0.4$ in what follows. We set up the\n", "`actionAngleStaeckel` object with `c=True` to use the fast C implementation:" ] }, { "cell_type": "code", "execution_count": 4, "id": "b8c9d0e1f2a3", "metadata": { "execution": { "iopub.execute_input": "2026-05-16T13:40:43.779781Z", "iopub.status.busy": "2026-05-16T13:40:43.779579Z", "iopub.status.idle": "2026-05-16T13:40:43.782979Z", "shell.execute_reply": "2026-05-16T13:40:43.782103Z" }, "papermill": { "duration": 0.009289, "end_time": "2026-05-16T13:40:43.783777+00:00", "exception": false, "start_time": "2026-05-16T13:40:43.774488+00:00", "status": "completed" }, "tags": [] }, "outputs": [], "source": [ "aAS = actionAngleStaeckel(pot=MWPotential2014, delta=0.4, c=True)" ] }, { "cell_type": "markdown", "id": "c9d0e1f2a3b4", "metadata": { "papermill": { "duration": 0.00429, "end_time": "2026-05-16T13:40:43.792499+00:00", "exception": false, "start_time": "2026-05-16T13:40:43.788209+00:00", "status": "completed" }, "tags": [] }, "source": [ "## Computing actions\n", "\n", "Calling the object returns $(J_R, L_z, J_z)$:" ] }, { "cell_type": "code", "execution_count": 5, "id": "d0e1f2a3b4c5", "metadata": { "execution": { "iopub.execute_input": "2026-05-16T13:40:43.802200Z", "iopub.status.busy": "2026-05-16T13:40:43.801964Z", "iopub.status.idle": "2026-05-16T13:40:43.811304Z", "shell.execute_reply": "2026-05-16T13:40:43.809868Z" }, "papermill": { "duration": 0.01513, "end_time": "2026-05-16T13:40:43.812007+00:00", "exception": false, "start_time": "2026-05-16T13:40:43.796877+00:00", "status": "completed" }, "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "J_R = 0.013636, L_z = 1.100000, J_z = 0.000463\n" ] } ], "source": [ "jr, jphi, jz = aAS(1.0, 0.1, 1.1, 0.0, 0.05)\n", "print(f\"J_R = {jr.item():.6f}, L_z = {jphi.item():.6f}, J_z = {jz.item():.6f}\")" ] }, { "cell_type": "markdown", "id": "e1f2a3b4c5d6", "metadata": { "papermill": { "duration": 0.004225, "end_time": "2026-05-16T13:40:43.822556+00:00", "exception": false, "start_time": "2026-05-16T13:40:43.818331+00:00", "status": "completed" }, "tags": [] }, "source": [ "For the orbit with higher vertical excursion:" ] }, { "cell_type": "code", "execution_count": 6, "id": "f2a3b4c5d6e7", "metadata": { "execution": { "iopub.execute_input": "2026-05-16T13:40:43.832140Z", "iopub.status.busy": "2026-05-16T13:40:43.831940Z", "iopub.status.idle": "2026-05-16T13:40:43.838075Z", "shell.execute_reply": "2026-05-16T13:40:43.836950Z" }, "papermill": { "duration": 0.011905, "end_time": "2026-05-16T13:40:43.838806+00:00", "exception": false, "start_time": "2026-05-16T13:40:43.826901+00:00", "status": "completed" }, "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "J_R = 0.019222, L_z = 1.100000, J_z = 0.015277\n" ] } ], "source": [ "jr_hz, jphi_hz, jz_hz = aAS(o.R(), o.vR(), o.vT(), o.z(), o.vz())\n", "print(f\"J_R = {jr_hz.item():.6f}, L_z = {jphi_hz.item():.6f}, J_z = {jz_hz.item():.6f}\")" ] }, { "cell_type": "markdown", "id": "a3b4c5d6e7f8", "metadata": { "papermill": { "duration": 0.004286, "end_time": "2026-05-16T13:40:43.848853+00:00", "exception": false, "start_time": "2026-05-16T13:40:43.844567+00:00", "status": "completed" }, "tags": [] }, "source": [ "## Computing actions, frequencies, and angles\n", "\n", "The `actionsFreqs` method returns $(J_R, L_z, J_z, \\Omega_R, \\Omega_\\phi, \\Omega_z)$:" ] }, { "cell_type": "code", "execution_count": 7, "id": "b4c5d6e7f8a9", "metadata": { "execution": { "iopub.execute_input": "2026-05-16T13:40:43.858687Z", "iopub.status.busy": "2026-05-16T13:40:43.858505Z", "iopub.status.idle": "2026-05-16T13:40:43.864571Z", "shell.execute_reply": "2026-05-16T13:40:43.864052Z" }, "papermill": { "duration": 0.01252, "end_time": "2026-05-16T13:40:43.865947+00:00", "exception": false, "start_time": "2026-05-16T13:40:43.853427+00:00", "status": "completed" }, "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "J_R = 0.013636, L_z = 1.100000, J_z = 0.000463\n", "Omega_R = 1.1595, Omega_phi = 0.8682, Omega_z = 2.1432\n" ] } ], "source": [ "jr, jphi, jz, Or, Op, Oz = aAS.actionsFreqs(1.0, 0.1, 1.1, 0.0, 0.05)\n", "print(f\"J_R = {jr.item():.6f}, L_z = {jphi.item():.6f}, J_z = {jz.item():.6f}\")\n", "print(\n", " f\"Omega_R = {Or.item():.4f}, Omega_phi = {Op.item():.4f}, Omega_z = {Oz.item():.4f}\"\n", ")" ] }, { "cell_type": "markdown", "id": "c5d6e7f8a9b0", "metadata": { "papermill": { "duration": 0.004319, "end_time": "2026-05-16T13:40:43.875927+00:00", "exception": false, "start_time": "2026-05-16T13:40:43.871608+00:00", "status": "completed" }, "tags": [] }, "source": [ "The `actionsFreqsAngles` method additionally returns the angles\n", "$(\\theta_R, \\theta_\\phi, \\theta_z)$. Note that you **must** specify `phi`\n", "for the angles:" ] }, { "cell_type": "code", "execution_count": 8, "id": "d6e7f8a9b0c1", "metadata": { "execution": { "iopub.execute_input": "2026-05-16T13:40:43.885730Z", "iopub.status.busy": "2026-05-16T13:40:43.885559Z", "iopub.status.idle": "2026-05-16T13:40:43.892642Z", "shell.execute_reply": "2026-05-16T13:40:43.890832Z" }, "papermill": { "duration": 0.013492, "end_time": "2026-05-16T13:40:43.893892+00:00", "exception": false, "start_time": "2026-05-16T13:40:43.880400+00:00", "status": "completed" }, "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "J_R = 0.013636\n", "Omega_R = 1.1595, Omega_phi = 0.8682, Omega_z = 2.1432\n", "theta_R = 0.4687, theta_phi = 6.1767, theta_z = 6.0535\n" ] } ], "source": [ "jr, jphi, jz, Or, Op, Oz, wr, wp, wz = aAS.actionsFreqsAngles(\n", " 1.0, 0.1, 1.1, 0.0, 0.05, 0.0\n", ")\n", "print(f\"J_R = {jr.item():.6f}\")\n", "print(\n", " f\"Omega_R = {Or.item():.4f}, Omega_phi = {Op.item():.4f}, Omega_z = {Oz.item():.4f}\"\n", ")\n", "print(\n", " f\"theta_R = {wr.item():.4f}, theta_phi = {wp.item():.4f}, theta_z = {wz.item():.4f}\"\n", ")" ] }, { "cell_type": "markdown", "id": "e7f8a9b0c1d2", "metadata": { "papermill": { "duration": 0.004319, "end_time": "2026-05-16T13:40:43.903507+00:00", "exception": false, "start_time": "2026-05-16T13:40:43.899188+00:00", "status": "completed" }, "tags": [] }, "source": [ "