%模擬電荷法計算電場程序
function?csm（rLRn）
L=13;r=0.05;E=13000;n=10;d=1;m=10;
%默認值
%下面為求模擬電電荷的坐標
%-------------------
MQ1=[0;L];??%點電荷坐標
for?j=1:n-1
????A（j）=80/81*r;
end
for?j=1:n-1
????B（j）=L-L*j^2/（n）^2;
end
MQ2=[A;B];%環(huán)電荷坐標
MQ=[MQ1?MQ2]%模擬電荷坐標
%-------------------
%下面為匹配點的坐標
%---------------------
K1=[r;L];%對應(yīng)點電荷的匹配點坐標
for?j=1:n-1
????C（j）=r;
end
K2=[C;B];%對應(yīng)環(huán)電荷的匹配點坐標
K=[K1?K2]%匹配點坐標
%--------------------
%下面為校核點的坐標
%-------------------
AP1=[0;L+r];
for?i=1:n-2
????D（i）=r;
end
for?i=1:n-2
????G（i）=（B（i）+B（i+1））/2;
end
AP2=[D;G];
AP=[AP1?AP2]%校核點AP的坐標
%---------------------
%下面為求系數(shù)矩陣P（ij）
%-------------------
pi=3.1415926;
ebcl=1/（36*pi）*10^-9;%空氣介電的常數(shù)，單位為:安秒/伏.米
gyz=1/（4*pi*ebcl）;%公因子1?/?（4*pi*ebcl）
%分離r和z方向
MQ1r=MQ1（1）;
MQ1z=MQ1（2）;
MQr=MQ（1:）;
MQz=MQ（2:）;
Kr=K（1:）;
Kz=K（2:）;
APr=AP（1:）;
APz=AP（2:）;
%下面為求點電荷的矩陣系數(shù)P（in）
for?i=1:n
????PQ1（i）=1/sqrt（（Kr（i）-MQ1r）^2+（Kz（i）-MQ1z）^2）-1/sqrt（（Kr（i）-MQ1r）^2+（Kz（i）+MQ1z）^2）;
end
%下面為求環(huán)電荷的矩陣系數(shù)P
for?i=1:n
????for?j=1:n-1
????????aph1（ij）=sqrt（（Kr（i）+MQr（j+1））^2+（Kz（i）-MQz（j+1））^2）;
????end
end????
for?i=1:n
????for?j=1:n-1
????????aph2（ij）=sqrt（（Kr（i）+MQr（j+1））^2+（Kz（i）+MQz（j+1））^2）;
????end
end
for?i=1:n
?????for?j=1:n-1
????????XD1（ij）=2*sqrt（MQr（j+1）*Kr（i））/aph1（ij）;
????end
end
for?i=1:n
?????for?j=1:n-1
????????XD2（ij）=2*sqrt（MQr（j+1）*Kr（i））/aph2（ij）;
????end
end
for?i=1:n
?????for?j=1:n-1
?????????syms?a
???????XXD1（ij）=int（1/sqrt（1-XD1（ij）*XD1（ij）*sin（a）*sin（a））0pi/2）;
????end
end
for?i=1:n
?????for?j=1:n-1
?????????syms?a
???????XXD2（ij）=int（1/sqrt（1-XD2（ij）*XD2（ij）*sin（a）*sin（a））0pi/2）;
????end
end
for?i=1:n
?????for?j=1:n-1
????????PQ2（ij）=2/pi*（XXD1（ij）/aph1（ij）-XXD2（ij）/aph2（ij））;
????end
end
%系數(shù)矩陣P（ij）
PP=-gyz*[PQ2PQ1‘];
P=double（PP）
%---------------------
%下面為求模擬電荷量
%---------------------
%E為均勻電場場強
for?i=1:n
????V（i）=-E*Kz（i）;
end
Q=inv（P）*V‘
%下面為計算誤差
%---------------
for?i=1:n-1
????APQ1（i）=1/sqrt（（APr（i）-MQ1r）^2+（APz（i）-MQ1z）^2）-1/sqrt（（APr（i）-MQ1r）^2+（APz（i）+MQ1z）^2）;
end
for?i=1:n-1
????for?j=1:n-1
????????Aaph1（ij）=sqrt（（APr（i）+MQr（j+1））^2+（APz（i）-MQz（j+1））^2）;?
????end
end????
for?i=1:n-1
????for?j=1:n-1
????????Aaph2（ij）=sqrt（（APr（i）+MQr（j+1））^2+（APz（i）+MQz（j+1））^2）;
????end
end?
for?i=1:n-1
?????for?j=1:n-1
????????AXD1（ij）=2*sqrt（MQr（j+1）*APr（i））/Aaph1（ij）;
????end
end