/
stationary_charge.py
58 lines (48 loc) · 2.02 KB
/
stationary_charge.py
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
#!/usr/bin/env python
"""Module calculates and plots the E field components from a stationary charge.
The x, y, and z components generated by a stationary point charge at the origin
(0, 0, 0) are plotted along the x-y plane between -10 nm to 10 nm. Each field
component is plotted individually on a logarithmic scale.
"""
# pragma pylint: disable=unexpected-keyword-arg
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
import pycharge as pc
# Create charge and simulation objects
source = pc.StationaryCharge(position=(0, 0, 0))
simulation = pc.Simulation(source)
# Create meshgrid in x-y plane between -10 nm to 10 nm at z=0
lim = 10e-9
npoints = 1000 # Number of grid points
coordinates = np.linspace(-lim, lim, npoints) # grid from -lim to lim
x, y, z = np.meshgrid(coordinates, coordinates, 0, indexing='ij') # z=0
# Calculate E field components at t=0
E_x, E_y, E_z = simulation.calculate_E(t=0, x=x, y=y, z=z)
# Plot E_x, E_y, and E_z fields
E_x_plane = E_x[:, :, 0] # Create 2D array at z=0 for plotting
E_y_plane = E_y[:, :, 0]
E_z_plane = E_z[:, :, 0]
# Create figs and axes, plot E components on log scale
fig, axs = plt.subplots(1, 3, sharey=True)
norm = mpl.colors.SymLogNorm(linthresh=1.01e6, linscale=1, vmin=-1e9, vmax=1e9)
extent = [-lim, lim, -lim, lim]
im_0 = axs[0].imshow(E_x_plane.T, origin='lower', norm=norm, extent=extent)
im_1 = axs[1].imshow(E_y_plane.T, origin='lower', norm=norm, extent=extent)
im_2 = axs[2].imshow(E_z_plane.T, origin='lower', norm=norm, extent=extent)
# Add labels
for ax in axs:
ax.set_xlabel('x (nm)')
axs[0].set_ylabel('y (nm)')
axs[0].set_title('E_x')
axs[1].set_title('E_y')
axs[2].set_title('E_z')
# Add colorbar to figure
Ecax = inset_axes(
axs[2], width="6%", height="100%", loc='lower left',
bbox_to_anchor=(1.05, 0., 1, 1), bbox_transform=axs[2].transAxes, borderpad=0
)
E_cbar = plt.colorbar(im_2, cax=Ecax) # right of im_2
E_cbar.ax.set_ylabel('E (N/C)', rotation=270, labelpad=12)
plt.show()