Address 5 review findings

1. (bug) progpot + streamTrack: save the pre-progpot potential as
   _orig_pot and use it for the track-progenitor integration, avoiding
   the MovingObjectPotential whose internal orbit is only valid on
   [-tdisrupt, 0].

2. (bug) vlos double-scaling: _helio_xv already returns km/s (it pre-
   multiplies by vo), so the 'vlos' branch of _scale now only attaches
   units without re-scaling.

3. (docs) Fix docstrings: track_time_range default (now auto-timerange),
   track_n_dense actual grid size, smoothing description (now GCV), add
   ntp kwarg docs.

4. (docs) Fix class docstring tp range description (not [-tdisrupt, 0]),
   add arm_sign and ntp to Parameters section.

5. (cleanup) plot spread band: remove dead s1 variable, gate physical
   scaling on self._physical only (consistent with _scale), drop the
   d1-in-cart_idx requirement since the band is only in the d2 direction.

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

Author: jobovy

galpy/df/streamTrack.py

@@ -120,9 +120,11 @@ class StreamTrack:
 
     A StreamTrack holds a smooth mean curve (and optional covariance) in
     galactic phase space, parameterized by a progenitor time coordinate
-    ``tp`` in ``[-tdisrupt, 0]``. It is constructed from a cloud of stream
-    particles (e.g. drawn via ``basestreamspraydf.sample``) plus the
-    progenitor's orbit.
+    ``tp``. ``tp=0`` is the progenitor today; for a leading arm ``tp > 0``
+    (future positions the progenitor has yet to reach) and for a trailing
+    arm ``tp < 0``. The ``tp`` range is determined by the data: it spans
+    the percentile-trimmed range of the closest-point assignments, bounded
+    by ``track_time_range`` (typically much smaller than ``tdisrupt``).
 
     Notes
     -----
