I get different results using galpy' internal units and physical units. Which result should I trust?

[jiweiliao6]

l calculate the closest distance of a star's orbit to sgr dsph's orbit. ButI l get different results using galpy' internal units and physical units. Which result should I trust?

Here is code:

import numpy as np
from numpy import sin,cos,deg2rad
from galpy.orbit import Orbit
import galpy.util.coords
from astropy import units
from galpy.potential import MWPotential2014

'''Cartesian Galactocentric coordinate system'''

R_sun, z_sun = 8.122, 20.8e-3   #sun position, units.kpc
def xyz(l,b,distance,rad=False):
    if rad==False:
        l,b = deg2rad(l), deg2rad(b)
    x = R_sun-distance*cos(b)*cos(l)
    y = -distance*cos(b)*sin(l)
    z = distance*sin(b)+z_sun
    return x,y,z

v_sun = np.array((11.1,245,7.25))    #sun velocity, units.km/s
def UVW(ra,dec,rv,dist,pmra,pmdec,rad=False):
    if rad==False:
        ra,dec = deg2rad(ra), deg2rad(dec)
    T = np.array((-0.06699, -0.87276, -0.48354, +0.49273, -0.45035, +0.74458, -0.86760, -0.18837, +0.46020)).reshape((3,3))
    A =np.array((cos(ra)*cos(dec), -sin(ra), -cos(ra)*sin(dec), +sin(ra)*cos(dec), +cos(ra), -sin(ra)*sin(dec), +sin(dec), 0, +cos(dec))).reshape((3,3))
    B = np.dot(T,A)
    V = np.array((rv,4.74047*dist*pmra,4.74047*dist*pmdec))
    uvw = np.dot(B,V)
    U,V,W = uvw[0]+v_sun[0],uvw[1]+v_sun[1],uvw[2]+v_sun[2]
    return U,V,W

'''Orbit initialization'''
def RvRvTzvzphi(XYZ,uvw):
    R = np.sqrt((XYZ[0])**2 + (XYZ[1])**2)
    vR = -uvw[0]*XYZ[0]/R - uvw[1]*XYZ[1]/R
    vT = uvw[0]*XYZ[1]/R - uvw[1]*XYZ[0]/R
    z = XYZ[2]
    vz = uvw[2]
    phi = np.arctan(XYZ[1]/XYZ[0])
    return R,vR,vT,z,vz,phi

'''Sgr dsph's orbit'''
s_ra,s_dec,s_dist,s_pmra,s_pmdec,s_rv = 283.83129,-30.5453,26.669,-2.692,-1.3589,139.99
s_lb = galpy.util.coords.radec_to_lb(s_ra,s_dec,degree=True)
s_XYZ = xyz(s_lb[0],s_lb[1],s_dist)
s_uvw = UVW(s_ra,s_dec,s_rv ,s_dist,s_pmra,s_pmdec)


s_R,s_vR,s_vT,s_z,s_vz,s_phi =  RvRvTzvzphi(s_XYZ,s_uvw)

s_vxvv = [s_R/8. ,s_vR/220. ,s_vT/220., s_z/8., s_vz/220., s_phi*180./np.pi]  #conversion galpy internal units
o_s = Orbit(s_vxvv,ro=8.,vo=220.)
ts= np.linspace(0.,-1.,10000)*units.Gyr
o_s.integrate(ts,MWPotential2014)
x_sgr,y_sgr,z_sgr = o_s.x(ts),o_s.y(ts),o_s.z(ts)

#use physical untis
o_s1 = Orbit([s_R*units.kpc,s_vR*units.km/units.s,s_vT*units.km/units.s,s_z*units.kpc,s_vz*units.km/units.s,s_phi*180./np.pi*units.deg])
o_s1.integrate(ts,MWPotential2014)
x_sgr1,y_sgr1,z_sgr1 = o_s1.x(ts),o_s1.y(ts),o_s1.z(ts)


'''a star's orbit'''
dist,pmra,pmdec = 1/0.395837,-46.914,-1.465
va,vb = 4.74047*dist*pmra,4.74047*dist*pmdec
l,b,ra,dec,rv = 13.466301829915965, 61.23809520105212, 220.85733333, 14.89341667,192
XYZ,uvw = xyz(l,b,dist),UVW(ra,dec,rv ,dist,pmra,pmdec)


R,vR,vT,z,vz,phi =  RvRvTzvzphi(XYZ,uvw)

vxvv = [R/8. ,vR/220. ,vT/220., z/8., vz/220., phi*180./np.pi]  #conversion galpy internal units
o = Orbit(vxvv,ro=8.,vo=220.)
o.integrate(ts,MWPotential2014)
o_x,o_y,o_z = o.x(ts),o.y(ts),o.z(ts)

#use physical untis
o1 = Orbit([R*units.kpc,vR*units.km/units.s,vT*units.km/units.s,z*units.kpc,vz*units.km/units.s,phi*180./np.pi*units.deg])
o1.integrate(ts,MWPotential2014)
o_x1,o_y1,o_z1 = o1.x(ts),o1.y(ts),o1.z(ts)
        

'''Calculate the closest distance of a star's orbit to sgr dsph's orbit '''

internal,physical = [],[]
for k in range(len(o_x)):
    rsgr_rstar = np.sqrt((o_x[k]-x_sgr[k])**2+(o_y[k]-y_sgr[k])**2+(o_z[k]-z_sgr[k])**2)
    internal.append(rsgr_rstar) 
    
    rsgr_rstar1 = np.sqrt((o_x1[k]-x_sgr1[k])**2+(o_y1[k]-y_sgr1[k])**2+(o_z1[k]-z_sgr1[k])**2)
    physical.append(rsgr_rstar1) 
    
print('galpy internal units,',np.min(internal))   #5.5kpc
print('use physical units,',np.min(physical))     #13.4kpc

[jiweiliao6]

@jobovy
Hi bovy.
I'm confused and don't know where is the problem, could you tell me?

[jobovy]

I am away without a laptop, so can only look at this closely next week. But I notice that you seem to convert phi to degree when you convert to internal units? Phi should be in rad in internal units.

[jiweiliao6]

Phi is the internal unit rad or deg, the result is the same.

Thank you for your reply, looking forward to your results

[jobovy]

I just looked at this more closely and it's exactly as I said before: the problem is that when you setup the orbit in internal units, you give phi in degree rather than in radian. That is, in these lines:

s_vxvv = [s_R/8. ,s_vR/220. ,s_vT/220., s_z/8., s_vz/220., s_phi*180./np.pi]  #conversion galpy internal units
...
vxvv = [R/8. ,vR/220. ,vT/220., z/8., vz/220., phi*180./np.pi]  #conversion galpy internal units

They should just be

s_vxvv = [s_R/8. ,s_vR/220. ,s_vT/220., s_z/8., s_vz/220., s_phi]  #conversion galpy internal units
...
vxvv = [R/8. ,vR/220. ,vT/220., z/8., vz/220., phi]  #conversion galpy internal units

Then both methods agree and the result is 13.4 kpc.

[jiweiliao6]

Sorry for my careless, the previous validation was for radians and degree in the physical units rather than the internal units. The internal unit is radians to get the same result 13.4kpc. Thank you so much for taking the time to check. I will express my gratitude to you in the acknowledgments