第1个回答 2008-10-26
clc;
clear;
%对空心球体的空间电位的仿真
%首先绘制两个球面，作为空心球体的模型
%假设内半径a=5m,外半径b=10m
%绘制内球面
a=5;b=10;
x=-5:0.5:5;
y=-5:0.5:5;
[X1,Y1]=meshgrid(x,y);
Z11=sqrt(a^2-X1.^2-Y1.^2);
index=find(X1.^2+Y1.^2>a^2);
Z11(index)=NaN;
plot3(X1,Y1,Z11,'LineStyle','.','Color','g');
hold on;
Z12=-sqrt(a^2-X1.^2-Y1.^2);
index=find(X1.^2+Y1.^2>a^2);
Z12(index)=NaN;
plot3(X1,Y1,Z12,'LineStyle','.','Color','g');
hold on;
%绘制外球面
x=-10:1:10;
y=-10:1:10;
[X2,Y2]=meshgrid(x,y);
Z21=sqrt(b^2-X2.^2-Y2.^2);
index=find(X2.^2+Y2.^2>b^2);
Z21(index)=NaN;
plot3(X2,Y2,Z21,'LineStyle','.','Color','r');
hold on;
index=find(X2.^2+Y2.^2>b^2);
Z22=-sqrt(b^2-X2.^2-Y2.^2);
Z22(index)=NaN;
plot3(X2,Y2,Z22,'LineStyle','.','Color','r');
%求解空间各点的电位，并在空间求导描述
%并且要区分内外空间的数值差异
%r=sqrt(x^2+y^2+z^2);
%在内半径a范围内，空间电位phi=Q/(pi*e_0)*(1/r+1/b-1/a);
%在内半径b之外，空间电位phi=Q/(pi*e_0)*(1/r);
%其中，e0为介电常数
%假设Q=10-10库仑(C);
Q=10-10;
phia=[];
e_0=8.85*1e-12;
r_inita=1;
r_tempa=r_inita;
r_intervala=0.1;
r_interval_consta=r_intervala;
%做内半径a范围内的空间电场分布图，并保存各点的电位数据于phia中
%设定半径步长为1，确保图像清晰，易分辨
while(r_tempa<=5)
[X,Y] = meshgrid(-r_tempa:r_intervala:r_tempa,-r_tempa:r_intervala:r_tempa);
Z=sqrt(r_tempa^2-X.^2-Y.^2);
Z1=-sqrt(r_tempa^2-X.^2-Y.^2);
index=find(X.^2+Y.^2>r_tempa^2);
Z(index)=NaN;
Z1(index)=NaN;
phi_temp1a=Q/(pi*e_0).*(1./(X.^2+Y.^2+Z.^2)+1/b-1/a);
phi_temp2a=Q/(pi*e_0).*(1./(X.^2+Y.^2+Z1.^2)+1/b-1/a);
phi_tempa=cat(1,phi_temp1a,phi_temp2a);
if(r_tempa==1)
phia=phi_tempa;
else
phia=cat(3,phia,phi_tempa);
end
[U,V,W] = surfnorm(X,Y,Z);
quiver3(X,Y,Z,-U,-V,-W,0.2,'k');
[U,V,W] = surfnorm(X,Y,Z1);
quiver3(X,Y,Z1,-U,-V,-W,0.2,'k');
hold on;
r_tempa=r_tempa+1;
r_intervala=r_intervala+r_interval_consta;
end
%做外半径b范围之外的空间电场分布图，并保存各点的电位数据于phib中
%同样设定半径步长为1，范围为10~12
phib=[];
r_initb=10;
r_tempb=r_initb;
r_intervalb=1;
r_interval_constb=0.1;
while(r_tempb<=12)
[X,Y] = meshgrid(-r_tempb:r_intervalb:r_tempb,-r_tempb:r_intervalb:r_tempb);
Z=sqrt(r_tempb^2-X.^2-Y.^2);
Z1=-sqrt(r_tempb^2-X.^2-Y.^2);
index=find(X.^2+Y.^2>r_tempb^2);
Z(index)=NaN;
Z1(index)=NaN;
phi_temp1b=Q/(pi*e_0).*(1./(X.^2+Y.^2+Z.^2));
phi_temp2b=Q/(pi*e_0).*(1./(X.^2+Y.^2+Z1.^2));
phi_tempb=cat(1,phi_temp1b,phi_temp2b);
if(r_tempb==10)
phib=phi_tempb;
else
phib=cat(3,phib,phi_tempb);
end
[U,V,W] = surfnorm(X,Y,Z);
quiver3(X,Y,Z,-U,-V,-W,0.2,'k');
[U,V,W] = surfnorm(X,Y,Z1);
quiver3(X,Y,Z1,-U,-V,-W,0.2,'k');
hold on;
r_tempb=r_tempb+1;
r_intervalb=r_intervalb+r_interval_constb;
end
grid on;
axis equal;
这个程序由于绘图的点比较多，所以运行起来会比较慢，在图像上显示空间的电场分布，所有的电位数据存储在phia,phib中本回答被提问者采纳