Skip to content

Initial conditions

.. ######################################################################## .. ######################################################################## .. # Copyright (c) 2013,2014 Svetlin Tassev .. # Princeton University,Harvard University .. # .. # This file is part of pyCOLA. .. # .. # pyCOLA is free software: you can redistribute it and/or modify .. # it under the terms of the GNU General Public License as published by .. # the Free Software Foundation, either version 3 of the License, or .. # (at your option) any later version. .. # .. # pyCOLA is distributed in the hope that it will be useful, .. # but WITHOUT ANY WARRANTY; without even the implied warranty of .. # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the .. # GNU General Public License for more details. .. # .. # You should have received a copy of the GNU General Public License .. # along with pyCOLA. If not, see http://www.gnu.org/licenses/. .. # .. ######################################################################## .. ########################################################################

ic_2lpt(cellsize, sx, sy, sz, sx_zoom=None, sy_zoom=None, sz_zoom=None, boxsize=100.0, ngrid_x_lpt=128, ngrid_y_lpt=128, ngrid_z_lpt=128, cellsize_zoom=0, offset_zoom=None, BBox_in=None, growth_2pt_calc=0.05, with_4pt_rule=False, factor_4pt=2.0)

:math:\vspace{-1mm}

Given a set of displacements calculated in the ZA at redshift zero, find the corresponding second order displacement. Works with a single grid of particles, as well as with one refined subvolume.

Arguments:

  • cellsize -- a float. The inter-particle spacing in Lagrangian space.

  • sx,sy,sz -- 3-dim NumPy arrays containing the components of the particle displacements today as calculated in the ZA. These particles should cover the whole box. If a refined subvolume is provided, the crude particles which reside inside that subvolume are discarded and replaced with the fine particles.

  • sx_zoom,sy_zoom,sz_zoom -- 3-dim NumPy arrays containing the components of the particle ZA displacements today for a refined subvolume (default: None).

  • boxsize -- a float (default: 100.0). Gives the size of the simulation box in Mpc/h.

  • ngrid_x_lpt,ngrid_y_lpt,ngrid_z_lpt -- integers (default: 128). Provide the size of the PM grid, which the algorithm uses to calculate the 2LPT displacements.

  • cellsize_zoom -- a float (default: 0). The inter-particle spacing in Lagrangian space for the refined subvolume, if such is provided. If not, cellsize_zoom must be set to zero (default), as that is used as a check for the presence of that subvolume.

  • offset_zoom -- a 3-vector of floats (default: None). Gives the physical coordinates of the origin of the refinement region relative to the the origin of the full box.

  • BBox_in -- a 3x2 array of integers (default: None). It has the form [[i0,i1],[j0,j1],[k0,k1]], which gives the bounding box for the refinement region in units of the crude particles Lagrangian index. Thus, the particles with displacements sx|sy|sz[i0:i1,j0:j1,k0:k1] are replaced with the fine particles with displacements sx_zoom|sy_zoom|sz_zoom.

  • growth_2pt_calc -- a float (default: 0.05). The linear growth factor used internally in the 2LPT calculation. A value of 0.05 gives excellent cross-correlation between the 2LPT field returned by this function, and the 2LPT returned using the usual fft tricks. Yet, some irrelevant short-scale noise is present, which one may decide to filter out. That noise is most probably due to lack of force accuracy for too low growth_2pt_calc. Experiment with this value as needed.

  • with_4pt_rule -- a boolean (default: False). See :func:ic.ic_2lpt_engine.

  • factor_4pt -- a float (default: 2.0). See :func:ic.ic_2lpt_engine.

Return:

  • If no refined subregion is supplied (indicated by cellsize_zoom=0), then return:

(sx2,sy2,sz2) -- 3-dim NumPy arrays containing the components of the second order particle displacements today as calculated in 2LPT.

  • If a refined subregion is supplied (indicated by cellsize_zoom>0), then return:

(sx2,sy2,sz2,sx2_zoom,sy2_zoom,sz2_zoom)

The first three arrays are as above. The last three give the second order displacements today for the particles in the refined subvolume.

Example:

Generate a realization for the displacement field in the ZA; calculate the corresponding second order displacement field; calculate the rms displacements; then show a projection of one of the components.

>>> from ic import ic_za,ic_2lpt
>>> sx,sy,sz=ic_za('camb_matterpower_z0.dat',npart=128)
Memory allocation done
Plans created
Power spectrum read.
Randoms done.
Nyquists fixed
sx fft ready
sy fft ready
sz fft ready
>>> sx2,sy2,sz2 = ic_2lpt( 100.0/float(sx.shape[0]),sx,sy,sz,
...                        growth_2pt_calc=0.1)
>>> ((sx**2+sy**2+sz**2).mean())**0.5/0.7    # ~10
11.605451188108798
>>> ((sx2**2+sy2**2+sz2**2).mean())**0.5/0.7 # ~2
2.3447627779313525
>>> import matplotlib.pyplot as plt # needs matplotlib to be installed!
>>> plt.imshow(sx.mean(axis=2))
<matplotlib.image.AxesImage object at 0x7fc4603697d0>
>>> plt.show()
>>> plt.imshow(sy2.mean(axis=2))
<matplotlib.image.AxesImage object at 0x7fc4603697d0>
>>> plt.show()

Algorithm:

This function issues a call to :func:ic.ic_2lpt_engine. See the Algorithm section of that function for details.

Source code in pycola3/ic.py
def ic_2lpt(
    cellsize,
    sx,
    sy,
    sz,
    sx_zoom=None,
    sy_zoom=None,
    sz_zoom=None,
    boxsize=100.00,
    ngrid_x_lpt=128,
    ngrid_y_lpt=128,
    ngrid_z_lpt=128,
    cellsize_zoom=0,
    offset_zoom=None,
    BBox_in=None,
    growth_2pt_calc=0.05,
    with_4pt_rule=False,
    factor_4pt=2.0,
):
    """
    :math:`\\vspace{-1mm}`


    Given a set of displacements calculated in the ZA at redshift
    zero, find the corresponding second order displacement. Works
    with a single grid of particles, as well as with one refined subvolume.

    **Arguments**:

    * ``cellsize`` --  a float. The inter-particle spacing in Lagrangian space.

    * ``sx,sy,sz`` --  3-dim NumPy arrays containing the
      components of the particle displacements today as calculated in
      the ZA. These particles should cover the whole box. If a refined
      subvolume is provided, the crude particles which reside inside
      that subvolume are discarded and replaced with the fine
      particles.

    * ``sx_zoom,sy_zoom,sz_zoom`` --  3-dim NumPy arrays containing the
      components of the particle
      ZA displacements today for a refined subvolume (default: ``None``).

    * ``boxsize`` -- a float (default: ``100.0``). Gives the size of the
      simulation box in Mpc/h.

    * ``ngrid_x_lpt,ngrid_y_lpt,ngrid_z_lpt`` -- integers
      (default: ``128``). Provide the size of the PM grid, which the algorithm
      uses to calculate the 2LPT displacements.

    * ``cellsize_zoom`` -- a float (default: ``0``). The inter-particle
      spacing in Lagrangian space for the refined subvolume, if such is
      provided. If not, ``cellsize_zoom`` must be set to zero
      (default), as that is used as a check for the presence of that
      subvolume.

    * ``offset_zoom`` -- a 3-vector of floats (default: ``None``). Gives the
      physical coordinates of the origin of the refinement region
      relative to the the origin of the full box.

    * ``BBox_in`` -- a 3x2 array of integers (default: ``None``). It has the
      form ``[[i0,i1],[j0,j1],[k0,k1]]``, which gives the bounding box
      for the refinement region in units of the crude particles
      Lagrangian index. Thus, the particles with displacements
      ``sx|sy|sz[i0:i1,j0:j1,k0:k1]`` are replaced with the fine
      particles with displacements ``sx_zoom|sy_zoom|sz_zoom``.

    * ``growth_2pt_calc`` --  a float (default: ``0.05``). The
      linear growth factor used internally in the 2LPT calculation. A
      value of 0.05 gives excellent cross-correlation between the 2LPT
      field returned by this function, and the 2LPT returned using the
      usual fft tricks. Yet, some irrelevant short-scale noise is
      present, which one may decide to filter out. That noise is most
      probably due to lack of force accuracy for too low
      ``growth_2pt_calc``. Experiment with this value as
      needed.

    * ``with_4pt_rule`` -- a boolean (default: ``False``). See :func:`ic.ic_2lpt_engine`.

    * ``factor_4pt`` -- a float (default: ``2.0``). See :func:`ic.ic_2lpt_engine`.



    **Return**:

    * If no refined subregion is supplied (indicated by
      ``cellsize_zoom=0``), then return:

      ``(sx2,sy2,sz2)`` --  3-dim NumPy
      arrays containing the components of the second order particle
      displacements today as calculated in 2LPT.

    * If a refined subregion is supplied (indicated by
      ``cellsize_zoom>0``), then return:

      ``(sx2,sy2,sz2,sx2_zoom,sy2_zoom,sz2_zoom)``

      The first three arrays are as
      above. The last three give the second order displacements today
      for the particles in the refined subvolume.

    **Example**:

    Generate a realization for the displacement field in the ZA;
    calculate the corresponding second order displacement field; calculate
    the rms displacements; then show a projection of one of the
    components.

        >>> from ic import ic_za,ic_2lpt
        >>> sx,sy,sz=ic_za('camb_matterpower_z0.dat',npart=128)
        Memory allocation done
        Plans created
        Power spectrum read.
        Randoms done.
        Nyquists fixed
        sx fft ready
        sy fft ready
        sz fft ready
        >>> sx2,sy2,sz2 = ic_2lpt( 100.0/float(sx.shape[0]),sx,sy,sz,
        ...                        growth_2pt_calc=0.1)
        >>> ((sx**2+sy**2+sz**2).mean())**0.5/0.7    # ~10
        11.605451188108798
        >>> ((sx2**2+sy2**2+sz2**2).mean())**0.5/0.7 # ~2
        2.3447627779313525
        >>> import matplotlib.pyplot as plt # needs matplotlib to be installed!
        >>> plt.imshow(sx.mean(axis=2))
        <matplotlib.image.AxesImage object at 0x7fc4603697d0>
        >>> plt.show()
        >>> plt.imshow(sy2.mean(axis=2))
        <matplotlib.image.AxesImage object at 0x7fc4603697d0>
        >>> plt.show()

    **Algorithm**:

    This function issues a call to :func:`ic.ic_2lpt_engine`. See the
    Algorithm section of that function for details.

    """

    # from ic import ic_2lpt_engine
    res = ic_2lpt_engine(
        sx,
        sy,
        sz,
        cellsize,
        ngrid_x_lpt,
        ngrid_y_lpt,
        ngrid_z_lpt,
        boxsize / float(ngrid_x_lpt),  # assumes cube
        with_2lpt=False,
        sx2_full=None,
        sy2_full=None,
        sz2_full=None,
        cellsize_zoom=cellsize_zoom,
        BBox_in=BBox_in,
        sx_full_zoom=sx_zoom,
        sy_full_zoom=sy_zoom,
        sz_full_zoom=sz_zoom,
        sx2_full_zoom=None,
        sy2_full_zoom=None,
        sz2_full_zoom=None,
        offset_zoom=offset_zoom,
        growth_2pt_calc=growth_2pt_calc,
    )
    if cellsize_zoom != 0:
        (
            sx_,
            sy_,
            sz_,
            sx2,
            sy2,
            sz2,
            sx_zoom_,
            sy_zoom_,
            sz_zoom_,
            sx2_zoom,
            sy2_zoom,
            sz2_zoom,
        ) = res
    else:
        sx_, sy_, sz_, sx2, sy2, sz2 = res

    del (
        sx_,
        sy_,
        sz_,
    )  # These have higher order corrections unlike the original *_full arrays, which are 'exact'. So, discard.
    if cellsize_zoom != 0:
        del sx_zoom_, sy_zoom_, sz_zoom_
        return sx2, sy2, sz2, sx2_zoom, sy2_zoom, sz2_zoom
    return sx2, sy2, sz2

