y = np.arange(10)
v = np.arange(100).reshape(10, 10)
i, j = 5, 5
r = np.sqrt((x[:, (None)] - x[i]) ** 2 + (y[(None), :] - y[j]) ** 2)
np.sum(v[r <= 10])
Matrix calculations on equidistant elements without for loops