Merge pull request #358 from jl3937/0131

Bug fix in FardalStreamDF.

CHANGES.rst

@@ -8,6 +8,11 @@ New Features
 Bug fixes
 ---------
 
+- Fixed the parameter values in the ``FardalStreamDF`` class to be consistent with
+  the values used in Fardal et al. (2015). Added an option ``gala_modified`` to the
+  class to enable using the new (correct) parameter values, but the default will
+  continue to use the Gala modified values (for backwards compatibility).
+
 
 API changes
 -----------

docs/dynamics/mockstreams.rst

@@ -77,9 +77,20 @@ argument values. We will also need to specify the progenitor mass, which is
 passed in to any ``StreamDF`` and is used to scale the particle release
 distribution::
 
-    >>> df = ms.FardalStreamDF()
+    >>> df = ms.FardalStreamDF(gala_modified=True)
     >>> prog_mass = 2.5E4 * u.Msun
 
+.. warning::
+
+    The parameter values of the FardalStreamDF have been updated (fixed) in v1.9 to
+    match the parameter values in the final published version of [fardal15]_. For now,
+    this class uses the Gala modified parameter values that have been adopted over the
+    last several years in Gala. In the future, the default behavior of this class will
+    use the [fardal15]_ parameter values instead, breaking backwards compatibility for
+    mock stream simulations. To use the [fardal15]_ parameters now, set
+    ``gala_modified=False``. To continue to use the Gala modified parameter values, set
+    ``gala_modified=True``.
+
 The final step before actually generating the stream is to create a
 `~gala.dynamics.mockstream.MockStreamGenerator` instance, which we will use to
 actually generate the stream. This takes the ``StreamDF`` and the external
@@ -115,7 +126,7 @@ Let's plot the stream::
     prog_w0 = gd.PhaseSpacePosition(pos=[10, 0, 0.] * u.kpc,
                                     vel=[0, 170, 0.] * u.km/u.s)
 
-    df = ms.FardalStreamDF()
+    df = ms.FardalStreamDF(gala_modified=True)
     prog_mass = 2.5E4 * u.Msun
 
     gen = ms.MockStreamGenerator(df, H)

docs/tutorials/mock-stream-heliocentric.rst

@@ -118,7 +118,7 @@ Here the negative timestep tells the stream generator to first integrate the orb
     pal5_mass = 2.5e4 * u.Msun
     pal5_pot = gp.PlummerPotential(m=pal5_mass, b=4*u.pc, units=galactic)
 
-    df = ms.FardalStreamDF()
+    df = ms.FardalStreamDF(gala_modified=True)
     gen_pal5 = ms.MockStreamGenerator(df, pot, progenitor_potential=pal5_pot)
     pal5_stream, _ = gen_pal5.run(pal5_w0, pal5_mass,
                                  dt=-1 * u.Myr, n_steps=4000)

gala/dynamics/mockstream/df.pxd

@@ -8,6 +8,9 @@ cdef class BaseStreamDF:
     cdef double _trail
     cdef public object random_state
 
