astropy / astropy

Astronomy and astrophysics core library
https://www.astropy.org
BSD 3-Clause "New" or "Revised" License
4.38k stars 1.75k forks source link

Velocities in cylindrical coordinates #6425

Closed jobovy closed 7 years ago

jobovy commented 7 years ago

I'm trying to compute velocities in the Galactocentric cylindrical coordinate frame with astropy v2, but as far as I can tell this isn't documented anywhere. The following seems to work

ra= 120.*u.deg
dec= 30.*u.deg
distance= 1.2*u.kpc
pmra= 5.*u.mas/u.yr
pmdec= -3.*u.mas/u.yr
vlos= 55.*u.km/u.s
c= apycoords.ICRS(ra=ra,dec=dec,
                  distance=distance,
                  pm_ra_cosdec=pmra,
                  pm_dec=pmdec,
                  radial_velocity=vlos)
# Define GC frame
v_sun = apycoords.CartesianDifferential([11.1,245.,7.25]*u.km/u.s)
gc_frame= apycoords.Galactocentric(galcen_distance=8.1*u.kpc,
                                   z_sun=25.*u.pc,
                                   galcen_v_sun=v_sun)
cg= c.transform_to(gc_frame)
cg.representation= 'cylindrical'
print("Cylindrical: v_r/(km/s)=",cg.d_rho.to(u.km/u.s),
      ", v_phi/(km/s)=",cg.d_phi.to(u.mas/u.yr).value*4.74047*cg.rho.to(u.kpc).value,
      " v_z/(km/s)=",cg.d_z.to(u.km/u.s))
Cylindrical: v_r/(km/s)= 17.213193237684614 km / s , v_phi/(km/s)= -213.59390592255582  v_z/(km/s)= 52.07396464035026 km / s

This gives the right answer, but maybe this isn't how these d_ attributes are supposed to be used? Also, having to manually convert d_phi from mas/yr to km/s must be wrong. Is there a way to get the velocity in the cylindrical representation out more directly that I am missing or if not, is it correct to use the d_ attributes in this way?

mhvk commented 7 years ago

@jobovy - thanks! This is all new and shiny, so things may not be as clear yet as they should be. I'm happy to see though that the values come out right -- I don't think we have many tests for cylindrical representations...

Here, I think the code works as planned: you ask it to represent the coordinates in cylindrical, i.e., a radius, angle, and height, and then the time derivatives will use the same units, so for the angle you get angle/time - it effectively is a proper motion, just as you would get if you went back to a spherical coordinate system (I wonder if it would be clearer if we called it pm_phi, though on the other hand I'm a bit wary of too many name overrides...)

As for the manual conversion, after multiplication with distance, the unit system needs to know how to get rid of the angle. Since the small angle approximation is a good one, we can treat angles in radians as dimensionless:

(cg.d_phi * cg.rho).to(u.km/u.s, equivalencies=u.dimensionless_angles())
# <Quantity -213.5939268065537 km / s>

(This is equivalent to cg.d_phi/u.rad*cg.rho)

jobovy commented 7 years ago

@mhvk -- thanks! I think I understand it a little better now. I guess I wasn't expecting to get an answer in mas/yr when converting a velocity to the Galactocentric frame, I think rad/yr would be more intuitive (nobody is measuring these angular velocities in mas/yr), because then you don't have to bother with unit equivalencies to get an answer that makes sense.

Feel free to close this issue, because I don't think there's anything actually incorrect here. I must say that after reading through the documentation (although without digging in) it wasn't clear to me how to get these velocities out of this representation. Are you sure you want to use the notation d_X for something that actually is a differential of X divided by time? v_X would seem to make more sense to me, or d_X_dt perhaps.

mhvk commented 7 years ago

@jobovy - you actually get mas because your original unit was mas... Internally, units get simplified only when specifically asked (indeed, even the internal representation is not actually changed when you set `cg.representation='cylindrical'). To see this,

In [13]: cg.data
Out[13]: 
<CartesianRepresentation (x, y, z) in kpc
    (-9.14615185, -0.20890406,  0.57440393)
 (has differentials w.r.t.: 's')>

In [14]: cg.data.differentials['s']
Out[14]: 
<CartesianDifferential (d_x, d_y, d_z) in kpc mas / (rad yr)
    (-4.65904331,  44.96287378,  10.98497819)>

Basically, what happened is that when you created a cartesian frame, the proper motion was converted to a velocity by multiplying by the distance and dividing by rad. Then, we you convert back to an angle, the number gets divided by a distance, and multiplied by rad, leaving mas/yr.

In general, I found one of the things that took longest is to stop worrying about what unit a number is in, and instead just convert when needed, counting on the unit system to keep things sane. But for looking at values, this is not always helpful, which is why the coordinates have "preferred units" and hence when you look at just cg, you get a velocity in km/s rather than kpc mas/(rad yr). The cylindrical system for galactocentric is sufficiently rare that, so far, nobody has bothered to set similar "preferred solutions units" (not sure everybody would agree on rad/yr, though...).

mhvk commented 7 years ago

I think the question of naming is a good one, and probably best addressed as a separate issue, so I opened #6434.