Phân bố điện tích trên hai quả cầu kim loại đặt gần nhau

function E_Field_Spheres_Q1_Q2
% Author: Tran Hai Cat
% Lecturer in Physics, HCM University of Technology and Education
% - Dai hoc Su pham Ky thuat Tp. Ho Chi Minh
% Created: 2019.08.05
clc;
clear variables;
close all;

%% CONST
ke = 1;

%% PARAMETERS:
R1 = 55e-2; % sphere radius
R2 = 55e-2;
x01 = -0.60;
y01 = 0;
z01 = 0;
x02 = 0.60;
y02 = 0;
z02 = 0;

q01 = 1; % charge per particle (left)
q02 = -1; % charge per particle (right)
N1 = 500; % Number of charges on the left sphere
N2 = 500; % Number of charges on the right sphere

ki = 5e5;

Nx = 100;
Ny = 50;

xmin = -2; xmax = 2; ymin = -1; ymax = 1; zmin = -1; zmax = 1;
%% CALCULATION
N = N1+N2;
q = zeros(1,N);
for i = 1:N1
q(i) = q01;
end
for i = N1+1:N
q(i) = q02;
end

x_e = zeros(1,N);
y_e = zeros(1,N);
z_e = zeros(1,N);
i = 0;
while i<N1
x = x01-R1+rand()*2*R1;
y = y01-R1+rand()*2*R1;
z = z01-R1+rand()*2*R1;
r = sqrt((x-x01).^2+(y-y01).^2+(z-z01).^2);
if r<R1
i = i+1;
x_e(i) = x;
y_e(i) = y;
z_e(i) = z;
end
end

while i<N
x = x02-R2+rand()*2*R2;
y = y02-R2+rand()*2*R2;
z = z02-R1+rand()*2*R2;
r = sqrt((x-x02).^2+(y-y02).^2+(z-z01).^2);
if r<R2
i = i+1;
x_e(i) = x;
y_e(i) = y;
z_e(i) = z;
end
end

figure('name','Electric Field','color','black','numbertitle','off');
hold on
% Draw spheres:
[x1, y1, z1] = sphere(24);
radius = R1;
x1 = x1(:)*radius+x01;
y1 = y1(:)*radius+y01;
z1 = z1(:)*radius+z01;
P = [x1 y1 z1];
P = unique(P,'rows');

shp = alphaShape(P,1.5);
plot(shp,'EdgeColor',[0.4 0.4 0.4])

[x1, y1, z1] = sphere(24);
radius = R1;
x1 = x1(:)*radius+x02;
y1 = y1(:)*radius+y02;
z1 = z1(:)*radius+z02;
P = [x1 y1 z1];
P = unique(P,'rows');

shp = alphaShape(P,1.5);
plot(shp,'EdgeColor',[0.4 0.4 0.4])
alpha(0.0);

e_plot = zeros(1,N);
for i = 1:N1
e_plot(i) = plot3(x_e(i),y_e(i),z_e(i),'ro','MarkerSize',2,'Markerfacecolor','r');
end
for i = N1+1:N
e_plot(i) = plot3(x_e(i),y_e(i),z_e(i),'go','MarkerSize',2,'Markerfacecolor','g');
end
xlabel('x [m]','fontsize',14);
ylabel('y [m]','fontsize',14);
% ht = title('Electric field');
set(gca,'fontsize',14,'color','k','xcolor','w','ycolor','w','zcolor','w');
axis([xmin xmax ymin ymax zmin zmax]);
axis equal
view(-30,30);
rotate3d on

dx_max = 1e-2;
dy_max = dx_max;
dz_max = dx_max;
for k = 1:ki
ii = randi([1 N]);
U_old = 0;
for i = 1:N
if i~=ii
dis = sqrt((x_e(ii)-x_e(i)).^2+(y_e(ii)-y_e(i)).^2+(z_e(ii)-z_e(i)).^2);
U_old = U_old+q(ii)*q(i)/dis;
end
end

r = 100;
while ((ii<=N1)&&(r>R1)) || ((ii>N1)&&(r>R2))
x_try = x_e(ii)-dx_max+rand()*2*dx_max;
y_try = y_e(ii)-dy_max+rand()*2*dy_max;
z_try = z_e(ii)-dz_max+rand()*2*dz_max;
if ii<=N1
r = sqrt((x_try-x01).^2+(y_try-y01).^2+(z_try-z01).^2);
else
r = sqrt((x_try-x02).^2+(y_try-y02).^2+(z_try-z02).^2);
end
end

U_new = 0;
for i = 1:N
if i~=ii
dis = sqrt((x_try-x_e(i)).^2+(y_try-y_e(i)).^2++(z_try-z_e(i)).^2);
U_new = U_new+q(ii)*q(i)/dis;
end
end

if U_new<U_old
x_e(ii) = x_try;
y_e(ii) = y_try;
z_e(ii) = z_try;

set(e_plot(ii),'xdata',x_try,'ydata',y_try,'zdata',z_try);
% set(ht,'string',sprintf('%0.0f',k));
fprintf('k=%d\n',k);
pause(0.002);
clc
end

end

x = linspace(xmin,xmax,Nx);
y = linspace(ymin,ymax,Ny);
V = zeros(Nx,Ny); % Potential (Voltage) matrix
for i=1:Nx
for j=1:Ny
V(i,j) = 0;
for k = 1:N
dis = sqrt((x(i)-x_e(k)).^2+(y(j)-y_e(k)).^2+z_e(k).^2);
V(i,j) = V(i,j)+ke*q(k)/dis;
end
end
end
% Take transpose for proper x-y orientation
V = V';
[Ex,Ey]=gradient(V);
Ex = -Ex;
Ey = -Ey;
E = sqrt(Ex.^2+Ey.^2);

%% FIGURES:
figure('name','Electric Field Magnitude','color','w','numbertitle','off');
hold on
pcolor(x,y,E)
axis image
shading interp;
colormap jet
xlabel('x [m]','fontsize',14);
ylabel('y [m]','fontsize',14);
title('Electric Field Magnitude');
set(gca,'fontsize',14);
axis equal

figure('name','Electric Field Magnitude','color','w','numbertitle','off');
hold on
pcolor(x,y,E)
axis image
shading interp;
colormap hot
xlabel('x [m]','fontsize',14);
ylabel('y [m]','fontsize',14);
title('Electric Field Magnitude');
set(gca,'fontsize',14);
axis equal

figure('name','Electric Field and Potential distribution','color','w','numbertitle','off');
pcolor(x,y,V)
axis image
shading interp;
colormap jet
hold on,
quiver(x,y,Ex,Ey,2)
colorbar('location','eastoutside','fontsize',14);
xlabel('x [m]','fontsize',14);
ylabel('y [m]','fontsize',14);
title('Electric field and Potential distribution');
set(gca,'fontsize',14);
axis equal