def laplacian(field, m, s, irho):
    """
    Spacial discretization for the isotropic acoustic wave equation. For a 4th
    order in time formulation, the 4th order time derivative is replaced by a
    double laplacian:
    H = (laplacian + s**2/12 laplacian(1/m*laplacian))

    Parameters
    ----------
    field : TimeFunction
        The computed solution.
    m : Function or float
        Square slowness.
    s : float or Scalar
        The time dimension spacing.
    """
    so = irho.space_order // 2
    Lap = sum([getattr(irho * getattr(field, 'd%s'%d.name)(x0=d + d.spacing/2, fd_order=so), 'd%s'% d.name)(x0=d - d.spacing/2, fd_order=so) for d in irho.dimensions]) # noqa

    return Lap