ic_2lpt_engine(sx_full, sy_full, sz_full, cellsize, ngrid_x, ngrid_y, ngrid_z, gridcellsize, growth_2pt_calc=0.05, with_4pt_rule=False, factor_4pt=2.0, with_2lpt=False, sx2_full=None, sy2_full=None, sz2_full=None, cellsize_zoom=0, BBox_in=None, sx_full_zoom=None, sy_full_zoom=None, sz_full_zoom=None, sx2_full_zoom=None, sy2_full_zoom=None, sz2_full_zoom=None, offset_zoom=None)

:math:\vspace{-1mm}

The same as :func:ic.ic_2lpt above, but calculates the 2LPT displacements for the particles in the COLA volume as generated by same particles displaced according to the 2LPT of the full box. (todo: expand this!) In fact, :func:ic.ic_2lpt works by making a call to this function.

Arguments:

  • sx_full,sy_full,sz_full -- 3-dim NumPy arrays containing the components of the particle displacements today as calculated in the ZA of the full box. These particles should cover the whole box. If a refined subvolume is provided, the crude particles which reside inside that subvolume are discarded and replaced with the fine particles.

  • cellsize -- a float. The inter-particle spacing in Lagrangian space.

  • ngrid_x,ngrid_y,ngrid_z -- integers. Provide the size of the PM grid, which the algorithm uses to calculate the 2LPT displacements.

  • gridcellsize --float. Provide the grid spacing of the PM grid, which the algorithm uses to calculate the 2LPT displacements.

  • growth_2pt_calc -- a float (default: 0.05). The linear growth factor used internally in the 2LPT calculation. A value of 0.05 gives excellent cross-correlation between the 2LPT field returned by this function, and the 2LPT returned using the usual fft tricks for a 100:math:\mathrm{Mpc}/h box. Yet, some irrelevant short-scale noise is present, which one may decide to filter out. That noise is probably due to lack of force accuracy for too low growth_2pt_calc. Experiment with this value as needed.

  • with_4pt_rule -- a boolean (default: False). Whether to use the 4-point force rule to evaluate the ZA and 2LPT displacements in the COLA region. See the Algorithm section below. If set to False, it uses the 2-point force rule.

  • factor_4pt -- a float, different from 1.0 (default: 2.0). Used for the 4-point force rule. See the Algorithm section below.

  • with_2lpt -- a boolean (default: False). Whether the second order displacement field over the full box is provided. One must provide it if the COLA volume is different from the full box. Only if they are the same (as in the case of ic_2lpt()) can one set with_2lpt=False.

  • sx2_full,sy2_full,sz2_full -- 3-dim NumPy float arrays giving the second order displacement field over the full box. Needs with_2lpt=True.

  • The rest of the input is as in :func:ic.ic_2lpt, with all LPT quantities provided for the whole box.

Return:

  • If no refined subregion is supplied (indicated by cellsize_zoom=0), then return:

(sx,sy,sz,sx2,sy2,sz2) -- 3-dim NumPy arrays containing the components of the first and second (s:sub:i\ 2) order particle displacements today as calculated in 2LPT in the COLA volume.

  • If a refined subregion is supplied (indicated by cellsize_zoom>0), then return:

(sx,sy,sz,sx2,sy2,sz2,sx_zoom,sy_zoom,sz_zoom,sx2_zoom,sy2_zoom,sz2_zoom)

The first 6 arrays are as above. The last 6 give the second order displacements today for the particles in the refined subvolume of the COLA volume.

Algorithm:

The first-order and second-order displacements, :math:\bm{s}^{(1)}_{\mathrm{COLA}} and :math:\bm{s}^{(2)}_{\mathrm{COLA}}, in the COLA volume at redshift zero are calculated according to the following 2-pt or 4-pt (denoted by subscript) equations:

.. math:: :nowrap:

\begin{eqnarray} \bm{s}{\mathrm{COLA},\mathrm{2pt}}^{(1)} & = & - \frac{1}{2g} \left[\bm{F}(g,\beta g^2)-\bm{F}(-g,\beta g^2)\right] \ \bm{s}{\mathrm{COLA},\mathrm{2pt}}^{(2)} & = & - \frac{\alpha}{2g^2} \left[\bm{F}(g,\beta g^2)+\bm{F}(-g,\beta g^2)\right] \ \bm{s}{\mathrm{COLA},\mathrm{4pt}}^{(1)} & = & - \frac{1}{2g} \frac{a^2}{a^2-1}\bigg[\bm{F}(g,\beta g^2)-\bm{F}(-g,\beta g^2)-\ & & \quad \quad \quad \quad \quad \quad - \frac{1}{a^3}\bigg(\bm{F}\left(a g,\beta a^2 g^2\right)-\bm{F}\left(-a g,\beta a^2 g^2\right)\bigg)\bigg] \ \bm{s}{\mathrm{COLA},\mathrm{4pt}}^{(2)} & = & - \frac{\alpha}{2g^2}\frac{a^2}{a^2-1}\bigg[\bm{F}(g,\beta g^2)+\bm{F}(-g,\beta g^2)-\ & & \quad \quad \quad \quad \quad \quad - \frac{1}{a^4}\bigg(\bm{F}\left(a g,\beta a^2 g^2\right)+\bm{F}\left(-a g,\beta a^2 g^2\right)\bigg)\bigg] \end{eqnarray}

where:

:math:a= factor_4pt

:math:g= growth_2pt_calc

if with_2lpt then:

  :math:`\beta=1` and :math:`\alpha=(3/10)\Omega_{m}^{1/143}`

else:

  :math:`\beta=0` and :math:`\alpha=(3/7)\Omega_{m}^{1/143}`

The factors of :math:\Omega_{m}^{1/143} (:math:\Omega_m being the matter density today) are needed to rescale the second order displacements to matter domination and are correct to :math:\mathcal{O}(\max(10^{-4},g^3/143)) in :math:\Lambda\mathrm{CDM}. The force :math:\bm{F}(g_1,g_2) is given by:

.. math:: :nowrap:

\begin{eqnarray} \bm{F}(g_1,g_2) = \bm{\nabla}\nabla^{-2}\delta\left[g_1\bm{s}{\mathrm{full}}^{(1)}+g_2\Omega{m}^{-1/143}\bm{s}_{\mathrm{full}}^{(2)}\right] \end{eqnarray}

where :math:\delta[\bm{s}] is the cloud-in-cell fractional overdensity calculated from a grid of particles displaced by the input displacement vector field :math:\bm{s}. Above, :math:\bm{s}_{\mathrm{full}}^{(1)}/\bm{s}_{\mathrm{full}}^{(2)} are the input first/second-order input displacement fields calculated in the full box at redshift zero.

It is important to note that implicitly for each particle at Lagrangian position :math:\bm{q}, the force :math:\bm{F}(g_1,g_2) is evaluated at the corresponding Eulerian position: :math:\bm{q}+g_1\bm{s}_{\mathrm{full}}^{(1)}+g_2\Omega_{m}^{-1/143}\bm{s}_{\mathrm{full}}^{(2)}.

