Discrepancy with Galpy when using agama.GalpyPotential

[brugalada]

Discrepancy with Galpy when using agama.GalpyPotential

[eugvas]

I think it's a confusion with units.
If you plot the circular-velocity curve of the original potential, plt.plot(r, (-r*galapot.Rforce(r, r*0))**0.5), it has a peak around 6, same as the red curve in your plot. Since the NFW potential has a peak circular velocity at ~2.2 times its scale length, this corresponds to 55 kpc, or ~6.2 in the "native" Galpy length units of 8 kpc. But if you call galapot.Rforce(r * u.kpc, r*0), it apparently perform the conversion from astropy dimensional units to its internal units.
Now it seems to me that even when you call galapot.plotRotcurve with a plain numpy array, it still assumes that the radii are provided in kpc, not in units of 8kpc, that's why the blue curve reaches the peak well beyond the right margin of the plot.
And if you create a pure Agama NFW potential with a scale length of 25 (presumably, kpc, though it is never specificed), its circular velocity will peak at 2.2 times that value, as expected - that's the green dashed line in your plot (mislabelled as "from Galpy", whereas it comes from a pure Agama potential).
I don't understand how Galpy deals with units when mixing astropy.units with its native ones, so it seems that the most straightforward (if awkward) approach is to stick to the length unit of 8kpc and velocity unit of 220 km/s (and not to forget to call agama.setUnits(length=8, velocity=220, mass=1) at the beginning of the script to instruct Agama to use the same units). Then if you interpret the x axis in terms of units of 8kpc, not 1 kpc, everything should match and peak at x=6.2

[brugalada]

First of all, yep, you are right: I put the labels in the wrong order.

Having said that, I think it is not worth losing any sleep over this since it is indeed a problem with units, and Galpy is too opaque for me to try and figure it out.

Thank you for the clarification! And sorry for the silly question jeje