k=9e9; q=1e-9; r0=0.1; u0=kq/r0; [X,Y,Z]=sphere(20); r0=0.1; x=r0X(:)'; y=r0Y(:)'; z=r0Z(:)'; x=[x;zeros(size(x))]; y=[y;zeros(size(y))]; z=[z;zeros(size(z))]; plot3(x,y,z);hold on; u=linspace(1,3,5)u0; [X,Y,Z]=sphere; r=kq./u; Z(X<0&Y<0)=nan; for i=1:5; surf(r(i)*X,r(i)*Y,r(i)*Z); end shading interp; 改写一下这个代码，不改变其功能，并添加注释

% 定义常量
k = 9e9; % 真空中的电磁力常数
q = 1e-9; % 电荷量
r0 = 0.1; % 初始距离

% 计算初始电势
u0 = k*q/r0;

% 生成球体的坐标
[X,Y,Z] = sphere(20);
r0 = 0.1;
x = r0*X(:)'; % x坐标
y = r0*Y(:)'; % y坐标
z = r0*Z(:)'; % z坐标

% 添加底部点坐标
x = [x; zeros(size(x))];
y = [y; zeros(size(y))];
z = [z; zeros(size(z))];

% 绘制底部点和球体
plot3(x, y, z);
hold on;

% 计算不同电势下的球体半径
u = linspace(1, 3, 5)*u0;
[X, Y, Z] = sphere;

% 计算球体半径
r = k*q./u;

% 将球体半径为负值的部分设为 NaN，使其不显示
Z(X<0 & Y<0) = nan;

% 绘制不同电势下的球体
for i = 1:5
surf(r(i)*X, r(i)*Y, r(i)*Z);
end

% 设置光照效果
shading interp;

% 添加注释
title('电势分布');
xlabel('X轴');
ylabel('Y轴');
zlabel('Z轴');
legend('底部点', '球体');
grid on;