{ "cells": [ { "cell_type": "markdown", "id": "a1b2c3d4e5f6", "metadata": { "papermill": { "duration": 0.00534, "end_time": "2026-06-30T01:07:32.977088+00:00", "exception": false, "start_time": "2026-06-30T01:07:32.971748+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-06-30T01:07:32.986892Z", "iopub.status.busy": "2026-06-30T01:07:32.986651Z", "iopub.status.idle": "2026-06-30T01:07:35.277480Z", "shell.execute_reply": "2026-06-30T01:07:35.276558Z" }, "papermill": { "duration": 2.296382, "end_time": "2026-06-30T01:07:35.278214+00:00", "exception": false, "start_time": "2026-06-30T01:07:32.981832+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.003721, "end_time": "2026-06-30T01:07:35.286069+00:00", "exception": false, "start_time": "2026-06-30T01:07:35.282348+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-06-30T01:07:35.298178Z", "iopub.status.busy": "2026-06-30T01:07:35.297789Z", "iopub.status.idle": "2026-06-30T01:07:35.308789Z", "shell.execute_reply": "2026-06-30T01:07:35.308048Z" }, "papermill": { "duration": 0.019942, "end_time": "2026-06-30T01:07:35.309922+00:00", "exception": false, "start_time": "2026-06-30T01:07:35.289980+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.00376, "end_time": "2026-06-30T01:07:35.317854+00:00", "exception": false, "start_time": "2026-06-30T01:07:35.314094+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-06-30T01:07:35.326667Z", "iopub.status.busy": "2026-06-30T01:07:35.326449Z", "iopub.status.idle": "2026-06-30T01:07:35.576714Z", "shell.execute_reply": "2026-06-30T01:07:35.575765Z" }, "papermill": { "duration": 0.255725, "end_time": "2026-06-30T01:07:35.577433+00:00", "exception": false, "start_time": "2026-06-30T01:07:35.321708+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.003869, "end_time": "2026-06-30T01:07:35.585412+00:00", "exception": false, "start_time": "2026-06-30T01:07:35.581543+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-06-30T01:07:35.595035Z", "iopub.status.busy": "2026-06-30T01:07:35.594753Z", "iopub.status.idle": "2026-06-30T01:07:35.598462Z", "shell.execute_reply": "2026-06-30T01:07:35.597531Z" }, "papermill": { "duration": 0.009501, "end_time": "2026-06-30T01:07:35.599097+00:00", "exception": false, "start_time": "2026-06-30T01:07:35.589596+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.003704, "end_time": "2026-06-30T01:07:35.607056+00:00", "exception": false, "start_time": "2026-06-30T01:07:35.603352+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-06-30T01:07:35.616031Z", "iopub.status.busy": "2026-06-30T01:07:35.615808Z", "iopub.status.idle": "2026-06-30T01:07:35.622598Z", "shell.execute_reply": "2026-06-30T01:07:35.621897Z" }, "papermill": { "duration": 0.012705, "end_time": "2026-06-30T01:07:35.623711+00:00", "exception": false, "start_time": "2026-06-30T01:07:35.611006+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.004007, "end_time": "2026-06-30T01:07:35.632989+00:00", "exception": false, "start_time": "2026-06-30T01:07:35.628982+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-06-30T01:07:35.642127Z", "iopub.status.busy": "2026-06-30T01:07:35.641924Z", "iopub.status.idle": "2026-06-30T01:07:35.647727Z", "shell.execute_reply": "2026-06-30T01:07:35.647094Z" }, "papermill": { "duration": 0.011535, "end_time": "2026-06-30T01:07:35.648547+00:00", "exception": false, "start_time": "2026-06-30T01:07:35.637012+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.00381, "end_time": "2026-06-30T01:07:35.657813+00:00", "exception": false, "start_time": "2026-06-30T01:07:35.654003+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-06-30T01:07:35.666869Z", "iopub.status.busy": "2026-06-30T01:07:35.666671Z", "iopub.status.idle": "2026-06-30T01:07:35.704953Z", "shell.execute_reply": "2026-06-30T01:07:35.704468Z" }, "papermill": { "duration": 0.044197, "end_time": "2026-06-30T01:07:35.706086+00:00", "exception": false, "start_time": "2026-06-30T01:07:35.661889+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.003873, "end_time": "2026-06-30T01:07:35.714261+00:00", "exception": false, "start_time": "2026-06-30T01:07:35.710388+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-06-30T01:07:35.723326Z", "iopub.status.busy": "2026-06-30T01:07:35.723089Z", "iopub.status.idle": "2026-06-30T01:07:35.728744Z", "shell.execute_reply": "2026-06-30T01:07:35.728274Z" }, "papermill": { "duration": 0.012008, "end_time": "2026-06-30T01:07:35.730190+00:00", "exception": false, "start_time": "2026-06-30T01:07:35.718182+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.003861, "end_time": "2026-06-30T01:07:35.739052+00:00", "exception": false, "start_time": "2026-06-30T01:07:35.735191+00:00", "status": "completed" }, "tags": [] }, "source": [ "