As noted above, with_2lpt=False is only allowed if the COLA volume covers the full box volume. In that case, :math:\bm{s}_{\mathrm{full}}^{(2)} is not needed as input since :math:\beta=0. Instead, the output :math:\bm{s}_{\mathrm{COLA}}^{(2)} can be used as a good approximation to :math:\bm{s}_{\mathrm{full}}^{(2)}. This fact is used in :func:ic.ic_2lpt to calculate :math:\bm{s}_{\mathrm{full}}^{(2)} from :math:\bm{s}_{\mathrm{full}}^{(1)}.

.. note:: If with_4pt_rule=False, then the first/second order displacements receive corrections at third/fourth order. If with_4pt_rule=True, then those corrections are fifth/sixth order. However, when using the 4-point rule instead of the 2-point rule, one must make two more force evaluations at a slightly different growth factor given by growth_2pt_calc*factor_4pt. Since the code is single precision and is using a simple PM grid to evaluate forces, one cannot make factor_4pt and growth_2pt_calc too small due to noise issues. Thus, when comparing the 2-pt and 4-pt rule, we should assume factor_4pt>1. And again due to numerical precision issues, one cannot choose factor_4pt to be too close to one; hence, the default value of 2.0.

Therefore, as the higher order corrections for the 4-pt rule are proportional to powers of growth_2pt_calc*factor_4pt, one may be better off using the 2-pt rule (the default) in this particular implementation. Yet for codes where force accuracy is not an issue, one may consider using the 4-pt rule. Thus, its inclusion in this code is mostly done as an illustration.

