模拟电荷法，电位系数

%模拟电荷法计算电场程序
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];%环电荷坐标
MQ=[MQ1 MQ2]%模拟电荷坐标
%-------------------
%下面为匹配点的坐标
%---------------------
K1=[r;L];%对应点电荷的匹配点坐标
for j=1:n-1
    C（j）=r;
end
K2=[C;B];%对应环电荷的匹配点坐标
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的坐标
%---------------------
%下面为求系数矩阵P（ij）
%-------------------
pi=3.1415926;
ebcl=1/（36*pi）*10^-9;%空气介电的常数，单位为:安秒/伏.米
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:）;
%下面为求点电荷的矩阵系数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
%下面为求环电荷的矩阵系数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
%系数矩阵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