% 清除工作区
clear all

% 用户输入
a = input('输入点电荷的位置如[1,0,1;2,0,2]表示位置在(1,0,1),(2,0,2): ');
Q = input('输入点电荷的电荷量，+-表示电性，如[1,-1]: ');
a1 = input('电场线角度间隔: '); % 角度间隔

% 设置绘图范围
xmin = min(a(:,1)) - 4;
xmax = max(a(:,1)) + 4;
ymin = min(a(:,2)) - 4;
ymax = max(a(:,2)) + 4;
zmin = min(a(:,3)) - 4;
zmax = max(a(:,3)) + 4;

% 创建网格
x = linspace(xmin, xmax, 20); % 减少网格密度以提高性能
y = linspace(ymin, ymax, 20);
z = linspace(zmin, zmax, 20);
[X, Y, Z] = meshgrid(x, y, z);

% 计算距离和电势
g = size(a, 1); % 电荷数量
U = zeros(size(X));
for i = 1:g
    R = sqrt((X - a(i,1)).^2 + (Y - a(i,2)).^2 + (Z - a(i,3)).^2);
    U = U + Q(i) ./ R; % 计算电势
end

% 计算电场强度
[Ex, Ey, Ez] = gradient(-U, x(2)-x(1), y(2)-y(1), z(2)-z(1));

% 创建图形
figure;
hold on;
axis equal;
grid on;
view(3);
xlabel('x');
ylabel('y');
zlabel('z');
title('三维点电荷系统的电场线', 'FontSize', 14);

% 绘制电荷
for i = 1:g
    if Q(i) > 0
        plot3(a(i,1), a(i,2), a(i,3), 'r+', 'MarkerSize', 16, 'LineWidth', 2);
    else
        plot3(a(i,1), a(i,2), a(i,3), 'bo', 'MarkerSize', 16, 'LineWidth', 2);
    end
end

% 生成电场线
r0 = 0.2; % 电场线起始半径
theta = deg2rad(a1:a1:180);   % 极角范围 (0° 到 180°)
phi = deg2rad(0:a1:360);      % 方位角范围 (0° 到 360°)

% 只从正电荷发出电场线
for i = 1:g
    if Q(i) > 0
        for j = 1:length(theta)
            for k = 1:length(phi)
                % 球坐标转直角坐标 (修正z坐标公式)
                x_start = a(i,1) + r0 * sin(theta(j)) * cos(phi(k));
                y_start = a(i,2) + r0 * sin(theta(j)) * sin(phi(k));
                z_start = a(i,3) + r0 * cos(theta(j));
                
                % 绘制电场线
                streamline(X, Y, Z, Ex, Ey, Ez, x_start, y_start, z_start);
            end
        end
    end
end

% 设置图形属性
axis([xmin xmax ymin ymax zmin zmax]);
light('Position',[1 1 1]);
lighting gouraud;
material shiny;
hold off;