import math
import numpy as np

K = 9.0e9

# Calculate the potential at point (x, y).

def potential_py(x, y, q, xq, yq):	# Python version

    phi = 0.
    for i in range(len(q)):

        dx = x - xq[i]
        dy = y - yq[i]
        dr = math.sqrt(dx*dx + dy*dy)

        #phi += K*q[i]/dr
        phi -= K*q[i]*math.log(dr)

    return phi

def potential_np(x, y, q, xq, yq):	# numpy version

    dr = np.sqrt((x-xq)**2+(y-yq)**2)
    #phi = K*q[dr>0]/dr[dr>0]
    phi = -K*q[dr>0]*np.log(dr)

    return phi.sum()

if __name__ == "__main__" :

    N = 1000
    q = np.random.uniform(-1.e-6, 1.e-6, N)
    xq = np.random.uniform(-1., 1., N)
    yq = np.random.uniform(-1., 1., N)

    err = 0.
    ne = 100
    for x, y in zip(np.random.uniform(-1., 1., ne),
                    np.random.uniform(-1., 1., ne)):
        pot1 = potential_py(x, y, q, xq, yq)
        pot2 = potential_np(x, y, q, xq, yq)
        dp = (pot2-pot1)/pot1
        print('pot =', pot2, 'dpot =', dp)
        err += dp**2

    print('err =', math.sqrt(err/ne))