+    # TODO: used only in the FardalStreamDF class
+    cdef int _gala_modified
+
     cdef void get_rj_vj_R(self, CPotential *cpotential, double G,
                           double *prog_x, double *prog_v,
                           double prog_m, double t,

gala/dynamics/mockstream/df.pyx

@@ -326,13 +326,40 @@ cdef class FardalStreamDF(BaseStreamDF):
 
     Parameters
     ----------
+    gala_modified : bool (optional)
+        If True, use the modified version of the Fardal method parameters used in Gala. If you would like to use the exact parameters from Fardal+2015, set this to False. Default: True.
     lead : bool (optional)
         Generate a leading tail. Default: True.
     trail : bool (optional)
         Generate a trailing tail. Default: True.
     random_state : `~numpy.random.RandomState` (optional)
         To control random number generation.
     """
+    def __init__(
+        self, gala_modified=None, lead=True, trail=True, random_state=None
+    ):
+        super().__init__(lead=lead, trail=trail, random_state=random_state)
+
+        if gala_modified is None:
+            from gala.util import GalaDeprecationWarning
+            import warnings
+            msg = (
+                "The parameter values of the FardalStreamDF have been updated (fixed) "
+                "to match the parameter values in the final published version of "
+                "Fardal+2015. For now, this class uses the Gala modified parameter "
+                "values that have been adopted over the last several years in Gala. "
+                "In the future, the default behavior of this class will use the "
+                "Fardal+2015 parameter values instead, breaking backwards "
+                "compatibility for mock stream simulations. To use the Fardal+2015 "
+                "parameters now, set gala_modified=False. To continue to use the Gala "
+                "modified parameter values, set gala_modified=True."
+            )
+            warnings.warn(msg, GalaDeprecationWarning)
+            gala_modified = True
+
+        self._gala_modified = int(gala_modified)
+
+
     cpdef _sample(self, potential,
                   double[:, ::1] prog_x, double[:, ::1] prog_v,
                   double[::1] prog_t, double[::1] prog_m, int[::1] nparticles):
@@ -360,14 +387,19 @@ cdef class FardalStreamDF(BaseStreamDF):
             CPotential cpotential = (<CPotentialWrapper>(potential.c_instance)).cpotential
             double G = potential.G
 
+        # TODO: support computing this, which requires knowing the peri/apo and values
+        # of Om**2 - d2Phi/dr2 at those points...
+        # kvt_fardal = min(0.15 * self.f_t**2 * Racc**(2/3), 0.4)
+        kvt_fardal = 0.4
+
         k_mean[0] = 2. # R
-        k_disp[0] = 0.5
+        k_disp[0] = 0.5 if self._gala_modified else 0.4
 
         k_mean[2] = 0. # z
         k_disp[2] = 0.5
 
         k_mean[4] = 0.3 # vt
-        k_disp[4] = 0.5
+        k_disp[4] = 0.5 if self._gala_modified else kvt_fardal
 
         k_mean[5] = 0. # vz
         k_disp[5] = 0.5
@@ -387,7 +419,9 @@ cdef class FardalStreamDF(BaseStreamDF):
                     kx = self.random_state.normal(k_mean[0], k_disp[0])
                     tmp_x[0] = kx * rj
                     tmp_x[2] = self.random_state.normal(k_mean[2], k_disp[2]) * rj
-                    tmp_v[1] = kx * self.random_state.normal(k_mean[4], k_disp[4]) * vj
+                    tmp_v[1] = self.random_state.normal(k_mean[4], k_disp[4]) * vj
+                    if self._gala_modified:  # for backwards compatibility
+                        tmp_v[1] *= kx
                     tmp_v[2] = self.random_state.normal(k_mean[5], k_disp[5]) * vj
                     particle_t1[j+k] = prog_t[i]
 
@@ -405,7 +439,9 @@ cdef class FardalStreamDF(BaseStreamDF):
                     kx = self.random_state.normal(k_mean[0], k_disp[0])
                     tmp_x[0] = kx * -rj
                     tmp_x[2] = self.random_state.normal(k_mean[2], k_disp[2]) * -rj
-                    tmp_v[1] = kx * self.random_state.normal(k_mean[4], k_disp[4]) * -vj
+                    tmp_v[1] = self.random_state.normal(k_mean[4], k_disp[4]) * -vj
+                    if self._gala_modified:  # for backwards compatibility
+                        tmp_v[1] *= kx
                     tmp_v[2] = self.random_state.normal(k_mean[5], k_disp[5]) * -vj
                     particle_t1[j+k] = prog_t[i]
 

gala/dynamics/mockstream/tests/test_df.py

@@ -5,52 +5,57 @@
 
 # Custom
 from ....integrate import DOPRI853Integrator
-from ....potential import (Hamiltonian, HernquistPotential, MilkyWayPotential,
-                           ConstantRotatingFrame)
+from ....potential import (
+    ConstantRotatingFrame,
+    Hamiltonian,
+    HernquistPotential,
+    MilkyWayPotential,
+)
 from ....units import galactic
 from ...core import PhaseSpacePosition
 
 # Project
-from ..df import StreaklineStreamDF, FardalStreamDF, LagrangeCloudStreamDF
-
+from ..df import FardalStreamDF, LagrangeCloudStreamDF, StreaklineStreamDF
 
 _DF_CLASSES = [StreaklineStreamDF, FardalStreamDF, LagrangeCloudStreamDF]
-_DF_KWARGS = [{}, {}, {'v_disp': 1*u.km/u.s}]
-_TEST_POTENTIALS = [HernquistPotential(m=1e12, c=5, units=galactic),
-                    MilkyWayPotential()]
+_DF_KWARGS = [{}, {"gala_modified": True}, {"v_disp": 1 * u.km / u.s}]
+_TEST_POTENTIALS = [
+    HernquistPotential(m=1e12, c=5, units=galactic),
+    MilkyWayPotential(),
+]
 
 
-@pytest.mark.parametrize('DF, DF_kwargs', zip(_DF_CLASSES, _DF_KWARGS))
-@pytest.mark.parametrize('pot', _TEST_POTENTIALS)
+@pytest.mark.parametrize("DF, DF_kwargs", zip(_DF_CLASSES, _DF_KWARGS))
+@pytest.mark.parametrize("pot", _TEST_POTENTIALS)
 def test_init_sample(DF, DF_kwargs, pot):
     H = Hamiltonian(pot)
 
-    orbit = H.integrate_orbit([10., 0, 0, 0, 0.2, 0], dt=1., n_steps=100)
+    orbit = H.integrate_orbit([10.0, 0, 0, 0, 0.2, 0], dt=1.0, n_steps=100)
     n_times = len(orbit.t)
 
     # Different ways to initialize successfully:
     df = DF(**DF_kwargs)
-    o = df.sample(orbit, 1e4*u.Msun)
+    o = df.sample(orbit, 1e4 * u.Msun)
     assert len(o.x) == 2 * n_times
 
     df = DF(lead=False, **DF_kwargs)
-    o = df.sample(orbit, 1e4*u.Msun)
+    o = df.sample(orbit, 1e4 * u.Msun)
     assert len(o.x) == n_times
 
     df = DF(trail=False, **DF_kwargs)
-    o = df.sample(orbit, 1e4*u.Msun)
+    o = df.sample(orbit, 1e4 * u.Msun)
     assert len(o.x) == n_times
 
     df1 = DF(random_state=np.random.RandomState(42), **DF_kwargs)
-    o1 = df1.sample(orbit, 1e4*u.Msun)
+    o1 = df1.sample(orbit, 1e4 * u.Msun)
     df2 = DF(random_state=np.random.RandomState(42), **DF_kwargs)
-    o2 = df2.sample(orbit, 1e4*u.Msun)
+    o2 = df2.sample(orbit, 1e4 * u.Msun)
     assert u.allclose(o1.xyz, o2.xyz)
     assert u.allclose(o1.v_xyz, o2.v_xyz)
     assert len(o1.x) == 2 * n_times
 
 
-@pytest.mark.parametrize('DF, DF_kwargs', zip(_DF_CLASSES, _DF_KWARGS))
+@pytest.mark.parametrize("DF, DF_kwargs", zip(_DF_CLASSES, _DF_KWARGS))
 def test_expected_failure(DF, DF_kwargs):
 
     # Expected failure:
@@ -62,33 +67,31 @@ def test_rotating_frame():
     DF = _DF_CLASSES[0]
     H_static = Hamiltonian(_TEST_POTENTIALS[0])
 
-    w0 = PhaseSpacePosition(pos=[10., 0, 0]*u.kpc,
-                            vel=[0, 220, 0.]*u.km/u.s,
-                            frame=H_static.frame)
-    int_kwargs = dict(w0=w0, dt=1, n_steps=100,
-                      Integrator=DOPRI853Integrator)
+    w0 = PhaseSpacePosition(
+        pos=[10.0, 0, 0] * u.kpc, vel=[0, 220, 0.0] * u.km / u.s, frame=H_static.frame
+    )
+    int_kwargs = dict(w0=w0, dt=1, n_steps=100, Integrator=DOPRI853Integrator)
 
     orbit_static = H_static.integrate_orbit(**int_kwargs)
 
-    rframe = ConstantRotatingFrame([0, 0, -40] * u.km/u.s/u.kpc,
-                                   units=galactic)
-    H_rotating = Hamiltonian(_TEST_POTENTIALS[0],
-                             frame=rframe)
+    rframe = ConstantRotatingFrame([0, 0, -40] * u.km / u.s / u.kpc, units=galactic)
+    H_rotating = Hamiltonian(_TEST_POTENTIALS[0], frame=rframe)
     orbit_rotating = H_rotating.integrate_orbit(**int_kwargs)
 
     _o = orbit_rotating.to_frame(H_static.frame)
-    assert u.allclose(_o.xyz, orbit_static.xyz, atol=1e-13*u.kpc)
-    assert u.allclose(_o.v_xyz, orbit_static.v_xyz, atol=1e-13*u.km/u.s)
+    assert u.allclose(_o.xyz, orbit_static.xyz, atol=1e-13 * u.kpc)
+    assert u.allclose(_o.v_xyz, orbit_static.v_xyz, atol=1e-13 * u.km / u.s)
 
     df_static = DF(trail=False)
-    xvt_static = df_static.sample(orbit_static, 1e6*u.Msun)
+    xvt_static = df_static.sample(orbit_static, 1e6 * u.Msun)
 
     df_rotating = DF(trail=False)
-    xvt_rotating = df_rotating.sample(orbit_rotating, 1e6*u.Msun)
-    xvt_rotating_static = xvt_rotating.to_frame(H_static.frame,
-                                                t=xvt_rotating.release_time)
-
-    assert u.allclose(xvt_static.xyz, xvt_rotating_static.xyz,
-                      atol=1e-9*u.kpc)
-    assert u.allclose(xvt_static.v_xyz, xvt_rotating_static.v_xyz,
-                      atol=1e-9*u.kpc/u.Myr)
+    xvt_rotating = df_rotating.sample(orbit_rotating, 1e6 * u.Msun)
+    xvt_rotating_static = xvt_rotating.to_frame(
+        H_static.frame, t=xvt_rotating.release_time
+    )
+
+    assert u.allclose(xvt_static.xyz, xvt_rotating_static.xyz, atol=1e-9 * u.kpc)
+    assert u.allclose(
+        xvt_static.v_xyz, xvt_rotating_static.v_xyz, atol=1e-9 * u.kpc / u.Myr
+    )