```
In [25]:
```import matplotlib.pyplot as plt
import numpy as np
import scipy.signal as signal

```
In [2]:
```#define grid
xmax=1 #extents of grid [m]
xpoints=256 #number of grid points to calculate
x=np.linspace(-xmax,xmax,xpoints)
y=np.linspace(-xmax,xmax,xpoints)
x2d,y2d=np.meshgrid(x,y)
r=np.sqrt(np.square(x2d)+np.square(y2d))

```
In [9]:
```gap=xmax/8 #defines gap
w=xmax/10
electrode1=np.multiply(r<gap/2,x2d<0) #half disc inside shielding ring
electrode2=np.multiply(r>gap,r<(gap+w)) #shielding ring
electrode3=np.sqrt(np.square(x2d-4*gap)+np.square(y2d))<gap #disc outside shielding
electrode4=np.sqrt(np.square(y2d-4*gap)+np.square(x2d))<gap #another disc outside shielding
#define the potential on each electrode
v1=10
v2=0 #assume shield is at gnd potential (arbitrary, but this is convention)
v3=1
v4=10
allelectrode=electrode1+electrode2+electrode3+electrode4
vac=allelectrode<1 #logic function defining the vacuum region where no electrodes are present
v=v1*electrode1+v2*electrode2+v3*electrode3+v4*electrode4 #potential boundary conditions=known constant potential on the electrodes
idelectrode=np.where(vac==0) #index for all points that contain electrodes
idvac=np.where(vac==1) #index for all points in vacuum

Show the electrode setup and their potential

```
In [37]:
```plt.imshow(allelectrode,extent=[-xmax,xmax,-xmax,xmax])
plt.axis('square')
plt.title('position of electrodes')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
plt.imshow(v,extent=[-xmax,xmax,-xmax,xmax])
plt.axis('square')
plt.title('potential on electrodes')
plt.xlabel('x')
plt.ylabel('y')
plt.colorbar()
plt.show()

```
```

```
In [38]:
```v[idvac]=np.mean(v[idelectrode]) #replace vacuum potential by the mean potential on the electrodes

```
In [39]:
```vn=v #starting condition is first rough estimate
kmax=1000 #number of itteration steps
kernel=np.array([[1/4,0,1/4],[0,0,0],[1/4,0,1/4]]) #approximative Laplace kernel replacing the central point by the average of 4 neighbouring points.
error=np.zeros(kmax)
for k in range(kmax):
vold=vn
vn=signal.convolve(vold,kernel,mode="same") #apply kernel to old estimate
vn[idelectrode]=v[idelectrode] #and reinforce the boundary conditions (on the electrodes we KNOW the potentials)
dif=vold-vn #difference between old and new estimate
error[k]=np.sum(np.square(dif)) #square difference as indicator of convergence
#imagesc(Vn); %see how the potential evolves while calculating (but
#takes a lot of time that should really go to the calculation and not
#to plotting
#pause(0.1); %to force refreshing everytime we come here, otherwise the image is not updated

```
In [40]:
```plt.semilogy(error)
plt.title('mean squared correction')
plt.xlabel('iteration')

```
Out[40]:
```

```
In [41]:
```plt.imshow(vn,extent=[-xmax,xmax,-xmax,xmax])
plt.colorbar()
plt.title('estimated V')
plt.axis('square')

```
Out[41]:
```

Now get the fields from E=-grad(V)

```
In [42]:
```ex,ey=np.gradient(-vn,x,y)
e=np.sqrt(np.square(ex)+np.square(ey))

```
In [43]:
```plt.imshow(e,extent=[-xmax,xmax,-xmax,xmax])
plt.colorbar()
plt.title('|E|')
plt.axis('square')

```
Out[43]:
```

```
In [45]:
```eps0=8.854e-12
dexdx,dexdy=np.gradient(ex,x,y)
deydx,deydy=np.gradient(ey,x,y)
div=dexdx+deydy #apparantly there is no divergence operator in numpy...so we write one ourselves
rho=eps0*div
plt.imshow(rho,extent=[-xmax,xmax,-xmax,xmax])
plt.colorbar()
plt.title('charge density')
plt.axis('square')

```
Out[45]:
```