question about obs argument in Orbit methods

[BrianTCook]

I am trying to create lists of dists/ras/decs for an orbit in the following way:

o_forward = Orbit(
vxvv=[R, vR, vT, z, vz, phi],
ro=8.23,
vo=249.44
)

o_forward.integrate(
ts,
MWPotential2014,
method='leapfrog'
)

o_backward = Orbit(
vxvv=[R, -vR, -vT, z, -vz, phi],
ro=8.23,
vo=249.44
)

o_backward.integrate(
ts,
MWPotential2014,
method='leapfrog'
)

x_unit = u.kpc
v_unit = u.km / u.second

#TODO: FIGURE OUT WHAT NEEDS TO GET PASSED IN
obs_sun_3d = [-8.23 * x_unit, 0. * x_unit, 0.0208 * x_unit]
obs_sun_6d = [-8.23 * x_unit, 0. * x_unit, 0.0208 * x_unit, -11.1 * v_unit, 12.24 * v_unit, 7.25 * v_unit]

dist_values_orbit = np.concatenate(
(np.asarray(o_backward.dist(ts, obs=obs_sun_3d, ro=8.23, vo=249.44)),
np.asarray(o_forward.dist(ts, obs=obs_sun_3d, ro=8.23, vo=249.44)))
)

ra_values_orbit = np.concatenate(
(np.asarray(o_backward.ra(ts, obs=obs_sun_3d, ro=8.23, vo=249.44)),
np.asarray(o_forward.ra(ts, obs=obs_sun_3d, ro=8.23, vo=249.44)))
)

dec_values_orbit = np.concatenate(
(np.asarray(o_backward.dec(ts, obs=obs_sun_3d, ro=8.23, vo=249.44)),
np.asarray(o_forward.dec(ts, obs=obs_sun_3d, ro=8.23, vo=249.44)))
)

where the [R, vR, vT, z, vz, phi] values are from a particular object of interest and the two Orbit objects are meant to find its past/future trajectories. I see in line 3565 of the Orbits.py script that the dec method accepts an optional (keyword?) argument that defines the position of the observer. I get my Solar phase values from schonrich2010 and leung2022. I followed the stack trace of the ra()/dec()/dist() methods to https://github.com/jobovy/galpy/blob/v1.8.1/galpy/util/coords.py and I'm having a hard time with the following things: determining when the obs list elements needs Astropy units attached to them, and when obs is [x,y,z,vx,vy,vz] vs. [ro,vo,zo,U,V,W]. Any guidance would be appreciated, thanks!

[jobovy]

Apologies for the slow reply, busy day... There are two answers to this question:

• I think what you are trying to do is to determine the (ra,dec.dist) along the forward and backward path of an orbit started at a given point? That you could do simply as

  import numpy
  from galpy.potential import MWPotential2014
  from galpy.orbit import Orbit
  import astropy.units as u
  R,vR,vT,z,vz,phi= 8.*u.kpc,20.*u.km/u.s,220.*u.km/u.s,0.*u.kpc,30.*u.km/u.s,0.*u.rad
  of= Orbit([R,vR,vT,z,vz,phi],ro=8.23,vo=249.44,zo=0.0208,solarmotion=[-11.1,24.,7.25])
  ts= numpy.linspace(0.,1.,1001)*u.Gyr
  of.integrate(ts,MWPotential2014)
  ob= of() # copy for backwards integration
  ob.integrate(-ts,MWPotential2014)
  # Concatenate the integrated orbits to easily get values along the forward and backward orbit
  o= Orbit([*of(ts),*ob(-ts)],ro=8.23,vo=249.44,zo=0.0208,solarmotion=[-11.1,24.,7.25]))
  o.ra() # has *all* the RA coordinates ordered forward first, then backwards
  o.dec() # all the Dec coordinates
  o.dist() # all the distances
  

  When just calling o.ra() etc. this way, the observer coordinates are taken from the Orbit instance and so are specified by the ro, vo, zo, and solarmotion parameters. I've listed them above, even though the zo and solarmotion ones are just the defaults. The solar motion is specified as the Sun's [-U,V,W] in the standard Galactic coordinates. For convenience, I've concatenated the forward and backward integrated orbits into a new Orbit instance that then has all forward and backward points as initial conditions, so you can get all ra as o.ra().

• In your code, you are trying to specify the observer coordinates directly when calling o.ra() etc. Exactly what obs should be is perhaps easiest to explain by giving the defaults taken from the Orbit instance: obs= [ro,0.,zo,solarmotion[0],solarmotion[1]+vo,solarmotion[2]]. It's essentially the Sun's position and velocity in the Galactocentric coordinate frame, although it's not actually that exactly because it doesn't account for the velocities for the fact that the Galactocentric and Galactic coordinate frames aren't exactly aligned in their z=0 planes. But that's a very small effect. When you only ask for positions, you can just specify obs=[X_sun,Y_sun,Z_sun], there's no need to specify the Sun's velocity. obs is assumed to have units of [kpc,kpc,kpc,km/s,km/s,km/s], but you can also explicitly use units, e.g., o.dist(obs=[8.*u.kpc,0.*u.kpc,20.8*u.pc]).

Hope this helps! I'm happy to follow-up with more explanation if necessary. It's generally easiest for me if you can send a full example that I can directly run (so has all imports and values specified).

[BrianTCook]

Thanks for this reply, I really appreciate it! I will implement your suggested changes and follow up with the result in the next few days. I can send a link to the relevant repo privately along with additional comments/context.