Source code in pycola3/ic.py
def ic_2lpt_engine(
    sx_full,
    sy_full,
    sz_full,
    cellsize,
    ngrid_x,
    ngrid_y,
    ngrid_z,
    gridcellsize,
    growth_2pt_calc=0.05,
    with_4pt_rule=False,
    factor_4pt=2.0,
    with_2lpt=False,
    sx2_full=None,
    sy2_full=None,
    sz2_full=None,
    cellsize_zoom=0,
    BBox_in=None,
    sx_full_zoom=None,
    sy_full_zoom=None,
    sz_full_zoom=None,
    sx2_full_zoom=None,
    sy2_full_zoom=None,
    sz2_full_zoom=None,
    offset_zoom=None,
):
    r"""
    :math:`\vspace{-1mm}`

    The same as :func:`ic.ic_2lpt` above, but calculates the 2LPT displacements for the particles in the
    COLA volume as generated by same particles displaced
    according to the 2LPT of the full box. (todo: *expand this!*) In fact, :func:`ic.ic_2lpt` works
    by making a call to this function.

    **Arguments**:

    * ``sx_full,sy_full,sz_full`` --  3-dim NumPy arrays containing the
      components of the particle displacements today as calculated in
      the ZA of the full box. These particles should cover the whole box. If a refined
      subvolume is provided, the crude particles which reside inside
      that subvolume are discarded and replaced with the fine
      particles.

    * ``cellsize`` --  a float. The inter-particle spacing in Lagrangian space.

    * ``ngrid_x,ngrid_y,ngrid_z`` -- integers. Provide the size of the
      PM grid, which the algorithm
      uses to calculate the 2LPT displacements.

    * ``gridcellsize`` --float. Provide the grid spacing of the PM
      grid, which the algorithm
      uses to calculate the 2LPT displacements.

    * ``growth_2pt_calc`` --  a float (default: ``0.05``). The
      linear growth factor used internally in the 2LPT calculation. A
      value of 0.05 gives excellent cross-correlation between the 2LPT
      field returned by this function, and the 2LPT returned using the
      usual fft tricks for a 100:math:`\mathrm{Mpc}/h` box. Yet, some
      irrelevant short-scale noise is present, which one may decide to
      filter out. That noise is probably due to lack of force accuracy
      for too low ``growth_2pt_calc``. Experiment with this value as
      needed.


    * ``with_4pt_rule`` -- a boolean (default: ``False``). Whether to use
      the 4-point force rule to evaluate the ZA and 2LPT displacements
      in the COLA region. See the Algorithm section below. If set to
      False, it uses the 2-point force rule.

    * ``factor_4pt`` -- a float, different from ``1.0`` (default:
      ``2.0``). Used for the 4-point
      force rule. See the Algorithm section below.

    * ``with_2lpt`` -- a boolean (default: ``False``). Whether the second
      order displacement field  over the full box is provided. One must
      provide it if the COLA volume is different from the full box.
      Only if they are the same (as in the case of ``ic_2lpt()``) can
      one set ``with_2lpt=False``.

    * ``sx2_full,sy2_full,sz2_full`` -- 3-dim NumPy float arrays giving the second
      order displacement field  over the full box. Needs ``with_2lpt=True``.

    * The rest of the input is as in :func:`ic.ic_2lpt`, with all LPT
      quantities provided for the whole box.


    **Return**:

    * If no refined subregion is supplied (indicated by
      ``cellsize_zoom=0``), then return:

      ``(sx,sy,sz,sx2,sy2,sz2)`` --  3-dim NumPy
      arrays containing the components of the first and second (``s``:sub:`i`\ ``2``)
      order particle displacements today as calculated in 2LPT in the
      COLA volume.

    * If a refined subregion is supplied (indicated by
      ``cellsize_zoom>0``), then return:

    ``(sx,sy,sz,sx2,sy2,sz2,sx_zoom,sy_zoom,sz_zoom,sx2_zoom,sy2_zoom,sz2_zoom)``

      The first 6 arrays are as
      above. The last 6 give the second order displacements today
      for the particles in the refined subvolume of the COLA volume.

    **Algorithm**:

    The first-order and second-order displacements,
    :math:`\bm{s}^{(1)}_{\mathrm{COLA}}` and
    :math:`\bm{s}^{(2)}_{\mathrm{COLA}}`, in the COLA volume at
    redshift zero are calculated according to the following 2-pt or
    4-pt (denoted by subscript) equations:

    .. math::
      :nowrap:

      \begin{eqnarray}
          \bm{s}_{\mathrm{COLA},\mathrm{2pt}}^{(1)}    & = & - \frac{1}{2g}                      \left[\bm{F}(g,\beta g^2)-\bm{F}(-g,\beta g^2)\right] \\
          \bm{s}_{\mathrm{COLA},\mathrm{2pt}}^{(2)}    & = & - \frac{\alpha}{2g^2}                 \left[\bm{F}(g,\beta g^2)+\bm{F}(-g,\beta g^2)\right] \\
          \bm{s}_{\mathrm{COLA},\mathrm{4pt}}^{(1)}    & = & - \frac{1}{2g}     \frac{a^2}{a^2-1}\bigg[\bm{F}(g,\beta g^2)-\bm{F}(-g,\beta g^2)-\\
                                                     &   & \quad \quad \quad \quad \quad \quad - \frac{1}{a^3}\bigg(\bm{F}\left(a g,\beta a^2 g^2\right)-\bm{F}\left(-a g,\beta a^2 g^2\right)\bigg)\bigg] \\
          \bm{s}_{\mathrm{COLA},\mathrm{4pt}}^{(2)}    & = & - \frac{\alpha}{2g^2}\frac{a^2}{a^2-1}\bigg[\bm{F}(g,\beta g^2)+\bm{F}(-g,\beta g^2)-\\
                                                     &   & \quad \quad \quad \quad \quad \quad - \frac{1}{a^4}\bigg(\bm{F}\left(a g,\beta a^2 g^2\right)+\bm{F}\left(-a g,\beta a^2 g^2\right)\bigg)\bigg]
      \end{eqnarray}

    where:

      :math:`a=` ``factor_4pt``

      :math:`g=` ``growth_2pt_calc``

      if ``with_2lpt`` then:

          :math:`\beta=1` and :math:`\alpha=(3/10)\Omega_{m}^{1/143}`

      else:

          :math:`\beta=0` and :math:`\alpha=(3/7)\Omega_{m}^{1/143}`

    The factors of :math:`\Omega_{m}^{1/143}` (:math:`\Omega_m` being
    the matter density today) are needed to rescale the second order
    displacements to matter domination and are correct to
    :math:`\mathcal{O}(\max(10^{-4},g^3/143))` in
    :math:`\Lambda\mathrm{CDM}`. The force :math:`\bm{F}(g_1,g_2)` is
    given by:

    .. math::
      :nowrap:

      \begin{eqnarray}
          \bm{F}(g_1,g_2) = \bm{\nabla}\nabla^{-2}\delta\left[g_1\bm{s}_{\mathrm{full}}^{(1)}+g_2\Omega_{m}^{-1/143}\bm{s}_{\mathrm{full}}^{(2)}\right]
      \end{eqnarray}

    where :math:`\delta[\bm{s}]` is the cloud-in-cell fractional
    overdensity calculated from a grid of particles displaced by the
    input displacement vector  field :math:`\bm{s}`. Above,
    :math:`\bm{s}_{\mathrm{full}}^{(1)}/\bm{s}_{\mathrm{full}}^{(2)}` are
    the input first/second-order input displacement fields calculated
    in the full box at redshift zero.

    It is important to note that implicitly for each particle at
    Lagrangian position :math:`\bm{q}`, the force
    :math:`\bm{F}(g_1,g_2)` is evaluated at the corresponding Eulerian position:
    :math:`\bm{q}+g_1\bm{s}_{\mathrm{full}}^{(1)}+g_2\Omega_{m}^{-1/143}\bm{s}_{\mathrm{full}}^{(2)}`.

    As noted above, ``with_2lpt=False`` is only allowed if the COLA
    volume covers the full box volume. In that case,
    :math:`\bm{s}_{\mathrm{full}}^{(2)}` is not needed as input since
    :math:`\beta=0`. Instead, the output
    :math:`\bm{s}_{\mathrm{COLA}}^{(2)}` can be used as a good
    approximation to :math:`\bm{s}_{\mathrm{full}}^{(2)}`. This fact
    is used in :func:`ic.ic_2lpt` to calculate
    :math:`\bm{s}_{\mathrm{full}}^{(2)}` from
    :math:`\bm{s}_{\mathrm{full}}^{(1)}`.


    .. note:: If ``with_4pt_rule=False``, then the first/second order
       displacements receive corrections at third/fourth order. If
       ``with_4pt_rule=True``, then those corrections are fifth/sixth
       order. However, when using the 4-point rule instead of the
       2-point rule, one must make two more force evaluations at a
       slightly different growth factor given by
       ``growth_2pt_calc*factor_4pt``. Since the code is single
       precision and is using a simple PM grid to evaluate forces, one
       cannot make ``factor_4pt`` and ``growth_2pt_calc`` too small due
       to noise issues. Thus, when comparing the 2-pt and 4-pt rule, we
       should assume ``factor_4pt>1``. And again due to numerical
       precision issues, one cannot choose ``factor_4pt`` to be too
       close to one; hence, the default value of ``2.0``.

       Therefore, as the higher order corrections for the 4-pt rule are
       proportional to powers of ``growth_2pt_calc*factor_4pt``, one
       may be better off using the 2-pt rule (the default) in this
       particular implementation. Yet for codes where force accuracy is
       not an issue, one may consider using the 4-pt rule. Thus, its
       inclusion in this code is mostly done as an illustration.

    """

    from numpy import float32, float64

    if cellsize_zoom != 0:
        cellsize_zoom = float32(cellsize_zoom)
        offset_zoom = offset_zoom.astype("float32")

    npart_x, npart_y, npart_z = sx_full.shape
    if cellsize_zoom != 0:
        npart_x_zoom, npart_y_zoom, npart_z_zoom = sx_full_zoom.shape

    from numpy import array, zeros

    sx = zeros((npart_x, npart_y, npart_z), dtype="float32")
    sy = zeros((npart_x, npart_y, npart_z), dtype="float32")
    sz = zeros((npart_x, npart_y, npart_z), dtype="float32")

    sx_minus = zeros((npart_x, npart_y, npart_z), dtype="float32")
    sy_minus = zeros((npart_x, npart_y, npart_z), dtype="float32")
    sz_minus = zeros((npart_x, npart_y, npart_z), dtype="float32")

    if cellsize_zoom != 0:
        sx_zoom = zeros((npart_x_zoom, npart_y_zoom, npart_z_zoom), dtype="float32")
        sy_zoom = zeros((npart_x_zoom, npart_y_zoom, npart_z_zoom), dtype="float32")
        sz_zoom = zeros((npart_x_zoom, npart_y_zoom, npart_z_zoom), dtype="float32")

        sx_minus_zoom = zeros(
            (npart_x_zoom, npart_y_zoom, npart_z_zoom), dtype="float32"
        )
        sy_minus_zoom = zeros(
            (npart_x_zoom, npart_y_zoom, npart_z_zoom), dtype="float32"
        )
        sz_minus_zoom = zeros(
            (npart_x_zoom, npart_y_zoom, npart_z_zoom), dtype="float32"
        )

    if with_4pt_rule:
        sx_4pt_minus = zeros((npart_x, npart_y, npart_z), dtype="float32")
        sy_4pt_minus = zeros((npart_x, npart_y, npart_z), dtype="float32")
        sz_4pt_minus = zeros((npart_x, npart_y, npart_z), dtype="float32")

        sx_4pt = zeros((npart_x, npart_y, npart_z), dtype="float32")
        sy_4pt = zeros((npart_x, npart_y, npart_z), dtype="float32")
        sz_4pt = zeros((npart_x, npart_y, npart_z), dtype="float32")

        if cellsize_zoom != 0:
            sx_4pt_minus_zoom = zeros(
                (npart_x_zoom, npart_y_zoom, npart_z_zoom), dtype="float32"
            )
            sy_4pt_minus_zoom = zeros(
                (npart_x_zoom, npart_y_zoom, npart_z_zoom), dtype="float32"
            )
            sz_4pt_minus_zoom = zeros(
                (npart_x_zoom, npart_y_zoom, npart_z_zoom), dtype="float32"
            )

            sx_4pt_zoom = zeros(
                (npart_x_zoom, npart_y_zoom, npart_z_zoom), dtype="float32"
            )
            sy_4pt_zoom = zeros(
                (npart_x_zoom, npart_y_zoom, npart_z_zoom), dtype="float32"
            )
            sz_4pt_zoom = zeros(
                (npart_x_zoom, npart_y_zoom, npart_z_zoom), dtype="float32"
            )
    else:
        sx_4pt_minus = 0.0
        sy_4pt_minus = 0.0
        sz_4pt_minus = 0.0

        sx_4pt = 0.0
        sy_4pt = 0.0
        sz_4pt = 0.0

        if cellsize_zoom != 0:
            sx_4pt_minus_zoom = 0.0
            sy_4pt_minus_zoom = 0.0
            sz_4pt_minus_zoom = 0.0

            sx_4pt_zoom = 0.0
            sy_4pt_zoom = 0.0
            sz_4pt_zoom = 0.0

    from .acceleration import grad_phi
    from .cic import CICDeposit_3
    from .potential import get_phi, initialize_density

    ###
    density, den_k, den_fft, phi_fft = initialize_density(ngrid_x, ngrid_y, ngrid_z)
    density.fill(0.0)
    Om = 0.274
    dd = Om ** (
        1.0 / 143.0
    )  # this is a good enough approximation at early times and is ~0.95
    if with_2lpt:
        cc = 3.0 / 10.0 * dd
        L2 = growth_2pt_calc * growth_2pt_calc / dd
    else:
        cc = 3.0 / 7.0 * dd
        L2 = 0.0

    gridcellsize = float32(gridcellsize)
    growth_2pt_calc = float32(growth_2pt_calc)
    L2 = float32(L2)

    offset = array([0.0, 0.0, 0.0], dtype="float32")
    if cellsize_zoom == 0:
        BBox_in = array([[0, 0], [0, 0], [0, 0]], dtype="int32")
    if not with_2lpt:
        sx2_full = zeros((0, 0, 0), dtype="float32")
        sy2_full = zeros((0, 0, 0), dtype="float32")
        sz2_full = zeros((0, 0, 0), dtype="float32")

    CICDeposit_3(
        sx_full,
        sy_full,
        sz_full,
        sx2_full,
        sy2_full,
        sz2_full,
        density,
        cellsize,
        gridcellsize,
        1,
        growth_2pt_calc,
        L2,
        BBox_in,
        offset,
        1,
    )

    if cellsize_zoom != 0:
        CICDeposit_3(
            sx_full_zoom,
            sy_full_zoom,
            sz_full_zoom,
            sx2_full_zoom,
            sy2_full_zoom,
            sz2_full_zoom,
            density,
            cellsize_zoom,
            gridcellsize,
            1,
            growth_2pt_calc,
            L2,
            array([[0, 0], [0, 0], [0, 0]], dtype="int32"),
            offset_zoom,
            1,
        )

    density -= 1.0

    # print "den ic",density.mean(dtype=float64)
    if with_4pt_rule:
        density *= -0.5 / growth_2pt_calc / (1.0 - 1.0 / factor_4pt / factor_4pt)
    else:
        density *= -0.5 / growth_2pt_calc

    get_phi(density, den_k, den_fft, phi_fft, ngrid_x, ngrid_y, ngrid_z, gridcellsize)
    phi = density  # density now holds phi, so rename it
    grad_phi(
        sx_full,
        sy_full,
        sz_full,
        sx2_full,
        sy2_full,
        sz2_full,
        sx,
        sy,
        sz,
        npart_x,
        npart_y,
        npart_z,
        phi,
        ngrid_x,
        ngrid_y,
        ngrid_z,
        cellsize,
        gridcellsize,
        growth_2pt_calc,
        L2,
        offset,
    )

    if cellsize_zoom != 0:
        grad_phi(
            sx_full_zoom,
            sy_full_zoom,
            sz_full_zoom,
            sx2_full_zoom,
            sy2_full_zoom,
            sz2_full_zoom,
            sx_zoom,
            sy_zoom,
            sz_zoom,
            npart_x_zoom,
            npart_y_zoom,
            npart_z_zoom,
            phi,
            ngrid_x,
            ngrid_y,
            ngrid_z,
            cellsize_zoom,
            gridcellsize,
            growth_2pt_calc,
            L2,
            offset_zoom,
        )

    #######
    density.fill(0.0)

    CICDeposit_3(
        sx_full,
        sy_full,
        sz_full,
        sx2_full,
        sy2_full,
        sz2_full,
        density,
        cellsize,
        gridcellsize,
        1,
        -growth_2pt_calc,
        L2,
        BBox_in,
        offset,
        1,
    )

    if cellsize_zoom != 0:
        CICDeposit_3(
            sx_full_zoom,
            sy_full_zoom,
            sz_full_zoom,
            sx2_full_zoom,
            sy2_full_zoom,
            sz2_full_zoom,
            density,
            cellsize_zoom,
            gridcellsize,
            1,
            -growth_2pt_calc,
            L2,
            array([[0, 0], [0, 0], [0, 0]], dtype="int32"),
            offset_zoom,
            1,
        )

    density -= 1.0

    # print "den ic",density.mean(dtype=float64)
    if with_4pt_rule:
        density *= -0.5 / growth_2pt_calc / (1.0 - 1.0 / factor_4pt / factor_4pt)
    else:
        density *= -0.5 / growth_2pt_calc

    get_phi(density, den_k, den_fft, phi_fft, ngrid_x, ngrid_y, ngrid_z, gridcellsize)
    phi = density  # density now holds phi, so rename it

    grad_phi(
        sx_full,
        sy_full,
        sz_full,
        sx2_full,
        sy2_full,
        sz2_full,
        sx_minus,
        sy_minus,
        sz_minus,
        npart_x,
        npart_y,
        npart_z,
        phi,
        ngrid_x,
        ngrid_y,
        ngrid_z,
        cellsize,
        gridcellsize,
        -growth_2pt_calc,
        L2,
        offset,
    )
    if cellsize_zoom != 0:
        grad_phi(
            sx_full_zoom,
            sy_full_zoom,
            sz_full_zoom,
            sx2_full_zoom,
            sy2_full_zoom,
            sz2_full_zoom,
            sx_minus_zoom,
            sy_minus_zoom,
            sz_minus_zoom,
            npart_x_zoom,
            npart_y_zoom,
            npart_z_zoom,
            phi,
            ngrid_x,
            ngrid_y,
            ngrid_z,
            cellsize_zoom,
            gridcellsize,
            -growth_2pt_calc,
            L2,
            offset_zoom,
        )

    ######
    ######
    ###### Two more force evaluations in the case of a 4pt rule.
    ######
    ######

    if with_4pt_rule:
        density.fill(0.0)
        CICDeposit_3(
            sx_full,
            sy_full,
            sz_full,
            sx2_full,
            sy2_full,
            sz2_full,
            density,
            cellsize,
            gridcellsize,
            1,
            growth_2pt_calc * factor_4pt,
            L2 * factor_4pt * factor_4pt,
            BBox_in,
            offset,
            1,
        )
        if cellsize_zoom != 0:
            CICDeposit_3(
                sx_full_zoom,
                sy_full_zoom,
                sz_full_zoom,
                sx2_full_zoom,
                sy2_full_zoom,
                sz2_full_zoom,
                density,
                cellsize_zoom,
                gridcellsize,
                1,
                growth_2pt_calc * factor_4pt,
                L2 * factor_4pt * factor_4pt,
                array([[0, 0], [0, 0], [0, 0]], dtype="int32"),
                offset_zoom,
                1,
            )

        density -= 1.0

        # print "den ic",density.mean(dtype=float64)

        density *= (
            -0.5
            / growth_2pt_calc
            / (1.0 - 1.0 / factor_4pt / factor_4pt)
            * (-1.0 / factor_4pt ** 3)
        )

        get_phi(
            density, den_k, den_fft, phi_fft, ngrid_x, ngrid_y, ngrid_z, gridcellsize
        )
        phi = density  # density now holds phi, so rename it
        grad_phi(
            sx_full,
            sy_full,
            sz_full,
            sx2_full,
            sy2_full,
            sz2_full,
            sx_4pt,
            sy_4pt,
            sz_4pt,
            npart_x,
            npart_y,
            npart_z,
            phi,
            ngrid_x,
            ngrid_y,
            ngrid_z,
            cellsize,
            gridcellsize,
            growth_2pt_calc * factor_4pt,
            L2 * factor_4pt * factor_4pt,
            offset,
        )
        if cellsize_zoom != 0:
            grad_phi(
                sx_full_zoom,
                sy_full_zoom,
                sz_full_zoom,
                sx2_full_zoom,
                sy2_full_zoom,
                sz2_full_zoom,
                sx_4pt_zoom,
                sy_4pt_zoom,
                sz_4pt_zoom,
                npart_x_zoom,
                npart_y_zoom,
                npart_z_zoom,
                phi,
                ngrid_x,
                ngrid_y,
                ngrid_z,
                cellsize_zoom,
                gridcellsize,
                growth_2pt_calc * factor_4pt,
                L2 * factor_4pt * factor_4pt,
                offset_zoom,
            )

        #######
        density.fill(0.0)

        CICDeposit_3(
            sx_full,
            sy_full,
            sz_full,
            sx2_full,
            sy2_full,
            sz2_full,
            density,
            cellsize,
            gridcellsize,
            1,
            -growth_2pt_calc * factor_4pt,
            L2 * factor_4pt * factor_4pt,
            BBox_in,
            offset,
            1,
        )
        if cellsize_zoom != 0:
            CICDeposit_3(
                sx_full_zoom,
                sy_full_zoom,
                sz_full_zoom,
                sx2_full_zoom,
                sy2_full_zoom,
                sz2_full_zoom,
                density,
                cellsize_zoom,
                gridcellsize,
                1,
                -growth_2pt_calc * factor_4pt,
                L2 * factor_4pt * factor_4pt,
                array([[0, 0], [0, 0], [0, 0]], dtype="int32"),
                offset_zoom,
                1,
            )

        density -= 1.0

        density *= (
            -0.5
            / growth_2pt_calc
            / (1.0 - 1.0 / factor_4pt / factor_4pt)
            * (-1.0 / factor_4pt ** 3)
        )

        get_phi(
            density, den_k, den_fft, phi_fft, ngrid_x, ngrid_y, ngrid_z, gridcellsize
        )
        phi = density  # density now holds phi, so rename it

        grad_phi(
            sx_full,
            sy_full,
            sz_full,
            sx2_full,
            sy2_full,
            sz2_full,
            sx_4pt_minus,
            sy_4pt_minus,
            sz_4pt_minus,
            npart_x,
            npart_y,
            npart_z,
            phi,
            ngrid_x,
            ngrid_y,
            ngrid_z,
            cellsize,
            gridcellsize,
            -growth_2pt_calc * factor_4pt,
            L2 * factor_4pt * factor_4pt,
            offset,
        )
        if cellsize_zoom != 0:
            grad_phi(
                sx_full_zoom,
                sy_full_zoom,
                sz_full_zoom,
                sx2_full_zoom,
                sy2_full_zoom,
                sz2_full_zoom,
                sx_4pt_minus_zoom,
                sy_4pt_minus_zoom,
                sz_4pt_minus_zoom,
                npart_x_zoom,
                npart_y_zoom,
                npart_z_zoom,
                phi,
                ngrid_x,
                ngrid_y,
                ngrid_z,
                cellsize_zoom,
                gridcellsize,
                -growth_2pt_calc * factor_4pt,
                L2 * factor_4pt * factor_4pt,
                offset_zoom,
            )

    ###### Done with the two more force evaluations in the case of a 4pt rule.

    del density, den_k, den_fft, phi, phi_fft

    if (
        cellsize_zoom != 0
    ):  # the variables *_4pt* (except factor_4pt) are init'd to zero if 4pt rule is not requested
        sx2_zoom = (
            (sx_zoom + sx_minus_zoom + (sx_4pt_zoom + sx_4pt_minus_zoom) / factor_4pt)
            * cc
            / growth_2pt_calc
        )
        sy2_zoom = (
            (sy_zoom + sy_minus_zoom + (sy_4pt_zoom + sy_4pt_minus_zoom) / factor_4pt)
            * cc
            / growth_2pt_calc
        )
        sz2_zoom = (
            (sz_zoom + sz_minus_zoom + (sz_4pt_zoom + sz_4pt_minus_zoom) / factor_4pt)
            * cc
            / growth_2pt_calc
        )

        sx_zoom += sx_4pt_zoom - (sx_minus_zoom + sx_4pt_minus_zoom)
        sy_zoom += sy_4pt_zoom - (sy_minus_zoom + sy_4pt_minus_zoom)
        sz_zoom += sz_4pt_zoom - (sz_minus_zoom + sz_4pt_minus_zoom)

    sx2 = (sx + sx_minus + (sx_4pt + sx_4pt_minus) / factor_4pt) * cc / growth_2pt_calc
    sy2 = (sy + sy_minus + (sy_4pt + sy_4pt_minus) / factor_4pt) * cc / growth_2pt_calc
    sz2 = (sz + sz_minus + (sz_4pt + sz_4pt_minus) / factor_4pt) * cc / growth_2pt_calc

    sx += sx_4pt - (sx_minus + sx_4pt_minus)
    sy += sy_4pt - (sy_minus + sy_4pt_minus)
    sz += sz_4pt - (sz_minus + sz_4pt_minus)

    del sx_minus, sy_minus, sz_minus
    if cellsize_zoom != 0:
        del sx_minus_zoom, sy_minus_zoom, sz_minus_zoom
    if with_4pt_rule:
        del sx_4pt_minus, sy_4pt_minus, sz_4pt_minus
        if cellsize_zoom != 0:
            del sx_4pt_minus_zoom, sy_4pt_minus_zoom, sz_4pt_minus_zoom
        del sx_4pt, sy_4pt, sz_4pt
        if cellsize_zoom != 0:
            del sx_4pt_zoom, sy_4pt_zoom, sz_4pt_zoom

    # The lines below fix the box-size irrotational condition for periodic boxes

    # Mx=sx.mean(axis=0,dtype=float64)
    # My=sy.mean(axis=1,dtype=float64)
    # Mz=sz.mean(axis=2,dtype=float64)
    # for i in range(npart_x):
    # sx[i,:,:]-=Mx
    # for i in range(npart_y):
    # sy[:,i,:]-=My
    # for i in range(npart_z):
    # sz[:,:,i]-=Mz

    # del Mx,My,Mz
    # Mx=sx2.mean(axis=0,dtype=float64)
    # My=sy2.mean(axis=1,dtype=float64)
    # Mz=sz2.mean(axis=2,dtype=float64)
    # for i in range(npart_x):
    # sx2[i,:,:]-=Mx
    # for i in range(npart_y):
    # sy2[:,i,:]-=My
    # for i in range(npart_z):
    # sz2[:,:,i]-=Mz

    # del Mx,My,Mz
    # Mx=sx_zoom.mean(axis=0,dtype=float64)
    # My=sy_zoom.mean(axis=1,dtype=float64)
    # Mz=sz_zoom.mean(axis=2,dtype=float64)
    # for i in range(npart_x_zoom):
    # sx_zoom[i,:,:]-=Mx
    # for i in range(npart_y_zoom):
    # sy_zoom[:,i,:]-=My
    # for i in range(npart_z_zoom):
    # sz_zoom[:,:,i]-=Mz

    # del Mx,My,Mz
    # Mx=sx2_zoom.mean(axis=0,dtype=float64)
    # My=sy2_zoom.mean(axis=1,dtype=float64)
    # Mz=sz2_zoom.mean(axis=2,dtype=float64)
    # for i in range(npart_x_zoom):
    # sx2_zoom[i,:,:]-=Mx
    # for i in range(npart_y_zoom):
    # sy2_zoom[:,i,:]-=My
    # for i in range(npart_z_zoom):
    # sz2_zoom[:,:,i]-=Mz
    # del Mx,My,Mz

    if cellsize_zoom != 0:
        return (
            sx,
            sy,
            sz,
            sx2,
            sy2,
            sz2,
            sx_zoom,
            sy_zoom,
            sz_zoom,
            sx2_zoom,
            sy2_zoom,
            sz2_zoom,
        )
    else:
        return sx, sy, sz, sx2, sy2, sz2

