def distances_squared(particles, other_particles):
    """
    Returns the distance squared from each particle in this set to each of the particles in the other set.

    :argument other_particles: the particles to which the distance squared is calculated
    
    >>> from amuse.datamodel import Particles
    >>> field_particles = Particles(2)
    >>> field_particles.x = [0.0, 2.0] | units.m
    >>> field_particles.y = [0.0, 0.0] | units.m
    >>> field_particles.z = [0.0, 0.0] | units.m
    >>> particles = Particles(3)
    >>> particles.x = [1.0, 3.0, 4.0] | units.m
    >>> particles.y = [0.0, 0.0, 0.0] | units.m
    >>> particles.z = [0.0, 0.0, 0.0] | units.m
    >>> particles.distances_squared(field_particles)
    quantity<[[1.0, 1.0], [9.0, 1.0], [16.0, 4.0]] m**2>
    """
    transposed_positions = particles.position.reshape((len(particles), 1, -1))
    dxdydz = transposed_positions - other_particles.position
    return (dxdydz**2).sum(-1)