Skip to content

Instantly share code, notes, and snippets.

@jonathansick
Forked from PBarmby/gal_radii_pb.py
Last active April 5, 2024 18:23
Show Gist options
  • Save jonathansick/9399842 to your computer and use it in GitHub Desktop.
Save jonathansick/9399842 to your computer and use it in GitHub Desktop.
Compute deprojected galactocentric distance using astropy.coordinates. (By @PBarmby)
from astropy.coordinates import ICRS, Distance, Angle
from astropy import units as u
import numpy as np
def correct_rgc(coord, glx_ctr=ICRS('00h42m44.33s +41d16m07.5s'),
glx_PA=Angle('37d42m54s'),
glx_incl=Angle('77.5d'),
glx_dist=Distance(783, unit=u.kpc)):
"""Computes deprojected galactocentric distance.
Inspired by: http://idl-moustakas.googlecode.com/svn-history/
r560/trunk/impro/hiiregions/im_hiiregion_deproject.pro
Parameters
----------
coord : :class:`astropy.coordinates.ICRS`
Coordinate of points to compute galactocentric distance for.
Can be either a single coordinate, or array of coordinates.
glx_ctr : :class:`astropy.coordinates.ICRS`
Galaxy center.
glx_PA : :class:`astropy.coordinates.Angle`
Position angle of galaxy disk.
glx_incl : :class:`astropy.coordinates.Angle`
Inclination angle of the galaxy disk.
glx_dist : :class:`astropy.coordinates.Distance`
Distance to galaxy.
Returns
-------
obj_dist : class:`astropy.coordinates.Distance`
Galactocentric distance(s) for coordinate point(s).
"""
# distance from coord to glx centre
sky_radius = glx_ctr.separation(coord)
avg_dec = 0.5 * (glx_ctr.dec + coord.dec).radian
x = (glx_ctr.ra - coord.ra) * np.cos(avg_dec)
y = glx_ctr.dec - coord.dec
# azimuthal angle from coord to glx -- not completely happy with this
phi = glx_PA - Angle('90d') \
+ Angle(np.arctan(y.arcsec / x.arcsec), unit=u.rad)
# convert to coordinates in rotated frame, where y-axis is galaxy major
# ax; have to convert to arcmin b/c can't do sqrt(x^2+y^2) when x and y
# are angles
xp = (sky_radius * np.cos(phi.radian)).arcmin
yp = (sky_radius * np.sin(phi.radian)).arcmin
# de-project
ypp = yp / np.cos(glx_incl.radian)
obj_radius = np.sqrt(xp ** 2 + ypp ** 2) # in arcmin
obj_dist = Distance(Angle(obj_radius, unit=u.arcmin).radian * glx_dist,
unit=glx_dist.unit)
# Computing PA in disk (unused)
obj_phi = Angle(np.arctan(ypp / xp), unit=u.rad)
# TODO Zero out very small angles, i.e.
# if np.abs(Angle(xp, unit=u.arcmin)) < Angle(1e-5, unit=u.rad):
# obj_phi = Angle(0.0)
return obj_dist
@jonathansick
Copy link
Author

Updated @PBarmby's gist to deal with input coordinates that are arrays (use numpy for all math). A bit of PEP8 too :)

@Gabriel-p
Copy link

Shouldn't this line

obj_dist = Distance(Angle(obj_radius, unit=u.arcmin).radian * glx_dist, unit=glx_dist.unit)

be

obj_dist = Distance(np.tan(Angle(obj_radius, unit=u.arcmin).radian) * glx_dist, unit=glx_dist.unit)

?

That is, shouldn't the obj_radius angle be inside a tan?

See e.g. Eq 5 in The metallicity gradient as a tracer of history and structure: the Magellanic Clouds and M33 galaxies, Cioni 2009

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment