Điện thế hệ các điện tích điểm

Chương trình này tính toán sự phân bố điện thế bao quanh hệ các điện tích điểm, kết hợp mô tả các đường sức điện trường. Bức tranh điện trường trở nên sinh động hơn nhờ những lát cắt xen kẽ. Ta có thể thấy quy luật đan xen giữa đường sức điện trường và các mặt đẳng thế: chúng luôn vuông góc với nhau.

Kết quả tính toán thể hiện như hình dưới.

Code chương trình

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
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101

function Equipotential_lines_Particles_Field
% Author: Tran Hai Cat
% Created: 2019.02.02
clc;
clear all;
close all;
%% INPUT DATA
Q = [2 -2]*1e-9; % charges

xQ = [-1 1]; % coordinates of charges
yQ = [0 0];
zQ = [0 0];

xmin = -2; xmax = 2; ymin = -1; ymax = 1; zmin = -1; zmax = 1;
%% CALCULATION
N = length(Q);
Nx = 50;
Ny = 20;
Nz = 20;
x = linspace(xmin,xmax,Nx);
y = linspace(ymin,ymax,Ny);
z = linspace(zmin,zmax,Nz);
[X,Y,Z] = meshgrid(x,y,z);

ke = 9e9;
Ex = zeros(size(X));
Ey = zeros(size(Y));
Ez = zeros(size(Z));
V = zeros(size(X));
for iN = 1:N
r = sqrt((X-xQ(iN)).^2+(Y-yQ(iN)).^2+(Z-zQ(iN)).^2);

E = ke*Q(iN)./r./r;
Ex = Ex+E./r.*(X-xQ(iN)); % x-component of vector field
Ey = Ey+E./r.*(Y-yQ(iN)); % y-component of vector field
Ez = Ez+E./r.*(Z-zQ(iN)); % z-component of vector field

condition = find(r<0.1);
Ex(condition) = 0;
Ey(condition) = 0;
Ez(condition) = 0;

V = V+ke*Q(iN)./r;

condition = r<0.3;
V(condition) = 0;
end

%% FIGURES
figure('name','Equipotential Lines Particles Field','color','w','numbertitle','off');
hold on
xslice = 0;
yslice = 0;
zslice = 0;
slice(X,Y,Z,V,xslice,yslice,zslice);
xslice = -1;
slice(X,Y,Z,V,xslice,yslice,zslice);
xslice = 1;
slice(X,Y,Z,V,xslice,yslice,zslice);
colormap(jet(64));
colorbar('vertical');
shading interp;

% Draw streamline:
for iN = 1:N
[x1, y1, z1] = sphere(12);
radius = 0.1;
x1 = x1(:)*radius+xQ(iN);
y1 = y1(:)*radius;
z1 = z1(:)*radius;

hlines = streamline(X,Y,Z,Ex,Ey,Ez,x1,y1,z1);
set(hlines,'Color','m')
end

% Draw charged particles:
for iN = 1:N
[x1, y1, z1] = sphere(24);
radius = 0.1;
x1 = x1(:)*radius+xQ(iN);
y1 = y1(:)*radius;
z1 = z1(:)*radius;
P = [x1 y1 z1];
P = unique(P,'rows');

shp = alphaShape(P,1.5);
plot(shp)
end

%% Axis properties
axis equal
axis([xmin xmax ymin ymax zmin zmax]);
rotate3d on

xlabel('X, m');
ylabel('Y, m');
zlabel('Z, m');
view(-30,20);
box off
grid off
set(gca,'color','w')