Skip to content

A pretty Python Jupyter notebook for screenshots #60

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 1 commit into from
Jul 8, 2025
Merged

A pretty Python Jupyter notebook for screenshots #60

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion requirements.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,4 +41,5 @@ nbclient
torch
SQLAlchemy
psycopg2-binary
ipykernel
ipykernel
astropy
335 changes: 335 additions & 0 deletions workspaces/galactocentric-frame/galactocentric-frame.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
}