ic_za(pspec0, boxsize=100.0, npart=64, init_seed=1234)

Generates Gaussian initial conditions for cosmological simulations in the Zel'dovich appoximation (ZA) -- the first order in Lagrangian Perturbation Theory (LPT).

Arguments:

  • pspec0 -- a string or callable function. If a string, this gives the filename for the plain text file containing the matter power spectrum at redshift zero from CAMB <http://www.camb.info/>_. For an example, see the included :download:camb_matterpower_z0.dat <./camb_matterpower_z0.dat>.

  • boxsize -- a float (default: 100.0). Gives the size of the simulation box in Mpc/h.

  • npart -- an integer (default: 64). The total number of particles is npart:sup:3.

  • init_seed -- an integer (default: 1234). The seed for the random number generator.

Return:

  • (sx,sy,sz) -- a tuple, where s:sub:i is a 3-dim single precision NumPy array containing the i-th component (i = x, y, z) of the particle displacements today as calculated in the ZA.

Example:

Generate a realization for the displacement field; calculate the rms displacements; and show a projection of one of the components.

>>> from ic import ic_za
>>> sx,sy,sz=ic_za('camb_matterpower_z0.dat')
Memory allocation done
Plans created
Power spectrum read.
Randoms done.
Nyquists fixed
sx fft ready
sy fft ready
sz fft ready
>>> ((sx**2+sy**2+sz**2).mean())**0.5/0.7 # O(10) for our universe
10.346006222040094
>>> import matplotlib.pyplot as plt # needs matplotlib to be installed!
>>> plt.imshow(sx.mean(axis=2))
<matplotlib.image.AxesImage object at 0x7fc4603697d0>
>>> plt.show()

