def pillar_distances_sqr(grid, xy, ref_k0 = 0, kp = 0, horizon_points = None):
    """Returns array of the square of the distances of primary pillars in x,y plane to point xy.

    arguments:
       xy (float pair): the xy coordinate to compute the pillar distances to
       ref_k0 (int, default 0): the horizon layer number to use
       horizon_points (numpy array, optional): if present, should be array as returned by
          horizon_points() method; pass for efficiency in case of multiple calls
    """

    # note: currently works with unmasked data and using primary pillars only
    pe_j = grid.extent_kji[1] + 1
    pe_i = grid.extent_kji[2] + 1
    if horizon_points is None:
        horizon_points = grid.horizon_points(ref_k0 = ref_k0, kp = kp)
    pillar_xy = horizon_points[:, :, 0:2]
    dxy = pillar_xy - xy
    dxy2 = dxy * dxy
    return (dxy2[:, :, 0] + dxy2[:, :, 1]).reshape((pe_j, pe_i))