@@ -173,15 +175,21 @@ def __init__(
         track_t_grid : array, shape (M,)
             The dense time grid on which ``track_prog_cart`` is evaluated.
             Used for the closest-point projection.
+        arm_sign : int, optional
+            ``+1`` for leading arm (tp >= 0), ``-1`` for trailing (tp <= 0).
+            Controls the sign constraint on the closest-point search window.
+        ntp : int, optional
+            Number of binning nodes. Default ``sqrt(N)`` clipped to
+            ``[21, 201]``.
         ninterp : int, optional
             Resolution of the fine tp grid on which the public track is
             stored.
         smoothing : None, float, or dict, optional
             Smoothing parameter for the mean spline per Cartesian coordinate.
-            ``None`` estimates s from the pairwise noise variance of the
-            sorted particle sample (Reinsch-like target s = N * sigma^2).
-            A float sets s for all six coordinates; a dict keyed by
-            ``'x','y','z','vx','vy','vz'`` sets per-coordinate values.
+            ``None`` (default) uses GCV auto-tuning via
+            ``scipy.interpolate.make_smoothing_spline``. A float sets an
+            explicit ``s`` for ``UnivariateSpline``; a dict keyed by
+            ``'x','y','z','vx','vy','vz'`` sets per-coordinate ``s`` values.
         niter : int, optional
             Iterations beyond the initial fit. Each iteration reassigns each
             particle to the closest point on the current track and refits.
@@ -410,8 +418,12 @@ def _scale(self, val, kind):
             return val
         if kind == "length":
             return val * self._ro * (units.kpc if _APY_UNITS else 1)
-        if kind == "velocity" or kind == "vlos":
+        if kind == "velocity":
             return val * self._vo * (units.km / units.s if _APY_UNITS else 1)
+        if kind == "vlos":
+            # vlos from _vrpmllpmbb is already in km/s (helio_xv pre-
+            # multiplies by vo); only attach units, don't re-scale.
+            return val * (units.km / units.s if _APY_UNITS else 1)
         if kind == "angle":
             return val * (units.rad if _APY_UNITS else 1)
         if kind == "degree":
@@ -674,21 +686,18 @@ def plot(self, d1="x", d2="y", spread=0, n=None, **kwargs):
         line = pyplot.plot(v1, v2, **kwargs)
         if spread > 0 and self._cov_xyz is not None:
             cart_idx = {"x": 0, "y": 1, "z": 2, "vx": 3, "vy": 4, "vz": 5}
-            if d1 in cart_idx and d2 in cart_idx:
-                i1 = cart_idx[d1]
+            if d2 in cart_idx:
                 i2 = cart_idx[d2]
                 cov = self.cov(tp)  # (n_eval, 6, 6)
-                s1 = numpy.sqrt(numpy.maximum(cov[:, i1, i1], 0.0))
                 s2 = numpy.sqrt(numpy.maximum(cov[:, i2, i2], 0.0))
-                # Apply physical scaling if active
-                if self._physical and self._roSet and d1 in ("x", "y", "z"):
-                    s1 = s1 * self._ro
-                if self._physical and self._roSet and d2 in ("x", "y", "z"):
-                    s2 = s2 * self._ro
-                if self._physical and self._voSet and d1 in ("vx", "vy", "vz"):
-                    s1 = s1 * self._vo
-                if self._physical and self._voSet and d2 in ("vx", "vy", "vz"):
-                    s2 = s2 * self._vo
+                # Scale covariance sigma the same way _scale handles the
+                # coordinate values (check _physical only, consistent with
+                # the coordinate accessors).
+                if self._physical:
+                    if d2 in ("x", "y", "z"):
+                        s2 = s2 * self._ro
+                    elif d2 in ("vx", "vy", "vz"):
+                        s2 = s2 * self._vo
                 color = line[0].get_color() if line else None
                 pyplot.fill_between(
                     v1,

galpy/df/streamspraydf.py

@@ -133,6 +133,7 @@ def __init__(
         else:
             self._center = None
         if progpot is not None:
+            self._orig_pot = self._pot  # save pre-progpot for streamTrack
             progtrajpot = MovingObjectPotential(
                 orbit=self._progenitor,
                 pot=progpot,
@@ -248,18 +249,25 @@ def streamTrack(
             the value set at initialization.
         track_time_range : float or Quantity, optional
             Half-range (symmetric about tp=0) of the finely-integrated
-            progenitor orbit used for closest-point matching and as the
-            public tp grid. Default is ``0.1*tdisrupt``. Must be larger than
-            the stream's spatial extent in orbital-time units.
+            progenitor orbit used for closest-point matching. Default is
+            data-driven: ``8 * d_max / |v_prog|`` clamped to ``[1,
+            tdisrupt]``, where ``d_max`` is the farthest particle's
+            distance from the progenitor.
         track_n_dense : int, optional
             Number of time points on the finely-integrated progenitor
-            orbit. Default 10001.
+            orbit. The actual grid has ``2 * (track_n_dense+1)//2 - 1``
+            points (forward + backward with shared t=0). Default 10001.
+        ntp : int, optional
+            Number of binning nodes. Default ``sqrt(N)`` clipped to
+            ``[21, 201]``.
         ninterp : int, optional
             Resolution of the public fine-grid track arrays. Default 1001.
         smoothing : None, float, or dict, optional
-            Smoothing spline ``s`` parameter for the mean-track fit.
-            ``None`` estimates s from the pairwise noise variance of the
-            sorted particle sample per coordinate.
+            Smoothing parameter. ``None`` (default) uses GCV auto-tuning
+            via ``scipy.interpolate.make_smoothing_spline``. A float sets
+            an explicit ``s`` value for ``UnivariateSpline``; a dict
+            keyed by ``'x','y','z','vx','vy','vz'`` sets per-coordinate
+            ``s`` values.
         niter : int, optional
             Iterations beyond the initial fit. Each iteration reassigns
             particles to the closest point on the current track.
@@ -323,15 +331,22 @@ def streamTrack(
         # Build a finely-sampled progenitor phase-space array spanning
         # [-T, +T] around the present day. Integrate forward and backward
         # separately from the progenitor's present-day state, then combine.
-        half_dense = max(2, int(track_n_dense) // 2)
+        # Use the BASE potential (without progpot's MovingObjectPotential)
+        # because the MovingObjectPotential's internal progenitor was only
+        # integrated on [-tdisrupt, 0] and would give wrong or erroring
+        # results for positive times.
+        _track_pot = self._pot
+        if hasattr(self, "_orig_pot"):
+            _track_pot = self._orig_pot
+        half_dense = (int(track_n_dense) + 1) // 2
         t_back = numpy.linspace(0.0, -track_time_range, half_dense)
         t_fwd = numpy.linspace(0.0, track_time_range, half_dense)
         prog_back = self._orig_progenitor()
         prog_back.turn_physical_off()
-        prog_back.integrate(t_back, self._pot)
+        prog_back.integrate(t_back, _track_pot)
         prog_fwd = self._orig_progenitor()
         prog_fwd.turn_physical_off()
-        prog_fwd.integrate(t_fwd, self._pot)
+        prog_fwd.integrate(t_fwd, _track_pot)
         # Combine, skipping the t=0 duplicate
         track_t_grid = numpy.concatenate([t_back[::-1], t_fwd[1:]])
         track_prog_cart = numpy.column_stack(