diff --git a/requirements.txt b/requirements.txt index 050e7d4..0a5c013 100644 --- a/requirements.txt +++ b/requirements.txt @@ -41,4 +41,5 @@ nbclient torch SQLAlchemy psycopg2-binary -ipykernel \ No newline at end of file +ipykernel +astropy diff --git a/workspaces/galactocentric-frame/galactocentric-frame.ipynb b/workspaces/galactocentric-frame/galactocentric-frame.ipynb new file mode 100644 index 0000000..99bf156 --- /dev/null +++ b/workspaces/galactocentric-frame/galactocentric-frame.ipynb @@ -0,0 +1,335 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "487b2dde", + "metadata": {}, + "source": [ + "# Transforming positions and velocities to and from a Galactocentric frame" + ] + }, + { + "cell_type": "markdown", + "id": "69e7f209", + "metadata": {}, + "source": [ + "Source: https://docs.astropy.org/en/latest/coordinates/example_gallery_plot_galactocentric_frame.html" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "78b091ce", + "metadata": {}, + "outputs": [], + "source": [ + "def plot_positions_and_velocities(gc_rings):\n", + " fig, axes = plt.subplots(1, 2, figsize=(12, 6), dpi=200)\n", + " axes[0].plot(gc_rings.x.T, gc_rings.y.T, marker=\"None\", linewidth=3)\n", + " axes[0].text(-8.0, 0, r\"$\\odot$\", fontsize=20)\n", + " axes[0].set_xlim(-30, 30)\n", + " axes[0].set_ylim(-30, 30)\n", + " axes[0].set_xlabel(\"$x$ [kpc]\")\n", + " axes[0].set_ylabel(\"$y$ [kpc]\")\n", + " axes[0].set_title(\"Positions\")\n", + " axes[1].plot(gc_rings.v_x.T, gc_rings.v_y.T, marker=\"None\", linewidth=3)\n", + " axes[1].set_xlim(-250, 250)\n", + " axes[1].set_ylim(-250, 250)\n", + " axes[1].set_xlabel(f\"$v_x$ [{(u.km / u.s).to_string('latex_inline')}]\")\n", + " axes[1].set_ylabel(f\"$v_y$ [{(u.km / u.s).to_string('latex_inline')}]\")\n", + " axes[1].set_title(\"Velocities\")\n", + " fig.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "id": "83278590", + "metadata": {}, + "source": [ + "This example shows a few examples of how to use and customize the\n", + "`Galactocentric` frame to transform Heliocentric sky\n", + "positions, distance, proper motions, and radial velocities to a Galactocentric,\n", + "Cartesian frame, and the same in reverse.\n", + "\n", + "The main configurable parameters of the `Galactocentric`\n", + "frame control the position and velocity of the solar system barycenter within\n", + "the Galaxy. These are specified by setting the ICRS coordinates of the\n", + "Galactic center, the distance to the Galactic center (the sun-galactic center\n", + "line is always assumed to be the x-axis of the Galactocentric frame), and the\n", + "Cartesian 3-velocity of the sun in the Galactocentric frame. We will first\n", + "demonstrate how to customize these values, then show how to set the solar motion\n", + "instead by inputting the proper motion of Sgr A*.\n", + "\n", + "Note that, for brevity, we may refer to the solar system barycenter as just \"the\n", + "sun\" in the examples below.\n", + "\n", + "Let's first define a barycentric coordinate and velocity in the ICRS frame.\n", + "We will use the data for the star HD 39881 from the\n", + "[Simbad](https://simbad.unistra.fr/simbad/) database:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fabf0fa4", + "metadata": {}, + "outputs": [], + "source": [ + "import astropy.coordinates as coord\n", + "from astropy import units as u\n", + "c1 = coord.SkyCoord(\n", + " ra=89.014303 * u.degree,\n", + " dec=13.924912 * u.degree,\n", + " distance=(37.59 * u.mas).to(u.pc, u.parallax()),\n", + " pm_ra_cosdec=372.72 * (u.mas / u.yr),\n", + " pm_dec=-483.69 * (u.mas / u.yr),\n", + " radial_velocity=0.37 * (u.km / u.s),\n", + " frame=\"icrs\",\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "58097f73", + "metadata": {}, + "source": [ + "This is a high proper-motion star; suppose we'd like to transform its position\n", + "and velocity to a Galactocentric frame to see if it has a large 3D velocity\n", + "as well. To use the Astropy default solar position and motion parameters, we\n", + "can do the following:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b2af242d", + "metadata": {}, + "outputs": [], + "source": [ + "gc1 = c1.transform_to(coord.Galactocentric)" + ] + }, + { + "cell_type": "markdown", + "id": "e8ddeafe", + "metadata": {}, + "source": [ + "From here, we can access the components of the resulting\n", + "`Galactocentric` instance to see the 3D Cartesian\n", + "velocity components:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dd0b6eba", + "metadata": {}, + "outputs": [], + "source": [ + "print(gc1.v_x, gc1.v_y, gc1.v_z)\n" + ] + }, + { + "cell_type": "markdown", + "id": "7a2aae39", + "metadata": {}, + "source": [ + "The default parameters for the `Galactocentric` frame\n", + "are detailed in the linked documentation, but we can modify the most commonly\n", + "changed values using the keywords ``galcen_distance``, ``galcen_v_sun``, and\n", + "``z_sun`` which set the sun-Galactic center distance, the 3D velocity vector\n", + "of the sun, and the height of the sun above the Galactic midplane,\n", + "respectively. The velocity of the sun can be specified as an\n", + "`Quantity` object with velocity units and is interpreted as a\n", + "Cartesian velocity, as in the example below. Note that, as with the positions,\n", + "the Galactocentric frame is a right-handed system (i.e., the Sun is at negative\n", + "x values) so ``v_x`` is opposite of the Galactocentric radial velocity:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2d7bef03", + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "v_sun = [11.1, 244, 7.25] * (u.km / u.s) # [vx, vy, vz]\n", + "gc_frame = coord.Galactocentric(\n", + " galcen_distance=8 * u.kpc, galcen_v_sun=v_sun, z_sun=0 * u.pc\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "08791399", + "metadata": {}, + "source": [ + "We can then transform to this frame instead, with our custom parameters:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d8717a8b", + "metadata": {}, + "outputs": [], + "source": [ + "gc2 = c1.transform_to(gc_frame)\n", + "print(gc2.v_x, gc2.v_y, gc2.v_z)" + ] + }, + { + "cell_type": "markdown", + "id": "c0b353c6", + "metadata": {}, + "source": [ + "It is sometimes useful to specify the solar motion using the\n", + "[proper motion of Sgr A*](https://arxiv.org/abs/astro-ph/0408107)\n", + "instead of Cartesian velocity components. With an assumed distance, we can convert\n", + "proper motion components to Cartesian velocity components using `units`:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bccb4ce4", + "metadata": {}, + "outputs": [], + "source": [ + "galcen_distance = 8 * u.kpc\n", + "pm_gal_sgrA = [-6.379, -0.202] * (u.mas / u.yr) # from Reid & Brunthaler 2004\n", + "vy, vz = -(galcen_distance * pm_gal_sgrA).to(u.km / u.s, u.dimensionless_angles())" + ] + }, + { + "cell_type": "markdown", + "id": "b8853799", + "metadata": {}, + "source": [ + "We still have to assume a line-of-sight velocity for the Galactic center,\n", + "which we will again take to be 11 km/s:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cb4ac416", + "metadata": {}, + "outputs": [], + "source": [ + "vx = 11.1 * (u.km / u.s)\n", + "v_sun2 = u.Quantity([vx, vy, vz]) # List of Quantity -> a single Quantity\n", + "gc_frame2 = coord.Galactocentric(\n", + " galcen_distance=galcen_distance, galcen_v_sun=v_sun2, z_sun=0 * u.pc\n", + ")\n", + "gc3 = c1.transform_to(gc_frame2)\n", + "print(gc3.v_x, gc3.v_y, gc3.v_z)" + ] + }, + { + "cell_type": "markdown", + "id": "4e42feb4", + "metadata": {}, + "source": [ + "The transformations also work in the opposite direction. This can be useful\n", + "for transforming simulated or theoretical data to observable quantities. As\n", + "an example, we will generate 4 theoretical circular orbits at different\n", + "Galactocentric radii with the same circular velocity, and transform them to\n", + "Heliocentric coordinates:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a52595f5", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import astropy.coordinates as coord\n", + "from astropy import units as u\n", + "ring_distances = np.arange(10, 26, 5) * u.kpc\n", + "circ_velocity = 220 * (u.km / u.s)\n", + "phi_grid = np.linspace(90, 270, 512) * u.degree # grid of azimuths\n", + "ring_rep = coord.CylindricalRepresentation(\n", + " rho=ring_distances[:, np.newaxis],\n", + " phi=phi_grid[np.newaxis],\n", + " z=np.zeros_like(ring_distances)[:, np.newaxis],\n", + ")\n", + "angular_velocity = (-circ_velocity / ring_distances).to(\n", + " u.mas / u.yr, u.dimensionless_angles()\n", + ")\n", + "ring_dif = coord.CylindricalDifferential(\n", + " d_rho=np.zeros(phi_grid.shape)[np.newaxis] * (u.km / u.s),\n", + " d_phi=angular_velocity[:, np.newaxis],\n", + " d_z=np.zeros(phi_grid.shape)[np.newaxis] * (u.km / u.s),\n", + ")\n", + "ring_rep = ring_rep.with_differentials(ring_dif)\n", + "gc_rings = coord.SkyCoord(ring_rep, frame=coord.Galactocentric)" + ] + }, + { + "cell_type": "markdown", + "id": "de802b0a", + "metadata": {}, + "source": [ + "First, let's visualize the geometry in Galactocentric coordinates. Here are\n", + "the positions and velocities of the rings; note that in the velocity plot,\n", + "the velocities of the 4 rings are identical and thus overlaid under the same\n", + "curve:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b9795d77", + "metadata": {}, + "outputs": [], + "source": [ + "plot_positions_and_velocities(gc_rings)" + ] + }, + { + "cell_type": "markdown", + "id": "c83d570e", + "metadata": {}, + "source": [ + "Now we can transform to Galactic coordinates and visualize the rings in\n", + "observable coordinates:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e570027d", + "metadata": {}, + "outputs": [], + "source": [ + "gal_rings = gc_rings.transform_to(coord.Galactic)\n", + "fig, ax = plt.subplots(1, 1, figsize=(8, 6), dpi=200)\n", + "for i in range(len(ring_distances)):\n", + " ax.plot(\n", + " gal_rings[i].l.degree,\n", + " gal_rings[i].pm_l_cosb.value,\n", + " label=str(ring_distances[i]),\n", + " marker=\"None\",\n", + " linewidth=3,\n", + " )\n", + "ax.set_xlim(360, 0)\n", + "ax.set_xlabel(\"$l$ [deg]\")\n", + "ax.set_ylabel(rf'$\\mu_l \\, \\cos b$ [{(u.mas/u.yr).to_string(\"latex_inline\")}]')\n", + "ax.legend()\n", + "plt.draw()" + ] + } + ], + "metadata": { + "language_info": { + "name": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}