{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e7f2fa69",
   "metadata": {
    "papermill": {
     "duration": 0.002152,
     "end_time": "2026-06-30T01:21:28.111630+00:00",
     "exception": false,
     "start_time": "2026-06-30T01:21:28.109478+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-06-30T01:21:28.116605Z",
     "iopub.status.busy": "2026-06-30T01:21:28.116358Z",
     "iopub.status.idle": "2026-06-30T01:21:30.340248Z",
     "shell.execute_reply": "2026-06-30T01:21:30.339357Z"
    },
    "papermill": {
     "duration": 2.227383,
     "end_time": "2026-06-30T01:21:30.341016+00:00",
     "exception": false,
     "start_time": "2026-06-30T01:21:28.113633+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.001455,
     "end_time": "2026-06-30T01:21:30.344607+00:00",
     "exception": false,
     "start_time": "2026-06-30T01:21:30.343152+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-06-30T01:21:30.348993Z",
     "iopub.status.busy": "2026-06-30T01:21:30.348371Z",
     "iopub.status.idle": "2026-06-30T01:21:30.352342Z",
     "shell.execute_reply": "2026-06-30T01:21:30.351627Z"
    },
    "papermill": {
     "duration": 0.006945,
     "end_time": "2026-06-30T01:21:30.353041+00:00",
     "exception": false,
     "start_time": "2026-06-30T01:21:30.346096+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-06-30T01:21:30.357144Z",
     "iopub.status.busy": "2026-06-30T01:21:30.356965Z",
     "iopub.status.idle": "2026-06-30T01:21:30.360157Z",
     "shell.execute_reply": "2026-06-30T01:21:30.359441Z"
    },
    "papermill": {
     "duration": 0.006039,
     "end_time": "2026-06-30T01:21:30.360760+00:00",
     "exception": false,
     "start_time": "2026-06-30T01:21:30.354721+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.001458,
     "end_time": "2026-06-30T01:21:30.363901+00:00",
     "exception": false,
     "start_time": "2026-06-30T01:21:30.362443+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-06-30T01:21:30.367920Z",
     "iopub.status.busy": "2026-06-30T01:21:30.367738Z",
     "iopub.status.idle": "2026-06-30T01:21:30.371122Z",
     "shell.execute_reply": "2026-06-30T01:21:30.370335Z"
    },
    "papermill": {
     "duration": 0.006224,
     "end_time": "2026-06-30T01:21:30.371688+00:00",
     "exception": false,
     "start_time": "2026-06-30T01:21:30.365464+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-06-30T01:21:30.375701Z",
     "iopub.status.busy": "2026-06-30T01:21:30.375517Z",
     "iopub.status.idle": "2026-06-30T01:21:30.378657Z",
     "shell.execute_reply": "2026-06-30T01:21:30.377997Z"
    },
    "papermill": {
     "duration": 0.006175,
     "end_time": "2026-06-30T01:21:30.379509+00:00",
     "exception": false,
     "start_time": "2026-06-30T01:21:30.373334+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.001564,
     "end_time": "2026-06-30T01:21:30.382761+00:00",
     "exception": false,
     "start_time": "2026-06-30T01:21:30.381197+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.001505,
     "end_time": "2026-06-30T01:21:30.386025+00:00",
     "exception": false,
     "start_time": "2026-06-30T01:21:30.384520+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.001578,
     "end_time": "2026-06-30T01:21:30.389280+00:00",
     "exception": false,
     "start_time": "2026-06-30T01:21:30.387702+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.001488,
     "end_time": "2026-06-30T01:21:30.392369+00:00",
     "exception": false,
     "start_time": "2026-06-30T01:21:30.390881+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.14"
  },
  "papermill": {
   "default_parameters": {},
   "duration": 3.750726,
   "end_time": "2026-06-30T01:21:31.010737+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-06-30T01:21:27.260011+00:00",
   "version": "2.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}