Algorithm:

Implemented in the usual fft way.

.. warning:: This function has been tested but not at the level of trusting it for doing research. Use at your own risk.

Source code in pycola3/ic.py
def ic_za(pspec0, boxsize=100.0, npart=64, init_seed=1234):
    """
    Generates Gaussian initial conditions for cosmological
    simulations in the Zel'dovich appoximation (ZA) -- the first order
    in Lagrangian Perturbation Theory (LPT).

    **Arguments**:

    * ``pspec0`` --  a string or callable function. If a string, this gives
      the filename for the plain text file containing the matter power
      spectrum at redshift zero from `CAMB <http://www.camb.info/>`_.
      For an example, see the included
      :download:`camb_matterpower_z0.dat <./camb_matterpower_z0.dat>`.

    * ``boxsize`` -- a float (default: ``100.0``). Gives the size of the
      simulation box in Mpc/h.

    * ``npart`` --   an integer (default: ``64``). The total number of
      particles is ``npart``:sup:`3`.

    * ``init_seed`` --   an integer (default: ``1234``). The seed for the
      random number generator.

    **Return**:

    * ``(sx,sy,sz)`` --   a tuple, where ``s``:sub:`i` is a 3-dim single
      precision NumPy array containing the ``i``-th component
      (``i`` = ``x``, ``y``, ``z``) of the particle displacements today as calculated
      in the ZA.


    **Example**:

    Generate a realization for the displacement field; calculate the
    rms displacements; and show a projection of one of the components.

        >>> from ic import ic_za
        >>> sx,sy,sz=ic_za('camb_matterpower_z0.dat')
        Memory allocation done
        Plans created
        Power spectrum read.
        Randoms done.
        Nyquists fixed
        sx fft ready
        sy fft ready
        sz fft ready
        >>> ((sx**2+sy**2+sz**2).mean())**0.5/0.7 # O(10) for our universe
        10.346006222040094
        >>> import matplotlib.pyplot as plt # needs matplotlib to be installed!
        >>> plt.imshow(sx.mean(axis=2))
        <matplotlib.image.AxesImage object at 0x7fc4603697d0>
        >>> plt.show()



    **Algorithm**:

        Implemented in the usual fft way.

    .. warning:: This function has been tested but not at the level of
       trusting it for doing research. Use at your own risk.

    """
    delta = 2.0 * np.pi / boxsize
    nyq = npart // 2

    nalign = pyfftw.simd_alignment

    npart_pad = 2 * (npart // 2 + 1)

    sx_pad = pyfftw.n_byte_align_empty((npart, npart, npart_pad), nalign, "float32")
    sy_pad = pyfftw.n_byte_align_empty((npart, npart, npart_pad), nalign, "float32")
    sz_pad = pyfftw.n_byte_align_empty((npart, npart, npart_pad), nalign, "float32")

    sx = sx_pad[:, :, :npart]
    sy = sy_pad[:, :, :npart]
    sz = sz_pad[:, :, :npart]

    sx_k = sx_pad.view("complex64")
    sy_k = sy_pad.view("complex64")
    sz_k = sz_pad.view("complex64")

    print("Memory allocation done")

    nthreads = cpu_count()

    sx_fft = pyfftw.FFTW(
        sx_k,
        sx,
        axes=(0, 1, 2),
        direction="FFTW_BACKWARD",
        threads=nthreads,
        flags=("FFTW_ESTIMATE", "FFTW_DESTROY_INPUT"),
    )
    sy_fft = pyfftw.FFTW(
        sy_k,
        sy,
        axes=(0, 1, 2),
        direction="FFTW_BACKWARD",
        threads=nthreads,
        flags=("FFTW_ESTIMATE", "FFTW_DESTROY_INPUT"),
    )
    sz_fft = pyfftw.FFTW(
        sz_k,
        sz,
        axes=(0, 1, 2),
        direction="FFTW_BACKWARD",
        threads=nthreads,
        flags=("FFTW_ESTIMATE", "FFTW_DESTROY_INPUT"),
    )

    print("Plans created")

    # Determine whether power spectrum is function or filename
    if isinstance(pspec0, str):
        p_of_k = _power_spectrum(pspec0)
    else:
        assert callable(pspec0), "pspec0 must be a filename or a callable function"
        p_of_k = pspec0

    print("Power spectrum read.")

    np.random.seed(int(init_seed))

    x, y, z = np.indices((nyq + 1, nyq + 1, nyq + 1), dtype="float32")

    c = [(npart - i) % npart for i in range(nyq + 1)]

    w2 = x * x + y * y + z * z

    w2[0, 0, 0] = 1.0  # irrelevant but a zero crashes

    amp = np.sqrt(p_of_k(np.sqrt(w2) * delta) / boxsize ** 3)
    amp *= boxsize / (2.0 * np.pi)  # fix dimensions of 1/k

    phi = np.exp(
        1j * 2.0 * np.pi * np.random.uniform(0.0, 1.0, (nyq + 1, nyq + 1, nyq + 1))
    )
    phi *= np.random.normal(0.0, 1.0, (nyq + 1, nyq + 1, nyq + 1))

    sx_k[0 : nyq + 1, 0 : nyq + 1, 0 : nyq + 1] = x * phi / w2 * amp
    sy_k[0 : nyq + 1, 0 : nyq + 1, 0 : nyq + 1] = y * phi / w2 * amp
    sz_k[0 : nyq + 1, 0 : nyq + 1, 0 : nyq + 1] = z * phi / w2 * amp

    del phi

    phi = np.exp(
        1j * 2.0 * np.pi * np.random.uniform(0.0, 1.0, (nyq + 1, nyq + 1, nyq + 1))
    )
    phi *= np.random.normal(0.0, 1.0, (nyq + 1, nyq + 1, nyq + 1))

    sx_k[c, 0 : nyq + 1, 0 : nyq + 1] = -x * phi / w2 * amp
    sy_k[c, 0 : nyq + 1, 0 : nyq + 1] = y * phi / w2 * amp
    sz_k[c, 0 : nyq + 1, 0 : nyq + 1] = z * phi / w2 * amp

    del phi

    phi = np.exp(
        1j * 2.0 * np.pi * np.random.uniform(0.0, 1.0, (nyq + 1, nyq + 1, nyq + 1))
    )
    phi *= np.random.normal(0.0, 1.0, (nyq + 1, nyq + 1, nyq + 1))

    sx_k[0 : nyq + 1, c, 0 : nyq + 1] = x * phi / w2 * amp
    sy_k[0 : nyq + 1, c, 0 : nyq + 1] = -y * phi / w2 * amp
    sz_k[0 : nyq + 1, c, 0 : nyq + 1] = z * phi / w2 * amp

    del phi

    phi = np.exp(
        1j * 2.0 * np.pi * np.random.uniform(0.0, 1.0, (nyq + 1, nyq + 1, nyq + 1))
    )
    phi *= np.random.normal(0.0, 1.0, (nyq + 1, nyq + 1, nyq + 1))

    tmp = -x * phi / w2 * amp
    sx_k[npart - 1 : nyq - 1 : -1, 0, :] = tmp[1 : nyq + 1, 0, :]
    sx_k[0, npart - 1 : nyq - 1 : -1, :] = tmp[0, 1 : nyq + 1, :]
    sx_k[0, 0, :] = tmp[0, 0, :]

    del tmp

    tmp = -y * phi / w2 * amp
    sy_k[npart - 1 : nyq - 1 : -1, 0, :] = tmp[1 : nyq + 1, 0, :]
    sy_k[0, npart - 1 : nyq - 1 : -1, :] = tmp[0, 1 : nyq + 1, :]
    sy_k[0, 0, :] = tmp[0, 0, :]

    del tmp

    tmp = z * phi / w2 * amp
    sz_k[npart - 1 : nyq - 1 : -1, 0, :] = tmp[1 : nyq + 1, 0, :]
    sz_k[0, npart - 1 : nyq - 1 : -1, :] = tmp[0, 1 : nyq + 1, :]
    sz_k[0, 0, :] = tmp[0, 0, :]

    del phi, w2, amp, tmp

    print("Randoms done.")

    sx_k[0, 0, 0] = 0
    sy_k[0, 0, 0] = 0
    sz_k[0, 0, 0] = 0

    sx_k[npart - 1 : 0 : -1, npart - 1 : nyq - 1 : -1, [0, nyq]] = (
        sx_k[1:npart, 1 : nyq + 1, [0, nyq]]
    ).conjugate()
    sx_k[0, npart - 1 : nyq - 1 : -1, [0, nyq]] = (
        sx_k[0, 1 : nyq + 1, [0, nyq]]
    ).conjugate()
    sx_k[npart - 1 : 0 : -1, 0, [0, nyq]] = (sx_k[1:npart, 0, [0, nyq]]).conjugate()

    sy_k[npart - 1 : 0 : -1, npart - 1 : nyq - 1 : -1, [0, nyq]] = (
        sy_k[1:npart, 1 : nyq + 1, [0, nyq]]
    ).conjugate()
    sy_k[0, npart - 1 : nyq - 1 : -1, [0, nyq]] = (
        sy_k[0, 1 : nyq + 1, [0, nyq]]
    ).conjugate()
    sy_k[npart - 1 : 0 : -1, 0, [0, nyq]] = (sy_k[1:npart, 0, [0, nyq]]).conjugate()

    sz_k[npart - 1 : 0 : -1, npart - 1 : nyq - 1 : -1, [0, nyq]] = (
        sz_k[1:npart, 1 : nyq + 1, [0, nyq]]
    ).conjugate()
    sz_k[0, npart - 1 : nyq - 1 : -1, [0, nyq]] = (
        sz_k[0, 1 : nyq + 1, [0, nyq]]
    ).conjugate()
    sz_k[npart - 1 : 0 : -1, 0, [0, nyq]] = (sz_k[1:npart, 0, [0, nyq]]).conjugate()

    idx = np.where((x % nyq) + (y % nyq) + (z % nyq) == 0)
    del x, y, z

    sx_k[idx] = sx_k[idx].real
    sy_k[idx] = sy_k[idx].real
    sz_k[idx] = sz_k[idx].real

    del idx

    print("Nyquists fixed")

    sx_fft(normalise_idft=False)
    # sX=sx.copy()
    del sx_k, sx_fft, sx_pad
    print("sx fft ready")

    sy_fft(normalise_idft=False)
    # sY=sy.copy()
    del sy_k, sy_fft, sy_pad
    print("sy fft ready")

    sz_fft(normalise_idft=False)
    # sZ=sz.copy()
    del sz_k, sz_fft, sz_pad
    print("sz fft ready")

    return sx, sy, sz

import_music_snapshot(hdf5_filename, boxsize, level0='09', level1=None)

:math:\vspace{-1mm}

Import a MUSIC snapshot calculated in the ZA.

Arguments:

  • hdf5_filename -- a string. Gives the filename for the HDF5 <http://www.hdfgroup.org/HDF5/>_ file, which MUSIC outputs.

  • boxsize -- a float. The size of the full simulation box in :math:\mathrm{Mpc}/h.

  • level0 -- a two-character string (default: '09'). A MUSIC level covering the whole box. With the settings below, it should equal levelmin from the MUSIC configuration file for the finest such level.

  • level1 -- a two-character string (default: None). A fine MUSIC level covering the refined subvolume. With the settings below, it should equal levelmax from the MUSIC configuration file for the finest such level.

Return:

  • if level1 is None: (sx,sy,sz) -- a tuple, where s:sub:i is a 3-dim single precision NumPy array containing the i-th component (i = x, y, z) of the particle displacements today as calculated in the ZA. s:sub:i are the displacements for the level0 particles.

  • if level1 is not None: (sx,sy,sz,sx_zoom,sy_zoom,sz_zoom,offset) -- a tuple, where:

  • s:sub:i and s:sub:i\ _zoom are 3-dim single precision NumPy arrays containing the i-th component (i = x, y, z) of the particle displacements today as calculated in the ZA. s:sub:i are the displacements for the crude level (level0) particles; while s:sub:i\ _zoom are the displacements for the fine level (level1) particles in the refined subvolume.

  • offset -- a list of three integers giving the crude-grid index coordinates of the origin of the fine grid.

.. note:: pyCOLA requires specific values for some keywords in the MUSIC configuration file. Those are::

   zstart = 0
   align_top = yes
   use_2LPT = no
   format = generic

Also, if level1 is not None, pyCOLA assumes that one uses only one (usually the finest) refinement level (level1) on the subvolume of interest. Then the following needs to hold::

   levelmin<levelmax
   ref_extent!=1.0,1.0,1.0

See the included :download:ics.conf <./ics.conf> for an example.

Source code in pycola3/ic.py
def import_music_snapshot(hdf5_filename, boxsize, level0="09", level1=None):
    """
    :math:`\\vspace{-1mm}`

    Import a MUSIC snapshot calculated in the ZA.

    **Arguments**:

    * ``hdf5_filename`` --  a string. Gives the filename for the `HDF5 <http://www.hdfgroup.org/HDF5/>`_
      file, which MUSIC outputs.

    * ``boxsize`` -- a float. The size of the full simulation box in :math:`\mathrm{Mpc}/h`.

    * ``level0`` --   a two-character string (default: ``'09'``). A
      MUSIC level covering the whole box. With the settings below, it should
      equal ``levelmin`` from the MUSIC configuration file for the
      finest such level.

    * ``level1`` --   a two-character string (default: ``None``). A fine
      MUSIC level covering the refined subvolume. With the settings below, it
      should equal ``levelmax`` from the MUSIC configuration file for
      the finest such level.


    **Return**:

    * if ``level1`` is ``None``: ``(sx,sy,sz)`` --   a tuple, where ``s``:sub:`i`
      is a 3-dim single precision NumPy array containing the ``i``-th
      component (``i`` = ``x``, ``y``, ``z``) of the particle
      displacements today as calculated in the ZA. ``s``:sub:`i` are the
      displacements for the ``level0`` particles.

    * if ``level1`` is not ``None``:
      ``(sx,sy,sz,sx_zoom,sy_zoom,sz_zoom,offset)`` --   a tuple, where:

      - ``s``:sub:`i` and ``s``:sub:`i`\ ``_zoom`` are 3-dim single precision NumPy arrays containing
        the ``i``-th component (``i`` = ``x``, ``y``, ``z``) of the
        particle displacements today as calculated in the ZA. ``s``:sub:`i` are
        the displacements for the crude level (``level0``) particles;
        while ``s``:sub:`i`\ ``_zoom`` are the displacements for the fine level
        (``level1``) particles in the refined subvolume.

      - ``offset`` -- a list of three integers giving the crude-grid
        index coordinates of the origin of the fine grid.


    .. note::
       pyCOLA requires specific values for some keywords in the MUSIC
       configuration file. Those are::

           zstart = 0
           align_top = yes
           use_2LPT = no
           format = generic

       Also, if ``level1`` is not ``None``, pyCOLA assumes that one
       uses only one (usually the finest) refinement level (``level1``) on the subvolume
       of interest. Then the following needs to hold::

           levelmin<levelmax
           ref_extent!=1.0,1.0,1.0

       See the included :download:`ics.conf <./ics.conf>` for an example.



    """
    print("Starting import ...")
    ss = h5py.File(hdf5_filename, "r")

    # for some reason MUSIC pads with 4 elements the displacement arrays
    # when `format = generic` in the MUSIC conf file.
    # In my checks, this did not depend on the settings for the
    # padding and overlap keywords. So, hardwiring this number...

    sx = ss["level_0" + level0 + "_DM_dx"].value[4:-4, 4:-4, 4:-4] * boxsize
    sy = ss["level_0" + level0 + "_DM_dy"].value[4:-4, 4:-4, 4:-4] * boxsize
    sz = ss["level_0" + level0 + "_DM_dz"].value[4:-4, 4:-4, 4:-4] * boxsize
    if not (level1 is None):
        offset = [
            ss["header"]["grid_off_x"].value[-1],
            ss["header"]["grid_off_y"].value[-1],
            ss["header"]["grid_off_z"].value[-1],
        ]
        sx_zoom = ss["level_0" + level1 + "_DM_dx"].value[4:-4, 4:-4, 4:-4] * boxsize
        sy_zoom = ss["level_0" + level1 + "_DM_dy"].value[4:-4, 4:-4, 4:-4] * boxsize
        sz_zoom = ss["level_0" + level1 + "_DM_dz"].value[4:-4, 4:-4, 4:-4] * boxsize
        del ss
        print("... done")
        return sx, sy, sz, sx_zoom, sy_zoom, sz_zoom, offset
    else:
        del ss
        print("... done")
        return sx, sy, sz

initial_positions(sx, sy, sz, sx2, sy2, sz2, cellsize, growth_factor, growth_factor_2lpt, ngrid_x, ngrid_y, ngrid_z, gridcellsize, offset=[0.0, 0.0, 0.0])

:math:\vspace{-1mm}

Add the Lagrangian particle position to the 2LPT displacement to obtain the Eulerian position. Periodic boundary conditions are assumed.

Arguments:

  • sx,sy,sz -- 3-dim NumPy arrays containing the components of the particle displacements today as calculated in the ZA.

  • sx2,sy2,sz2 -- 3-dim NumPy arrays containing the components of the second order particle displacements today as calculated in 2LPT.

  • cellsize -- a float. The inter particle spacing in Lagrangian space.

  • growth_factor -- a float. The linear growth factor for the redshift for which the Eulerian positions are requested.

  • growth_factor_2lpt -- a float. The second order growth factor for the redshift for which the Eulerian positions are requested.

  • ngrid_x, ngrid_y, ngrid_z -- integers. The grid size of the box. Only used together with gridcellsize below to find the physical size of the box, which is needed to apply the periodic boundary conditions.

  • gridcellsize -- a float. The grid spacing of the box.

  • offset -- a list of three floats (default: [0.0,0.0,0.0]). Offset the Eulerian particle positions by this amount. Useful for placing refined subregions at their proper locations inside a bigger box.

Return:

  • (px,py,pz) -- a tuple, where p:sub:i is a 3-dim single precision NumPy array containing the i-th component (i = x, y, z) of the particle Eulerian position.

Example:

In this example we generate the initial conditions in 2LPT, and then plot a slice through the 2LPT realization at redshift of zero.

>>> from pycola.ic import ic_za,ic_2lpt,initial_positions
>>> sx,sy,sz=ic_za('camb_matterpower_z0.dat',npart=128)
Memory allocation done
Plans created
Power spectrum read.
Randoms done.
Nyquists fixed
sx fft ready
sy fft ready
sz fft ready
>>> sx2,sy2,sz2 = ic_2lpt( 100.0/float(sx.shape[0]),sx,sy,sz,
...                        growth_2pt_calc=0.1)
>>> px,py,pz = initial_positions(sx,sy,sz,sx2,sy2,sz2,100./128.,1.0,1.0,
...            1, # only ngrid_i*gridcellsize=boxsize is relevant here
...            1,
...            1,
...            100.0)
>>> import matplotlib.pyplot as plt # needs matplotlib to be installed
>>> import numpy as np
>>> ind=np.where(pz<3)
>>> px_slice=px[ind]
>>> py_slice=py[ind]
>>> plt.figure(figsize=(10,10))
<matplotlib.figure.Figure object at 0x7f21044d2e10>
>>> plt.scatter(px_slice,py_slice,marker='.',alpha=0.03,color='r')
<matplotlib.collections.PathCollection object at 0x7f2102cd3290>
>>> plt.show()
Source code in pycola3/ic.py
def initial_positions(
    sx,
    sy,
    sz,
    sx2,
    sy2,
    sz2,
    cellsize,
    growth_factor,
    growth_factor_2lpt,
    ngrid_x,
    ngrid_y,
    ngrid_z,
    gridcellsize,
    offset=[0.0, 0.0, 0.0],
):
    """
    :math:`\\vspace{-1mm}`

    Add the Lagrangian particle position to the 2LPT displacement to
    obtain the Eulerian position. Periodic boundary conditions are assumed.


    **Arguments**:

    * ``sx,sy,sz`` --  3-dim NumPy arrays containing the
      components of the particle
      displacements today as calculated in the ZA.

    * ``sx2,sy2,sz2`` --  3-dim NumPy arrays containing the
      components of the second order particle
      displacements today as calculated in 2LPT.


    * ``cellsize`` --  a float. The inter particle spacing in Lagrangian space.

    * ``growth_factor`` --  a float. The linear growth factor for the
      redshift for which the Eulerian positions are requested.

    * ``growth_factor_2lpt`` --  a float. The second order growth factor for the
      redshift for which the Eulerian positions are requested.

    * ``ngrid_x``, ``ngrid_y``, ``ngrid_z`` -- integers. The grid size
      of the box. Only used together with ``gridcellsize`` below to find the
      physical size of the box, which is needed to apply the periodic
      boundary conditions.

    * ``gridcellsize`` -- a float. The grid spacing of the box.

    * ``offset`` -- a list of three floats (default: ``[0.0,0.0,0.0]``).
      Offset the Eulerian particle positions by this amount. Useful for
      placing refined subregions at their proper locations inside a
      bigger box.

    **Return**:

    * ``(px,py,pz)`` --   a tuple, where ``p``:sub:`i`
      is a 3-dim single precision NumPy array containing the ``i``-th
      component (``i`` = ``x``, ``y``, ``z``) of the particle
      Eulerian position.

    **Example**:

    In this example we generate the initial conditions in 2LPT, and
    then plot a slice through the 2LPT realization at redshift of zero.

        >>> from pycola.ic import ic_za,ic_2lpt,initial_positions
        >>> sx,sy,sz=ic_za('camb_matterpower_z0.dat',npart=128)
        Memory allocation done
        Plans created
        Power spectrum read.
        Randoms done.
        Nyquists fixed
        sx fft ready
        sy fft ready
        sz fft ready
        >>> sx2,sy2,sz2 = ic_2lpt( 100.0/float(sx.shape[0]),sx,sy,sz,
        ...                        growth_2pt_calc=0.1)
        >>> px,py,pz = initial_positions(sx,sy,sz,sx2,sy2,sz2,100./128.,1.0,1.0,
        ...            1, # only ngrid_i*gridcellsize=boxsize is relevant here
        ...            1,
        ...            1,
        ...            100.0)
        >>> import matplotlib.pyplot as plt # needs matplotlib to be installed
        >>> import numpy as np
        >>> ind=np.where(pz<3)
        >>> px_slice=px[ind]
        >>> py_slice=py[ind]
        >>> plt.figure(figsize=(10,10))
        <matplotlib.figure.Figure object at 0x7f21044d2e10>
        >>> plt.scatter(px_slice,py_slice,marker='.',alpha=0.03,color='r')
        <matplotlib.collections.PathCollection object at 0x7f2102cd3290>
        >>> plt.show()

    """
    npart_x, npart_y, npart_z = sx.shape

    px, py, pz = np.indices((npart_x, npart_y, npart_z), dtype="float32")

    px *= cellsize
    py *= cellsize
    pz *= cellsize

    px += float(ngrid_x) * gridcellsize + offset[0]
    py += float(ngrid_y) * gridcellsize + offset[1]
    pz += float(ngrid_z) * gridcellsize + offset[2]

    px += sx * growth_factor
    py += sy * growth_factor
    pz += sz * growth_factor

    px += sx2 * growth_factor_2lpt
    py += sy2 * growth_factor_2lpt
    pz += sz2 * growth_factor_2lpt

    px %= float(ngrid_x) * gridcellsize
    py %= float(ngrid_y) * gridcellsize
    pz %= float(ngrid_z) * gridcellsize

    return px, py, pz

Initial conditions

First order initial conditions for pyCOLA can be calculated using either MUSIC [MUSIC]_, or internally. With MUSIC, however, one can do refinements on a region, which is not supported internally.

The second-order displacement field is generated using a novel algorithm using force evaluations. See the Algorithm section of :func:ic_2lpt_engine for details.

.. warning:: As of MUSIC rev. 116353436ee6 <https://bitbucket.org/ohahn/music/commits/116353436ee6ff009abda2e51cd4792a209fa63b>, the second-order displacement field returned by MUSIC gets unphysical large-scale deviations when a refined subvolume is requested (seems to be fine for single grid). Until that problem is fixed, use the function :func:ic.ic_2lpt to get the second order displacements from the first order result. Update: MUSIC received a fix with rev. ed51fcaffee <https://bitbucket.org/ohahn/music/commits/ed51fcaffeec08d676b8a3435a0c4009f4024220>, which supposedly fixes the problem.

Initial displacements at first order ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

.. autofunction:: ic.import_music_snapshot

.. autofunction:: ic.ic_za

Initial displacements at second order ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

.. autofunction:: ic.ic_2lpt

.. autofunction:: ic.ic_2lpt_engine

Obtaining the Eulerian positions ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

.. autofunction:: ic.initial_positions

.. automodule:: ic :undoc-members: :show-inheritance: