add python astropy notebook

edits

clear outputs

add source

fix name typo

Author: seeM

requirements.txt

@@ -41,4 +41,5 @@ nbclient
 torch
 SQLAlchemy
 psycopg2-binary
-ipykernel
+ipykernel
+astropy

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
+}