import tensorflow as tf

    A = tf.reduce_sum(phi_sampling(A0/(4*y), D))

    B0 = euclidean_norm_squared(tf.subtract(tf.expand_dims(Y, 0), tf.expand_dims(Y, 1)), axis=2)
    B = tf.reduce_sum(phi_sampling(B0/(4*y), D))

    C0 = euclidean_norm_squared(tf.subtract(tf.expand_dims(X, 0), tf.expand_dims(Y, 1)), axis=2)
    C = tf.reduce_sum(phi_sampling(C0/(4*y), D))

    return T*(A + B - 2*C)
