1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
#!/usr/bin/env python
# -*- coding: utf8 -*-

import pylab
import threading

class MyThread(threading.Thread):
    def __init__(self,y_lower,y_upper):
        threading.Thread.__init__(self)
        self.y_lower= y_lower
        self.y_upper = y_upper

    def run(self):
        for xi in range(len(x)):
            for yi in range(self.y_lower,self.y_upper):
                P[yi,xi] = pylab.sum(
                    pylab.select(
                        [(x[xi]-X)*(x[yi]-Y)==0,],
                        [0,],
                        Q/pylab.sqrt((x[xi]-X)**2+(x[yi]-Y)**2)
                    )
                )

dx = .05
x = pylab.arange(-5,5,dx)
X,Y = pylab.meshgrid(x,x)

# Homogeneous charge density
Q = pylab.zeros_like(X)
Q = pylab.select([abs(X+2)<4*dx,],[pylab.select([abs(Y)<4*dx,],[-1,],0),],Q)
#Q += pylab.select([abs(X+.5)<.3,],[pylab.select([abs(Y)<.33],[-1,],0),],Q)

# Sharp rectangular charge density
#Q = pylab.zeros_like(X)
#Q[40,31:70] = 1
#Q[41:45,30] = 1
#Q[41:45,70] = 1
#Q[45,31:70] = 1

#Q[60,46:54] = -1
#Q[61:62,45] = -1
#Q[61:62,54] = -1
#Q[62,46:54] = -1

# Sharp circular charge density
#Q = pylab.zeros_like(X)
Q = pylab.select([abs((X-2)**2+(Y+2)**2-(1.5)**2)<pylab.sqrt(2*dx)],[1],Q)

P = pylab.zeros_like(X)

threads = []
dy = 4
for i in range(dy):
    y_lower = int(len(x)/dy*float(i))
    y_upper = int(len(x)/dy*float(i+1))
    t = MyThread(y_lower,y_upper)
    threads.append(t)
    t.start()

if y_upper < len(x):
    t = MyThread(y_upper,len(x))
    threads.append(t)
    t.start()

for t in threads:
    t.join()

#E = pylab.gradient(P)
pylab.gray()
pylab.pcolormesh(X,Y,Q)
pylab.spectral()
pylab.contourf(X,Y,P,levels=pylab.linspace(P.min(),P.max(),125),alpha=.8))
pylab.colorbar()
#pylab.quiver(X,Y,E[1],E[0],pivot='middle',alpha=.5,color='gray')
pylab.show()