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
#!/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+Q_shadow)/pylab.sqrt((x[xi]-X)**2+(x[yi]-Y)**2)
                    )
                )

dx = .05
x = pylab.arange(-3.4,3.4,dx)
X,Y = pylab.meshgrid(x,x)
Q_shadow = pylab.zeros_like(X)
Q = pylab.zeros_like(X)

# Sharp circular charge density
P_ini = pylab.zeros_like(X)
x0, y0, r =2, 0, 1
Q = pylab.select([abs((X-x0)**2+(Y-y0)**2-r**2)<(5*dx)**2],[5],[0])
Q_shadow = -1*pylab.fliplr(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()

P = pylab.select([X>0],[P],[0])

pylab.spectral()
pylab.contourf(X,Y,P,levels=pylab.linspace(P.min()-1e-2,P.max(),96),alpha=1)
pylab.colorbar()
pylab.show()