diff --git a/CHANGES.rst b/CHANGES.rst index 8bc8d9d98..be8f760cf 100644 --- a/CHANGES.rst +++ b/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 ----------- diff --git a/docs/dynamics/mockstreams.rst b/docs/dynamics/mockstreams.rst index 67a3989eb..cfcbf4d66 100644 --- a/docs/dynamics/mockstreams.rst +++ b/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) diff --git a/docs/tutorials/mock-stream-heliocentric.rst b/docs/tutorials/mock-stream-heliocentric.rst index a6ea0d12c..b830c2b07 100644 --- a/docs/tutorials/mock-stream-heliocentric.rst +++ b/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) diff --git a/gala/dynamics/mockstream/df.pxd b/gala/dynamics/mockstream/df.pxd index b96298714..f189bb9f6 100644 --- a/gala/dynamics/mockstream/df.pxd +++ b/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, diff --git a/gala/dynamics/mockstream/df.pyx b/gala/dynamics/mockstream/df.pyx index d092dd3ab..97f04d2a0 100644 --- a/gala/dynamics/mockstream/df.pyx +++ b/gala/dynamics/mockstream/df.pyx @@ -326,6 +326,8 @@ 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) @@ -333,6 +335,31 @@ cdef class FardalStreamDF(BaseStreamDF): 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 = ((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] diff --git a/gala/dynamics/mockstream/tests/test_df.py b/gala/dynamics/mockstream/tests/test_df.py index ac7fe1dfc..5b7751f2b 100644 --- a/gala/dynamics/mockstream/tests/test_df.py +++ b/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 + ) diff --git a/gala/dynamics/mockstream/tests/test_mockstream.py b/gala/dynamics/mockstream/tests/test_mockstream.py index 18c1095c5..f1d122be1 100644 --- a/gala/dynamics/mockstream/tests/test_mockstream.py +++ b/gala/dynamics/mockstream/tests/test_mockstream.py @@ -1,29 +1,35 @@ -import os import itertools +import os # Third-party import astropy.units as u import numpy as np import pytest +from gala.tests.optional_deps import HAS_H5PY + +from ....dynamics import Orbit, PhaseSpacePosition + # Custom -from ....potential import (Hamiltonian, NFWPotential, HernquistPotential, - ConstantRotatingFrame) -from ....dynamics import PhaseSpacePosition, Orbit +from ....potential import ( + ConstantRotatingFrame, + Hamiltonian, + HernquistPotential, + NFWPotential, +) from ....units import galactic from ...nbody import DirectNBody -from ..mockstream_generator import MockStreamGenerator from ..df import FardalStreamDF -from gala.tests.optional_deps import HAS_H5PY +from ..mockstream_generator import MockStreamGenerator def test_init(): - w0 = PhaseSpacePosition(pos=[15., 0., 0]*u.kpc, - vel=[0, 0, 0.13]*u.kpc/u.Myr) - potential = NFWPotential.from_circular_velocity(v_c=0.2, r_s=20., - units=galactic) + w0 = PhaseSpacePosition( + pos=[15.0, 0.0, 0] * u.kpc, vel=[0, 0, 0.13] * u.kpc / u.Myr + ) + potential = NFWPotential.from_circular_velocity(v_c=0.2, r_s=20.0, units=galactic) H = Hamiltonian(potential) - df = FardalStreamDF() + df = FardalStreamDF(gala_modified=True) with pytest.raises(TypeError): MockStreamGenerator(df="some df", hamiltonian=H) @@ -32,106 +38,115 @@ def test_init(): MockStreamGenerator(df=df, hamiltonian=H, progenitor_potential="stuff") # Test validating the input nbody - nbody_w0 = PhaseSpacePosition(pos=[25., 0., 0]*u.kpc, - vel=[0, 0, 0.13]*u.kpc/u.Myr) - potential2 = NFWPotential.from_circular_velocity(v_c=0.2, r_s=25., - units=galactic) - nbody = DirectNBody(w0=nbody_w0, external_potential=potential2, - particle_potentials=[None]) + nbody_w0 = PhaseSpacePosition( + pos=[25.0, 0.0, 0] * u.kpc, vel=[0, 0, 0.13] * u.kpc / u.Myr + ) + potential2 = NFWPotential.from_circular_velocity(v_c=0.2, r_s=25.0, units=galactic) + nbody = DirectNBody( + w0=nbody_w0, external_potential=potential2, particle_potentials=[None] + ) gen = MockStreamGenerator(df=df, hamiltonian=H) with pytest.raises(ValueError): gen._get_nbody(w0, nbody) - frame2 = ConstantRotatingFrame([0, 0, 25.]*u.km/u.s/u.kpc, units=galactic) - nbody = DirectNBody(w0=nbody_w0, external_potential=potential, frame=frame2, - particle_potentials=[None]) + frame2 = ConstantRotatingFrame([0, 0, 25.0] * u.km / u.s / u.kpc, units=galactic) + nbody = DirectNBody( + w0=nbody_w0, + external_potential=potential, + frame=frame2, + particle_potentials=[None], + ) with pytest.raises(ValueError): gen._get_nbody(w0, nbody) # we expect success! - nbody = DirectNBody(w0=nbody_w0, external_potential=potential, - particle_potentials=[None]) + nbody = DirectNBody( + w0=nbody_w0, external_potential=potential, particle_potentials=[None] + ) new_nbody = gen._get_nbody(w0, nbody) # noqa def test_run(): - potential = NFWPotential.from_circular_velocity(v_c=0.2, r_s=20., - units=galactic) + potential = NFWPotential.from_circular_velocity(v_c=0.2, r_s=20.0, units=galactic) H = Hamiltonian(potential) - w0 = PhaseSpacePosition(pos=[15., 0., 0]*u.kpc, - vel=[0, 0, 0.13]*u.kpc/u.Myr) + w0 = PhaseSpacePosition( + pos=[15.0, 0.0, 0] * u.kpc, vel=[0, 0, 0.13] * u.kpc / u.Myr + ) mass = 2.5e4 * u.Msun - prog_pot = HernquistPotential(mass, 4*u.pc, units=galactic) + prog_pot = HernquistPotential(mass, 4 * u.pc, units=galactic) # The basic run: - df = FardalStreamDF() + df = FardalStreamDF(gala_modified=True) gen = MockStreamGenerator(df=df, hamiltonian=H) - stream1, _ = gen.run(w0, mass, dt=-1., n_steps=100) + stream1, _ = gen.run(w0, mass, dt=-1.0, n_steps=100) # Expected errors: with pytest.raises(TypeError): - gen.run(w0, mass.value, dt=-1., n_steps=100) + gen.run(w0, mass.value, dt=-1.0, n_steps=100) # With self-gravity - gen = MockStreamGenerator(df=df, hamiltonian=H, - progenitor_potential=prog_pot) - stream2, _ = gen.run(w0, mass, dt=-1., n_steps=100) + gen = MockStreamGenerator(df=df, hamiltonian=H, progenitor_potential=prog_pot) + stream2, _ = gen.run(w0, mass, dt=-1.0, n_steps=100) assert not u.allclose(stream1.xyz, stream2.xyz) # Skipping release steps: gen = MockStreamGenerator(df=df, hamiltonian=H) - stream3, _ = gen.run(w0, mass, dt=-1., n_steps=100, - release_every=4, n_particles=4) - assert stream3.shape == ((100//4 + 1) * 4 * 2, ) + stream3, _ = gen.run(w0, mass, dt=-1.0, n_steps=100, release_every=4, n_particles=4) + assert stream3.shape == ((100 // 4 + 1) * 4 * 2,) # Custom n_particles: gen = MockStreamGenerator(df=df, hamiltonian=H) n_particles = np.random.randint(0, 4, size=101) - stream3, _ = gen.run(w0, mass, dt=-1., n_steps=100, - release_every=1, n_particles=n_particles) + stream3, _ = gen.run( + w0, mass, dt=-1.0, n_steps=100, release_every=1, n_particles=n_particles + ) assert stream3.shape[0] == 2 * n_particles.sum() # TODO: add nbody test @pytest.mark.parametrize( - 'dt, nsteps, output_every, release_every, n_particles, trail', - list(itertools.product([1, -1], [16, 17], - [1, 2], [1, 4], [1, 4], [True, False]))) -@pytest.mark.skipif(not HAS_H5PY, reason='h5py required for this test') -def test_animate(tmpdir, dt, nsteps, output_every, release_every, - n_particles, trail): + "dt, nsteps, output_every, release_every, n_particles, trail", + list(itertools.product([1, -1], [16, 17], [1, 2], [1, 4], [1, 4], [True, False])), +) +@pytest.mark.skipif(not HAS_H5PY, reason="h5py required for this test") +def test_animate(tmpdir, dt, nsteps, output_every, release_every, n_particles, trail): import h5py - potential = NFWPotential.from_circular_velocity(v_c=0.2, r_s=20., - units=galactic) + potential = NFWPotential.from_circular_velocity(v_c=0.2, r_s=20.0, units=galactic) H = Hamiltonian(potential) - w0 = PhaseSpacePosition(pos=[15., 0., 0]*u.kpc, - vel=[0, 0, 0.13]*u.kpc/u.Myr) + w0 = PhaseSpacePosition( + pos=[15.0, 0.0, 0] * u.kpc, vel=[0, 0, 0.13] * u.kpc / u.Myr + ) mass = 2.5e4 * u.Msun # The basic run: - df = FardalStreamDF(trail=trail) + df = FardalStreamDF(gala_modified=True, trail=trail) gen = MockStreamGenerator(df=df, hamiltonian=H) filename = os.path.join(str(tmpdir), "test.hdf5") - stream, _ = gen.run(w0, mass, dt=dt, n_steps=nsteps, - release_every=release_every, - n_particles=n_particles, - output_every=output_every, output_filename=filename, - overwrite=True) - - with h5py.File(filename, mode='r') as f: - stream_orbits = Orbit.from_hdf5(f['stream']) - nbody_orbits = Orbit.from_hdf5(f['nbody']) + stream, _ = gen.run( + w0, + mass, + dt=dt, + n_steps=nsteps, + release_every=release_every, + n_particles=n_particles, + output_every=output_every, + output_filename=filename, + overwrite=True, + ) + + with h5py.File(filename, mode="r") as f: + stream_orbits = Orbit.from_hdf5(f["stream"]) + nbody_orbits = Orbit.from_hdf5(f["nbody"]) noutput_times = 1 + nsteps // output_every if nsteps % output_every != 0: noutput_times += 1 tail_n_particles = (1 + int(trail)) * n_particles - expected_shape = (noutput_times, - tail_n_particles * (nsteps // release_every + 1)) + expected_shape = (noutput_times, tail_n_particles * (nsteps // release_every + 1)) assert stream_orbits.shape == expected_shape assert np.isfinite(stream_orbits[:, 0].xyz).all() diff --git a/gala/dynamics/mockstream/tests/test_mockstream_class.py b/gala/dynamics/mockstream/tests/test_mockstream_class.py index 77ed807fa..7b9b3e42d 100644 --- a/gala/dynamics/mockstream/tests/test_mockstream_class.py +++ b/gala/dynamics/mockstream/tests/test_mockstream_class.py @@ -13,9 +13,9 @@ def test_init(): vxyz = np.random.random(size=(3, 100)) * u.km / u.s t1 = np.random.random(size=100) * u.Myr - lead_trail = np.empty(100, dtype='U1') - lead_trail[::2] = 't' - lead_trail[1::2] = 'l' + lead_trail = np.empty(100, dtype="U1") + lead_trail[::2] = "t" + lead_trail[1::2] = "l" stream = MockStream(xyz, vxyz) stream = MockStream(xyz, vxyz, release_time=t1) @@ -31,23 +31,23 @@ def test_init(): def test_one_burst(): # Regression test: Tests a bug found by Helmer when putting all particles at # one timestep - import gala.potential as gp import gala.dynamics as gd + import gala.potential as gp from gala.dynamics import mockstream as ms from gala.units import galactic # NFW MW with v_c = 232.8 km/s @ r = 8.2 kpc - pot = gp.NFWPotential.from_circular_velocity(v_c=232.8*u.km/u.s, - r_s=15*u.kpc, - r_ref=8.2*u.kpc, - units=galactic) + pot = gp.NFWPotential.from_circular_velocity( + v_c=232.8 * u.km / u.s, r_s=15 * u.kpc, r_ref=8.2 * u.kpc, units=galactic + ) H = gp.Hamiltonian(pot) - prog_w0 = gd.PhaseSpacePosition(pos=[10, 0, 0.] * u.kpc, - vel=[0, 10, 0.] * u.km/u.s) + prog_w0 = gd.PhaseSpacePosition( + pos=[10, 0, 0.0] * u.kpc, vel=[0, 10, 0.0] * u.km / u.s + ) - dt = 1*u.Myr + dt = 1 * u.Myr nsteps = 100 orbit = H.integrate_orbit(prog_w0, dt=dt, n_steps=nsteps) @@ -57,16 +57,58 @@ def test_one_burst(): argmin = r[0:150].argmin() n_array[argmin] = 1000 - df = ms.FardalStreamDF() + df = ms.FardalStreamDF(gala_modified=True) - dt = 1*u.Myr - prog_mass = 2.5E4 * u.Msun - prog_pot = gp.PlummerPotential(m=prog_mass, b=4*u.pc, units=galactic) + dt = 1 * u.Myr + prog_mass = 2.5e4 * u.Msun + prog_pot = gp.PlummerPotential(m=prog_mass, b=4 * u.pc, units=galactic) gen = ms.MockStreamGenerator(df, H, progenitor_potential=prog_pot) stream, prog = gen.run( - prog_w0, prog_mass, - n_particles=n_array, - dt=dt, - n_steps=nsteps, progress=False) + prog_w0, prog_mass, n_particles=n_array, dt=dt, n_steps=nsteps, progress=False + ) + + +def test_Fardal_vs_GalaModified(): + """ + Regression test: Check that one can actually use the original Fardal parameter + values, and that makes a different stream than the Gala-modified values: + https://github.com/adrn/gala/pull/358 + """ + import gala.dynamics as gd + import gala.potential as gp + from gala.dynamics import mockstream as ms + from gala.units import galactic + + # NFW MW with v_c = 232.8 km/s @ r = 8.2 kpc + pot = gp.NFWPotential.from_circular_velocity( + v_c=232.8 * u.km / u.s, r_s=15 * u.kpc, r_ref=8.2 * u.kpc, units=galactic + ) + + H = gp.Hamiltonian(pot) + + prog_w0 = gd.PhaseSpacePosition( + pos=[10, 0, 0.0] * u.kpc, vel=[0, 300, 20.0] * u.km / u.s + ) + + with pytest.warns(DeprecationWarning, match="Fardal"): + ms.FardalStreamDF() + + df_false = ms.FardalStreamDF( + gala_modified=False, random_state=np.random.default_rng(seed=42) + ) + df_true = ms.FardalStreamDF( + gala_modified=True, random_state=np.random.default_rng(seed=42) + ) + + gen_false = ms.MockStreamGenerator(df_false, H) + gen_true = ms.MockStreamGenerator(df_true, H) + + prog_mass = 2.5e4 * u.Msun + stream_false, _ = gen_false.run( + prog_w0, prog_mass, dt=1, n_steps=128, progress=False + ) + stream_true, _ = gen_true.run(prog_w0, prog_mass, dt=1, n_steps=128, progress=False) + + assert not u.allclose(stream_false.xyz, stream_true.xyz) diff --git a/gala/util.py b/gala/util.py index fea7e0b95..8d637b073 100644 --- a/gala/util.py +++ b/gala/util.py @@ -6,7 +6,7 @@ # Third-party import numpy as np -__all__ = ['rolling_window', 'atleast_2d', 'assert_angles_allclose'] +__all__ = ["rolling_window", "atleast_2d", "assert_angles_allclose"] class ImmutableDict(Mapping): @@ -44,6 +44,7 @@ def __str__(self): def copy(self): import copy + return copy.deepcopy(self._dict) @@ -152,7 +153,7 @@ def atleast_2d(*arys, **kwargs): [array([[1]]), array([[1, 2]]), array([[1, 2]])] """ - insert_axis = kwargs.pop('insert_axis', 0) + insert_axis = kwargs.pop("insert_axis", 0) slc = [slice(None)] * 2 slc[insert_axis] = None slc = tuple(slc) @@ -177,13 +178,12 @@ def assert_angles_allclose(x, y, **kwargs): """ Like numpy's assert_allclose, but for angles (in radians). """ - c2 = (np.sin(x)-np.sin(y))**2 + (np.cos(x)-np.cos(y))**2 - diff = np.arccos((2.0 - c2)/2.0) # a = b = 1 + c2 = (np.sin(x) - np.sin(y)) ** 2 + (np.cos(x) - np.cos(y)) ** 2 + diff = np.arccos((2.0 - c2) / 2.0) # a = b = 1 assert np.allclose(diff, 0.0, **kwargs) -class GalaDeprecationWarning(Warning): +class GalaDeprecationWarning(DeprecationWarning): """ - A warning class to indicate a deprecated feature. Use this over the built-in - DeprecationWarning because those are silenced by default. + A warning class to indicate a deprecated feature. """