Can't validate resultes

[edan0314]

I'm new to Galpy and trying to validate the results from the following article:
https://arxiv.org/pdf/2007.02350.pdf
Table 4 shows the resultsobtained after integrating 13.5 Gyrs in the past.
I'm trying to reconstitute the results for the dwarf galaxy of Sagitarius I.
The author make orbit integration for 3 poptentials, teo of then been modified versions of Allen & Santillan Milk-Way potential and modified version of Navarro-Frenk-White and a third beeing the Allen and Santillan versio of Irrgang.
Since Galpy provides a model of Allen & Santillan modified by Irrgang as Irrgang13I, I've been trying to integrate over the same potential however the results are too different, and I do not know what is going on

Here's my code

#importing functions
from galpy.potential import (plotRotcurve, MWPotential as mwp, MWPotential2014 as mwp14, MiyamotoNagaiPotential, PlummerPotential,NFWPotential)
import matplotlib.pyplot as plt
from galpy.potential.mwpotentials import Irrgang13I as asi1, Irrgang13II as asi2, McMillan17 as mm17
import numpy
from galpy.util.conversion import get_physical
from astropy import units
from galpy.util.coords import XYZ_to_galcencyl, vxvyvz_to_galcencyl, cyl_to_rect, rect_to_cyl,rect_to_cyl_vec
from galpy.orbit import Orbit

#Position of SgrI w.r.t. Sun
X=25.2-8.33
Y=2.2
Z=-6.2+16/1000
R=numpy.sqrt((X**2)+(Y**2))
rgc=numpy.sqrt((X**2)+(Y**2)+(Z**2))
print(rgc)
print("Galatocentric rectangular position vector:",[X,Y,Z])
Rv=X*(X/R)+Y*(Y/R)
phiv=numpy.arctan(Y/X)
print("Galatocentric cylindrical position vector",[Rv,phiv,Z])
Rprime=rect_to_cyl(X,Y,Z)
print("posotion vector according to galpy's function:",Rprime)
print("#################################################################")

# velocity field of SgrI w.r.t. Sun
u0m=217.6
v0m=-270.7
w0m=201.7
#galatocentric rectangular velocity field of SgrI 
U=11.1+u0m
V=12.2+242+v0m
W=7.3+w0m
Pi=-U*(X/R)+V*(Y/R)
Theta=U*(Y/R)+V*(X/R)
print("cylindrical velocity field:",[Pi,Theta,W])
print("#################################################################")

timesSgrI= numpy.linspace(0,-13.5,100001)*units.Gyr
o_asi1_SgrI= Orbit([Rv*units.kpc,Pi*units.km/units.s,Theta*units.km/units.s,Z*units.kpc,W*units.km/units.s,phiv*units.deg],ro=8.4,vo=242)
o_asi1_SgrI.integrate(timesSgrI,asi1,method="rk4_c")
print("position vector:",[o_asi1_SgrI.R(-13.5*units.Gyr),o_asi1_SgrI.phi(-13.5*units.Gyr),o_asi1_SgrI.z(-13.5*units.Gyr)],"kpc")
print("distance from GC:",o_asi1_SgrI.r(-13.5*units.Gyr),"kpc")
print("velocity vector:",[o_asi1_SgrI.vR(-13.5*units.Gyr),o_asi1_SgrI.vT(-13.5*units.Gyr),o_asi1_SgrI.vz(-13.5*units.Gyr)],"km/s")
print("total velocity:",numpy.sqrt((o_asi1_SgrI.vR(-13.5*units.Gyr)**2)+(o_asi1_SgrI.vT(-13.5*units.Gyr)**2)+(o_asi1_SgrI.vz(-13.5*units.Gyr)**2)),"km/s")
o_asi1_SgrI.plot()

I know Galpy posses functions in which are possible to transform coordinates the handness of Galpy's coordiantes are different from the one used in the article, and since I'm trying to reproduce the results, I'm doing (or trying to do) the most similar way as possible according to the relations given in the article.

