模拟电荷法程序

disp('本程序中所有物理量均在国际单位制下');
h=input('输入球板距离h:');
r=input('输入球半径r(ru=input('输入球上电位U:');
y2=h;
y1=h-r*r/h/2;
p1=[0,y1];p2=[0,y2];p3=[0,-y1,];p4=[0,-y2];pa=[0,h-r];pb=[0,h+r];
uu=inline('9*10^9/((a1-a0)^2+(b1-b0)^2)^0.5','a1','a0','b1','b0');
Pa1=uu(p1(1,1),pa(1,1),p1(1,2),pa(1,2));
Pa2=uu(p2(1,1),pa(1,1),p2(1,2),pa(1,2));
Pa3=uu(p3(1,1),pa(1,1),p3(1,2),pa(1,2));
Pa4=uu(p4(1,1),pa(1,1),p4(1,2),pa(1,2));
Pb1=uu(p1(1,1),pb(1,1),p1(1,2),pb(1,2));
Pb2=uu(p2(1,1),pb(1,1),p2(1,2),pb(1,2));
Pb3=uu(p3(1,1),pb(1,1),p3(1,2),pb(1,2));
Pb4=uu(p4(1,1),pb(1,1),p4(1,2),pb(1,2));
P=[Pa1-Pa3,Pa2-Pa4;Pb1-Pb3,Pb2-Pb4];
q=inv(P)*[u;u];
disp('计算结果如下:');
fprintf('q1=%5i ，q2=%5i ，',q);
disp('q3=-q1 ，q4=-q2');
disp('q1,q2,q3,q4的放置位置分别为：');
fprintf('（0，%f），（0，%f），（0，%f），（0，%f)\n',y1,y2,-y1,-y2);
pc=[r,h];pd=[r/2^0.5,h+r/2^0.5];pe=[r/2^0.5,h-r/2^0.5];
Pc1=uu(p1(1,1),pc(1,1),p1(1,2),pc(1,2));
Pc2=uu(p2(1,1),pc(1,1),p2(1,2),pc(1,2));
Pc3=uu(p3(1,1),pc(1,1),p3(1,2),pc(1,2));
Pc4=uu(p4(1,1),pc(1,1),p4(1,2),pc(1,2));
Pd1=uu(p1(1,1),pd(1,1),p1(1,2),pd(1,2));
Pd2=uu(p2(1,1),pd(1,1),p2(1,2),pd(1,2));
Pd3=uu(p3(1,1),pd(1,1),p3(1,2),pd(1,2));
Pd4=uu(p4(1,1),pd(1,1),p4(1,2),pd(1,2));
Pe1=uu(p1(1,1),pe(1,1),p1(1,2),pe(1,2));
Pe2=uu(p2(1,1),pe(1,1),p2(1,2),pe(1,2));
Pe3=uu(p3(1,1),pe(1,1),p3(1,2),pe(1,2));
Pe4=uu(p4(1,1),pe(1,1),p4(1,2),pe(1,2));
P2=[Pc1-Pc3,Pc2-Pc4;Pd1-Pd3,Pd2-Pd4;Pe1-Pe3,Pe2-Pe4];
U2=P2*q;
fprintf('校验值:Uc=%f ，Ud=%f ，Ue=%f\n',U2);
E2=(U2-[u;u;u])/u*100;
fprintf('校验值误差:Ec=%f%% ，Ed=%f%% ，Ee=%f%%\n',E2)
fprintf('球板之间的电容值C=%i\n',(q(1,1)+q(2,1))/u );
disp('求场点P（x,y）的电场强度及电位');
l=inline('((a1-a0)^2+(b1-b0)^2)^0.5','a1','a0','b1','b0');
e=inline('9*10^9/((a1-a0)^2+(b1-b0)^2)','a1','a0','b1','b0');
x=input('输入场点横坐标：');
y=input('输入场点纵坐标：');
if x^2+(y-h)^2fprintf('场点在导体球内,场强E=0,电位U=%d\n',u);
else if y<0
disp('场点在板下,场强E=0,电位U=0');
else
e1=e(x,p1(1,1),y,p1(1,2));
e2=e(x,p2(1,1),y,p2(1,2));
e3=e(x,p3(1,1),y,p3(1,2));
e4=e(x,p4(1,1),y,p4(1,2));
Pp1=uu(p1(1,1),x,p1(1,2),y);
Pp2=uu(p2(1,1),x,p2(1,2),y);
Pp3=uu(p3(1,1),x,p3(1,2),y);
Pp4=uu(p4(1,1),x,p4(1,2),y);
Ex=[e1*(x-p1(1,1))/l(x,p1(1,1),y,p1(1,2))-e3*(x-p3(1,1))/l(x,p3(1,1),y,p3(1,2)),...
e2*(x-p2(1,1))/l(x,p2(1,1),y,p2(1,2))-e4*(x-p4(1,1))/l(x,p4(1,1),y,p4(1,2))]*q;
Ey=[e1*(y-p1(1,2))/l(x,p1(1,1),y,p1(1,2))-e3*(y-p3(1,2))/l(x,p3(1,1),y,p3(1,2)),...
e2*(y-p2(1,2))/l(x,p2(1,1),y,p2(1,2))-e4*(y-p4(1,2))/l(x,p4(1,1),y,p4(1,2))]*q;
Up=[Pp1-Pp3,Pp2-Pp4]*q;
fprintf('电场强度E=(%f)ex+(%f)ey\n电位U=%f\n',Ex,Ey,Up);
fprintf('场强大小|E|=%f\n',l(Ex,0,Ey,0));
end
end