{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e7f2fa69",
   "metadata": {
    "papermill": {
     "duration": 0.002338,
     "end_time": "2026-05-16T13:54:54.153253+00:00",
     "exception": false,
     "start_time": "2026-05-16T13:54:54.150915+00:00",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Using galpy Potentials in Other Codes\n",
    "\n",
    "galpy potentials can be exported to and used within several other stellar dynamics and N-body simulation frameworks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ae9d4a13",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-05-16T13:54:54.158070Z",
     "iopub.status.busy": "2026-05-16T13:54:54.157831Z",
     "iopub.status.idle": "2026-05-16T13:54:56.560647Z",
     "shell.execute_reply": "2026-05-16T13:54:56.559835Z"
    },
    "papermill": {
     "duration": 2.406813,
     "end_time": "2026-05-16T13:54:56.561964+00:00",
     "exception": false,
     "start_time": "2026-05-16T13:54:54.155151+00:00",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "\n",
    "warnings.filterwarnings(\"ignore\", category=FutureWarning)\n",
    "import numpy\n",
    "from galpy.potential import MiyamotoNagaiPotential, MWPotential2014"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6d0f1a4",
   "metadata": {
    "papermill": {
     "duration": 0.001543,
     "end_time": "2026-05-16T13:54:56.565824+00:00",
     "exception": false,
     "start_time": "2026-05-16T13:54:56.564281+00:00",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "## NEMO\n",
    "\n",
    "[NEMO](http://bima.astro.umd.edu/nemo/) is a toolbox for stellar dynamics. It provides programs for creating N-body realizations of galaxies, integrating their evolution, and analyzing the results. One of the most widely used programs within NEMO is gyrfalcON (Dehnen 2000, 2002), a fast, momentum-conserving N-body code.\n",
    "\n",
    "galpy potentials can be used as external gravitational fields in NEMO programs such as gyrfalcON. Two methods make this easy:\n",
    "\n",
    "- `nemo_accname()`: returns the NEMO acceleration name for the potential\n",
    "- `nemo_accpars(vo, ro)`: returns the NEMO acceleration parameters, where `vo` and `ro` are the velocity and distance scales in km/s and kpc\n",
    "\n",
    "For a single potential:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "17b84423",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-05-16T13:54:56.570301Z",
     "iopub.status.busy": "2026-05-16T13:54:56.569936Z",
     "iopub.status.idle": "2026-05-16T13:54:56.574032Z",
     "shell.execute_reply": "2026-05-16T13:54:56.573102Z"
    },
    "papermill": {
     "duration": 0.007192,
     "end_time": "2026-05-16T13:54:56.574717+00:00",
     "exception": false,
     "start_time": "2026-05-16T13:54:56.567525+00:00",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MiyamotoNagai\n"
     ]
    }
   ],
   "source": [
    "mp = MiyamotoNagaiPotential(a=0.5, b=0.0375, normalize=1.0)\n",
    "print(mp.nemo_accname())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "2957e76a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-05-16T13:54:56.579055Z",
     "iopub.status.busy": "2026-05-16T13:54:56.578874Z",
     "iopub.status.idle": "2026-05-16T13:54:56.582373Z",
     "shell.execute_reply": "2026-05-16T13:54:56.581543Z"
    },
    "papermill": {
     "duration": 0.006475,
     "end_time": "2026-05-16T13:54:56.583008+00:00",
     "exception": false,
     "start_time": "2026-05-16T13:54:56.576533+00:00",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0,592617.1107469305,4.0,0.3\n"
     ]
    }
   ],
   "source": [
    "print(mp.nemo_accpars(220.0, 8.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "823f6227",
   "metadata": {
    "papermill": {
     "duration": 0.001704,
     "end_time": "2026-05-16T13:54:56.586574+00:00",
     "exception": false,
     "start_time": "2026-05-16T13:54:56.584870+00:00",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "For a composite potential like `MWPotential2014`, these methods automatically combine the individual components:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4b090729",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-05-16T13:54:56.591121Z",
     "iopub.status.busy": "2026-05-16T13:54:56.590923Z",
     "iopub.status.idle": "2026-05-16T13:54:56.594387Z",
     "shell.execute_reply": "2026-05-16T13:54:56.593580Z"
    },
    "papermill": {
     "duration": 0.006826,
     "end_time": "2026-05-16T13:54:56.595304+00:00",
     "exception": false,
     "start_time": "2026-05-16T13:54:56.588478+00:00",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PowSphwCut+MiyamotoNagai+NFW\n"
     ]
    }
   ],
   "source": [
    "print(MWPotential2014.nemo_accname())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f4ccc5ec",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-05-16T13:54:56.599870Z",
     "iopub.status.busy": "2026-05-16T13:54:56.599686Z",
     "iopub.status.idle": "2026-05-16T13:54:56.603228Z",
     "shell.execute_reply": "2026-05-16T13:54:56.602349Z"
    },
    "papermill": {
     "duration": 0.006714,
     "end_time": "2026-05-16T13:54:56.603892+00:00",
     "exception": false,
     "start_time": "2026-05-16T13:54:56.597178+00:00",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0,1001.7912681044569,1.8,1.9#0,306770.4183858958,3.0,0.28#0,16.0,162.95824180864443\n"
     ]
    }
   ],
   "source": [
    "print(MWPotential2014.nemo_accpars(220.0, 8.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51379c23",
   "metadata": {
    "papermill": {
     "duration": 0.001836,
     "end_time": "2026-05-16T13:54:56.607692+00:00",
     "exception": false,
     "start_time": "2026-05-16T13:54:56.605856+00:00",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "You could then use these in a gyrfalcON command like:\n",
    "\n",
    "```bash\n",
    "gyrfalcON in.nemo out.nemo tstop=200 eps=0.05 step=0.01 kmax=8 \\\n",
    "  accname=MiyamotoNagai+NFW+MiyamotoNagai \\\n",
    "  accpars=<output of nemo_accpars>\n",
    "```\n",
    "\n",
    "galpy also includes a custom NEMO potential `PowSphwCut` in the `galpy/nemo/` directory, which can be compiled and used for potentials based on `PowerSphericalPotentialwCutoff`.\n",
    "\n",
    "The following galpy potentials support NEMO export:\n",
    "\n",
    "- `MiyamotoNagaiPotential`\n",
    "- `NFWPotential`\n",
    "- `PowerSphericalPotentialwCutoff`\n",
    "- `PlummerPotential`\n",
    "- `MN3ExponentialDiskPotential`\n",
    "- `LogarithmicHaloPotential`\n",
    "- Any list of supported potentials"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2492b1b",
   "metadata": {
    "papermill": {
     "duration": 0.001814,
     "end_time": "2026-05-16T13:54:56.611455+00:00",
     "exception": false,
     "start_time": "2026-05-16T13:54:56.609641+00:00",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "## AMUSE\n",
    "\n",
    "[AMUSE](http://amusecode.org/) (the Astrophysical Multipurpose Software Environment) is a Python framework for astrophysical simulations that combines existing codes for gravitational dynamics, stellar evolution, hydrodynamics, and radiative transfer.\n",
    "\n",
    "galpy potentials can be converted to an AMUSE-compatible representation using the `to_amuse` function:\n",
    "\n",
    "```python\n",
    "from galpy.potential import to_amuse, MWPotential2014\n",
    "mwp_amuse = to_amuse(MWPotential2014)\n",
    "```\n",
    "\n",
    "The `to_amuse` function accepts the following keyword arguments:\n",
    "\n",
    "- `ro=`, `vo=`: distance and velocity scales for unit conversion (default to galpy defaults)\n",
    "- `t=`: current time in AMUSE units (default: 0 Myr)\n",
    "- `tgalpy=`: current time in galpy natural units (alternative to `t=`)\n",
    "- `reverse=`: if `True`, reverse the sign of the potential (useful for certain AMUSE bridge configurations)\n",
    "\n",
    "Here is a complete example of setting up a Plummer-sphere cluster and evolving its N-body dynamics using an AMUSE ``BHTree`` in the external MWPotential2014 potential is:\n",
    "\n",
    "```python\n",
    "from amuse.lab import *\n",
    "from amuse.couple import bridge\n",
    "from amuse.datamodel import Particles\n",
    "from galpy.potential import to_amuse, MWPotential2014\n",
    "from galpy.util import plot as galpy_plot\n",
    "\n",
    "# Convert galpy MWPotential2014 to AMUSE representation\n",
    "mwp_amuse= to_amuse(MWPotential2014)\n",
    "\n",
    "# Set initial cluster parameters\n",
    "N= 1000\n",
    "Mcluster= 1000. | units.MSun\n",
    "Rcluster= 10. | units.parsec\n",
    "Rinit= [10.,0.,0.] | units.kpc\n",
    "Vinit= [0.,220.,0.] | units.km/units.s\n",
    "# Setup star cluster simulation\n",
    "tend= 100.0 | units.Myr\n",
    "dtout= 5.0 | units.Myr\n",
    "dt= 1.0 | units.Myr\n",
    "\n",
    "def setup_cluster(N,Mcluster,Rcluster,Rinit,Vinit):\n",
    "    converter= nbody_system.nbody_to_si(Mcluster,Rcluster)\n",
    "    stars= new_plummer_sphere(N,converter)\n",
    "    stars.x+= Rinit[0]\n",
    "    stars.y+= Rinit[1]\n",
    "    stars.z+= Rinit[2]\n",
    "    stars.vx+= Vinit[0]\n",
    "    stars.vy+= Vinit[1]\n",
    "    stars.vz+= Vinit[2]\n",
    "    return stars,converter\n",
    "\n",
    "# Setup cluster\n",
    "stars,converter= setup_cluster(N,Mcluster,Rcluster,Rinit,Vinit)\n",
    "cluster_code= BHTree(converter,number_of_workers=1) #Change number of workers depending no. of CPUs\n",
    "cluster_code.parameters.epsilon_squared= (3. | units.parsec)**2\n",
    "cluster_code.parameters.opening_angle= 0.6\n",
    "cluster_code.parameters.timestep= dt\n",
    "cluster_code.particles.add_particles(stars)\n",
    "\n",
    "# Setup channels between stars particle dataset and the cluster code\n",
    "channel_from_stars_to_cluster_code= stars.new_channel_to(cluster_code.particles,\n",
    "                                       attributes=[\"mass\", \"x\", \"y\", \"z\", \"vx\", \"vy\", \"vz\"])\n",
    "channel_from_cluster_code_to_stars= cluster_code.particles.new_channel_to(stars,\n",
    "                                       attributes=[\"mass\", \"x\", \"y\", \"z\", \"vx\", \"vy\", \"vz\"])\n",
    "\n",
    "# Setup gravity bridge\n",
    "gravity= bridge.Bridge(use_threading=False)\n",
    "# Stars in cluster_code depend on gravity from external potential mwp_amuse (i.e., MWPotential2014)\n",
    "gravity.add_system(cluster_code, (mwp_amuse,))\n",
    "# External potential mwp_amuse still needs to be added to system so it evolves with time\n",
    "gravity.add_system(mwp_amuse,)\n",
    "# Set how often to update external potential\n",
    "gravity.timestep= cluster_code.parameters.timestep/2.\n",
    "# Evolve\n",
    "time= 0.0 | tend.unit\n",
    "while time<tend:\n",
    "    gravity.evolve_model(time+dt)\n",
    "    # If you want to output or analyze the simulation, you need to copy\n",
    "    # stars from cluster_code\n",
    "    #channel_from_cluster_code_to_stars.copy()\n",
    "\n",
    "    # If you edited the stars particle set, for example to remove stars from the\n",
    "    # array because they have been kicked far from the cluster, you need to\n",
    "    # copy the array back to cluster_code:\n",
    "    #channel_from_stars_to_cluster_code.copy()\n",
    "\n",
    "    # Update time\n",
    "    time= gravity.model_time\n",
    "\n",
    "channel_from_cluster_code_to_stars.copy()\n",
    "gravity.stop()\n",
    "\n",
    "galpy_plot.plot(stars.x.value_in(units.kpc),stars.y.value_in(units.kpc),'.',\n",
    "                    xlabel=r'$X\\,(\\mathrm{kpc})$',ylabel=r'$Y\\,(\\mathrm{kpc})$')\n",
    "```\n",
    "\n",
    "This should give a plot like the following, which shows a cluster in the first stages of disruption:\n",
    "\n",
    "<img src=\"../../images/potential-amuse-example.png\" alt=\"AMUSE example\" width=\"75%\">\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d812a7a0",
   "metadata": {
    "papermill": {
     "duration": 0.001802,
     "end_time": "2026-05-16T13:54:56.615168+00:00",
     "exception": false,
     "start_time": "2026-05-16T13:54:56.613366+00:00",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "## AGAMA\n",
    "\n",
    "The [AGAMA](https://github.com/GalacticDynamics-Oxford/Agama) library (Action-based Galaxy Modelling Architecture) by [Vasiliev (2019)](https://ui.adsabs.harvard.edu/abs/2019MNRAS.482.1525V/abstract) supports importing galpy potentials. See the [AGAMA documentation](https://github.com/GalacticDynamics-Oxford/Agama) for details."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7404a268",
   "metadata": {
    "papermill": {
     "duration": 0.001866,
     "end_time": "2026-05-16T13:54:56.618907+00:00",
     "exception": false,
     "start_time": "2026-05-16T13:54:56.617041+00:00",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "## gala\n",
    "\n",
    "The [gala](https://gala.adrian.pw/) package by [Price-Whelan (2017)](https://ui.adsabs.harvard.edu/abs/2017JOSS....2..388P/abstract) supports converting galpy potentials. See the [gala documentation](https://gala.adrian.pw/en/latest/interop.html) for details on interoperability."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "py312",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.13"
  },
  "papermill": {
   "default_parameters": {},
   "duration": 3.965966,
   "end_time": "2026-05-16T13:54:57.237689+00:00",
   "environment_variables": {},
   "exception": null,
   "input_path": "doc/source/tutorials/potentials/other_codes.ipynb",
   "output_path": "doc/source/tutorials/potentials/other_codes.ipynb",
   "parameters": {},
   "start_time": "2026-05-16T13:54:53.271723+00:00",
   "version": "2.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}