[jobovy]

Happy to help here, but I'm not sure what the problem is. I ran your code and the results seem plausibly similar to the ones in the paper you link to. Which differences are you concerned about in particular? I think their Vtot is the one today, not 13.5 Gyr ago (Vtot today ~ 310 km/s as they list in their table). I get similar peri- and apo-center radii and similar energy. Lz is off by a factor of two, but I think that is probably because galpy uses a different default solar motion than they are using. Lz = R x vT and vT here ~ 13 km/s, so if the solar motion is different by ~13 km/s one could end up with twice as large Lz like what they have.

Note that Sgr is also part of the objects that you can call with Orbit.from_name:

o= Orbit.from_name('Sgr')

This will use slightly different current position/velocity than you are using, but you can use this to check, because present day coordinates and orbital properties should be quite similar.

[edan0314]

Table 4 says:
"Dwarf orbital properties. For each dwarf we quote values derived from orbits integrated for 13.5 Gyr
backward in three potentials of Models NFWBB, ASBB and ASI in the first, second, and third row, respectively."
I understood that the values in this table are the reults after the integration.
In this case, the results I obtain after the intagration:

timesSgrI= numpy.linspace(0,-13.5,100001)units.Gyr
o_asi1_SgrI= Orbit([Rvunits.kpc,Piunits.km/units.s,Thetaunits.km/units.s,Zunits.kpc,Wunits.km/units.s,phivunits.deg],ro=8.4,vo=242)
o_asi1_SgrI.integrate(timesSgrI,asi1,method="rk4_c")
print("position vector:",[o_asi1_SgrI.R(-13.5units.Gyr),o_asi1_SgrI.phi(-13.5units.Gyr),o_asi1_SgrI.z(-13.5units.Gyr)],"kpc")
print("distance from GC:",o_asi1_SgrI.r(-13.5units.Gyr),"kpc")
print("velocity vector:",[o_asi1_SgrI.vR(-13.5units.Gyr),o_asi1_SgrI.vT(-13.5units.Gyr),o_asi1_SgrI.vz(-13.5units.Gyr)],"km/s")
print("total velocity:",numpy.sqrt((o_asi1_SgrI.vR(-13.5units.Gyr)**2)+(o_asi1_SgrI.vT(-13.5units.Gyr)**2)+(o_asi1_SgrI.vz(-13.5*units.Gyr)**2)),"km/s")
o_asi1_SgrI.plot()

in this part, I do not get the same results as table 4.
But what I understood is that you're saying tha values from table 4 are vslues in the present?

[jobovy]

Hi,

Please remember that I did not write the paper that you linked to and I haven't even read it. The authors don't even use galpy I think and even if they did, I'm not responsible for knowing exactly what they did. I can't tell you what they mean when they write

For each dwarf we quote values derived from orbits integrated for 13.5 Gyr

but I interpret this as meaning that properties such as the eccentricity and the apo/pericenter radii are obtained from the orbit integration, while positions and velocities are just those today. And Vtot is exactly what it is today; it would be quite coincidental that it is also the value 13.5 Gyr ago. But I can't be sure about what they did. All I'm saying is that from my own orbit integration, all of the values they listed plausibly agree with the orbit integration that you did, modulo that I don't know what they used for the solar motion or the Sun's height above the mid-plane.

If your further questions relate to the interpretation of the table, then I suggest you contact the authors of the paper.

[edan0314]

I'm aware you're not responsible by the article, but I thougth you're knowledge would help, and helped a lot for me, because you confirmed that the results I just found and some previous are correct, if one interprets the values on that table of the velocity field as the values in the present. For a moment a was wondering what I did wrong. I've also tryed to do the same orbital integration for SgrII and PisII and the same happened, the values of the position and velocities are the same as the present.

[jobovy]

Okay, great that this was helpful, as I said, I'm happy to help! If my first answer essentially answered the question, please mark it as the answer to other people looking at the discussions can see that this was answered. Thanks!